Introduction

HCC is a major health burden, and its incidence is increasing worldwide. With the advancement of medical technology in the fields of medical diagnosis, imaging, and nonsurgical therapies, the practice of early diagnosis and surgical resection of HCC has become increasingly prevalent. However, it is disconcerting to note that a significant proportion of patients, exceeding 70%, are already in the advanced stage at their initial diagnosis, rendering them ineligible for surgical resection1. Nevertheless, the 5-year survival rate is lower than 70% in patients with early HCC who have received surgical resection. Despite the constrained accessibility of diagnostic tools such as histology and imaging, there exists a substantial demand for expeditious, precise, and convenient diagnostic approaches for HCC2. Consequently, it is imperative to advance novel diagnostic methods for monitoring cancer progression. While clinically pertinent biomarkers such as alpha-Fetoprotein (AFP) are frequently employed in HCC diagnosis, their application is surrounded with numerous controversies and limitations3,4.

Programmed cell death (PCD) is a meticulously orchestrated and systematic mechanism of intracellular self-degradation, which plays a vital role in preserving tissue integrity and governing the progression of pathological conditions. PCD assumes a significant role in a diverse range of diseases, encompassing tumors, neurodegenerative disorders, and cardiovascular ailments5. Dysregulation of PCD is associated with various diseases, including cancer. In the case of HCC, the most common type of liver cancer, dysregulation of PCD has been implicated in the initiation, progression, and treatment of the disease6,7. Extensive investigation has revealed that distinct characteristics exhibited by various PCD patterns hold significant implications for disease progression and therapeutic interventions. Apoptosis was once synonymous with PCD but now many other forms have been described, cell death subroutines can be categorized based on their discernible morphological and biochemical characteristics, encompassing apoptosis, autophagy, anoikis, lysosome-dependent cell death, immunogenic cell death, necroptosis, ferroptosis, NETosis, pyroptosis, disulfidptosis, entotic cell death, cuproptosis, parthanatos, netotic cell death, alkaliptosis, and oxeiptosis, each of which exhibits distinct molecular cascades and regulatory pathways8,9,10. Studies have shown that distinct PCD patterns are not mutually exclusive but share a coordinated system. The redundancy ensures that the process of cell death can be triggered before they cause unnecessary cell death even if one of the pathways is impaired or not functioning properly. The gradual emergence of supportive evidence has elucidated the principal molecular mechanisms underlying each subroutine, consequently presenting a range of potential targets for cancer therapy. However, the intricate interconnections among distinct cell death subroutines still require further clarification. The utilization of bioinformatics techniques, including transcriptomics, proteomics, and bio-network analysis, facilitates the elucidation of molecular mechanisms underlying PCD, identification of potential targets, and refinement of therapeutic approaches. In HCC, dysregulation of PCD pathways can occur through various genetic and epigenetic mechanisms, each involving different factors including alterations of gene expression involved in PCD, dysregulation of mitochondrial function, and changes in the tumor microenvironment (TME). For example, increased levels of anti-apoptotic proteins such as the BCL-2 family members (for example, Bcl-2, Bcl-XL, and Mcl-1) and decreased expression of pro-apoptotic proteins (BAX, BAK) have been described in HCC, leading to a reduction in PCD. Additionally, alterations in the expression of microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) have been shown to modulate PCD biological process in HCC cells6,7,11. In this context, the discovery of comprehensive PCD patterns will promote the potential benefits for patients through the utilization of small-molecule compounds that target classical PCD modalities. Notable advances have been accomplished in targeted cancer therapy for the past several decades, several targeted small-molecule drugs have been approved for the treatment of various types of tumor. Emerging evidence has elucidated that local immune suppression of the TME exhibits a strong association with PCD patterns in cancer growth, metastasis and even tumor immune escape12,13. This correlation encompasses the stimulation of inflammatory reactions, facilitation of cell phagocytosis, mediation of inflammatory responses, recruitment and activation of immune cells such as natural killer (NK) cells, dendritic cells (DCs), and CD8 + T cells, as well as the proliferation of cytokines14. The dysregulation of PCD in HCC has important implications for the development of novel therapeutic strategies. Targeting PCD pathways, such as apoptosis, autophagy, and necroptosis, has been explored as a potential therapeutic strategy for HCC. For example, several agents that induce apoptosis, such as BH3 mimetics and Smac mimetics, have shown promise in preclinical studies15. Several studies have developed a variety of nanoparticle drugs, which can initiate PCD progress and reprogram the immunosuppressive TME. When combined with an immune checkpoint blockade (αPD-1, αPD-L1), they can reduce tumor growth and enhance gemcitabine chemotherapy efficacy for cancer treatment16,17. Although immune- and targeted therapeutics have been approved for HCC treatment, biomarker-based stratification of patients for optimal response to therapy is an unmet need.

While increasing advances have been employed in understanding different PCD patterns in several cancer types, a comprehensive evaluation and integrated multi-omics analysis of different cell death subroutines are still essential. Insight into the molecular players and pathways involved in cell death are important for understanding the mechanisms of HCC. Moreover, the association between PCD patterns and HCC prognosis remains unclear, the dearth of research highlights the need for further investigation in this area. The integration of different modalities is necessary to elucidate potential underlying mechanisms, which were contingent upon coordinated changes occurring across different regulatory molecular layers. In our study, we defined a novel PCD index (PCDI) and developed a novel ML-based classifier for HCC patients. We then sought to validate the prognostic value of this PCDI across three independent cohorts. In addition, we performed a scheme of integrative analyses of multi-omics data, including TME and gene mutation landscape, single-cell RNA-seq, and small targeted therapeutic molecules sensitivity analysis, to gain a greater in-depth cellular and molecular understanding of PCD, the heterogeneous nature of HCC will become more apparent. Finally, we validated the expression level of DLAT through human tissues including HCC tissue and adjacent hepatic tissues in a real-world cohort. These results provide new insights informing the patient stratification strategies for targeted therapeutic approaches for HCC. Overall, the dysregulation of PCD in HCC is a complex process that involves multiple mechanisms. Further exploration is still needed to better investigate the PCD process in the development and progression of HCC and to develop more effective targeted therapies for this devastating disease.

Materials and methods

Study subjects and biospecimens of HCC patients

A schematic illustration of our research design was displayed in Fig. 1. 983 individuals with HCC were enrolled from three public RNA-seq cohorts (TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; ICGC, International Cancer Genome Consortium) and a proteome cohort (PDC000198)18. The detailed demographic characteristics and clinical information of the study subjects were obtained from corresponding raw records. Samples without complete overall survival or recurrence-free survival data were removed. Single‑cell RNA‑seq data were downloaded from GSE14961419. Detailed demographic and clinical characteristics of participants at baseline were presented in Supplementary Table S1.

Fig. 1
figure 1

The flowchart for comprehensive analysis of 16 PCD patterns in HCC.

Full size image

Data processing

For bulk RNA-seq data, the HTSeq generated read counts were transformed into TPM (per kilobase million) values by dividing raw counts by gene length, then scaling with a factor of 1e6and dividing by total read counts in the sample. For proteome analysis, data normalization median centering was performed across the total proteins for each sample to correct sample loading differences using the NormalyzerDE package20. In normalized samples, these reporter intensities should have a log TMT ratio centered at zero. For single-cell RNA-seq analysis, the UMI counts of each cell were normalized using the NormalizeData () function and identified highly variable features for each sample using the FindVariableFeatures function with 10,000 scale factors in R package Seurat21. We applied dimensionality reduction and unsupervised cell clustering to identify each cell subpopulation, including principal component analysis (PCA) and t-stochastic neighbor embedding (t-SNE), after batch effect correction in R package Harmony22.

Construction and validation of a prognostic PCD index

Several machine learning algorithms, including random forest (RF), XGBoost, Boruta, and Elastic net proportional hazards models, were employed to identify and integrate gene signatures linked to the overall survival (OS) of HCC patients23,24. The model efficacy was evaluated through cross-validation and compared against other clinical and genomic predictors. We used a ten-fold nested cross-validation strategy to assess the robustness and generalizability of the PCD index (PCDI) in three independent HCC cohorts. The performance of the model was assessed using various metrics such as receiver operating characteristic (ROC) curves, area under the curve (AUC) scores, and Kaplan–Meier survival curves.

Estimation of immune infiltration and matrix composition

In order to estimate the profiling of immune cells in the TME and abundances of immune cells from the bulk RNA-seq, we employed CIBERSORT to calculate the proportions of 22 human leukocyte cell subsets defined in the CIBERSORT package for each bulk RNA-seq sample25.

Chemotherapy sensitivity prediction

Drug responses of cancer cell lines were represented by the half-maximal inhibitory concentration (IC50), the slope of the dose response curve, and area under the dose response curve (AUC), and all screened cell line/ drug combinations were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC)26. We employed a ML model to predict the clinical response to drug response in TCGA patients by training a multinomial logistic regression classifier on IC50 values on 982 cell lines from the GDSC, in the R package pRRophetic27.

Molecular subtyping of PCD signatures by NMF

The candidate signatures with high median absolute deviation (MAD) value (MAD ≥ 0.5) across all samples were adopted for unsupervised decomposition and clustering analysis. Based on the maximum cophenetic correlation score and the average silhouette width of the consensus membership matrix (non-negative matrix factorization (NMF) method), the optimal number of clusters (optimal factorization rank k, 2 to 10) was determined. The NMF clustering was performed on 100 runs based on the 2000 most variable genes by the R package NMF, as determined in the RNA-seq analysis. To conduct a comprehensive prognostic assessment of three NMF clusters, KM survival curves were generated, with estimated survival probabilities computed using a product limit formula and log-rank tests.

Bioinformatic analysis

Unless otherwise noted, most of the data processing pipelines described in our study (data cleaning, feature selection, significance testing, correlation evaluation, supervised or unsupervised clustering, and plots producting) were performed as previously described in R (https://www.r-project.org/). Differentially expressed genes (DEGs) were identified using the moderated (empirical Bayesian) t-test implemented in the R package limma. The genes with an absolute log2 fold change > 1 and adjusted P value ≤ 0.05 were identified as DEGs. Gene set variation analysis (GSVA) or gene set enrichment analysis (GSEA) was used to identify the most significant pathways using the Molecular Signatures Database (MSigDB) including Kyoto encyclopedia of genes and genomes (KEGG), Reactome, WikiPathways and other customized gene sets that have been reported between different subgroups18,28,29. The maftools package in R was used to summarize, analyze, annotate, and visualize the variants studied. Significant differences in independent samples were assessed using Fisher exact test for categorical variables, and t-test or Mann–Whitney test for continuous variables as appropriate. All P values were two-sided and P-value < 0.05 was considered statistically significant.

Sample collection and quantitative PCR

A total of five fresh HCC tissue samples and matched normal para-cancerous tissues were obtained from HCC patients who have undergone surgical resection at the First Affiliated of Bengbu Medical University (Bengbu, China). Participants or their legal guardians signed a written informed consent agreeing on the use of tumor tissues for research, approved by the Bengbu Medical University Ethics Committee, including any relevant details. RNA extraction, cDNA synthesis and real-time quantitative PCR were performed as previously described18. All PCR amplification procedures were performed in triplicates using quantitative PCR kit (Takara) and ABI7500 Real-Time PCR System (Applied Biosystem Company). Sequence-specific primers for DLAT and GAPDH were as follows: DLAT forward primer (5′-AGGAACTTCGGCATTTCATCG-3′) and DLAT reverse primer (5′-TGTCCCGGTATTGTAGTCCCA-3′); GAPDH forward primer (5′-ACAACTTTGGTATCGTGGAAGG-3′), and GAPDH reverse primer (5′-GCCATCACGCCACAGTTTC-3′).

Western blotting and immunohistochemical staining (IHC)

In clinically matched HCC tissues, DLAT protein levels were measured by Western blot and immunohistochemical (IHC) staining. Proteins were extracted in RIPA buffer supplemented with complete, EDTA-free protease inhibitor cocktail (Beyotime Biotechnology, China) and phosphatase inhibitor cocktail 3 (Sigma). Protein extracts were separated on sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE, 8–12% gels) and transferred onto 0.22 μm PVDF membranes (Merck Millipore Ltd., Massachusetts, USA). Membranes were blocked with 5% milk in TBST (tris-buffered saline with 1% Triton X-100, 10 mM Tris–HCl pH 8.0, 150 mM NaCl, and 0.1% Tween-20), or 5% bovine serum albumin (BSA) in TBST, for 30 min at room temperature. Then the membranes were washed with TBST, and incubated with primary antibodies (anti-DLAT, 68,303–1-Ig, Proteintech; anti-GAPDH, AP0066, Bioworld) in 5% BSA in TBST overnight at 4°℃, then probed with HRP-conjugated secondary antibodies (goat anti-mouse IgG-HRP, sc-2005; goat anti-rabbit IgG-HRP, sc-2004, Santa Cruz Biotechnology). For IHC analysis, slides were incubated with primary antibodies for 1 h at room temperature and washed with TRIS–HCL buffer after antigen retrieval, then sections were incubated with HRP-conjugated secondary antibodies for 15 min followed by visualizing with a detection kit (DAB). Detailed methods for WB and IHC also refer to the reported work. Unprocessed images of gels and membranes are presented in our Supplementary Materials. All methods in our study were carried out in accordance with relevant guidelines and regulations.

Results

Identification of PCD-associated DEGs between HCC and non‑malignant liver samples

The overarching purpose of the literature survey was to integrate previously reported signatures that have been implicated in PCD-associated biological processes. A total of 16 PCD patterns yielded 1599 nonredundant gene signatures were identified. Specifically, 580 apoptosis genes, 367 autophagy genes, 338 anoikis genes, 220 lysosome-dependent cell death genes, 117 immunogenic cell death genes, 101 necroptosis genes, 88 ferroptosis genes, 69 NETosis genes, 52 pyroptosis genes, 16 disulfidptosis genes, 15 entotic cell death genes, 14 cuproptosis genes, 9 parthanatos genes, 8 netotic cell death genes, 7 alkaliptosis genes, and 5 oxeiptosis genes (Supplementary Table S2). We next checked the expression patterns of PCD-associated genes between HCC tumor tissues and adjacent normal liver tissues. This analysis identified 950 and 69 genes that were significantly up- or down-regulated, respectively (log2FC > 0.5, false discovery rate (FDR) < 0.05). RNA-seq data in TCGA was analyzed and represented as heatmap (Fig. 2A) and volcano plot (Fig. 2B). Furthermore, three independent cohorts also demonstrated that PCD-related signatures were dysregulated in HCC (Supplementary Fig. S1). T-SNE analysis for the expression of 1109 DEGs showed that two groups were separated, indicating that they were well distinguished based on the expression of 1109 DEGs (Fig. 2C). Enrichment analysis revealed several significantly enriched GO terms were also related to programmed cell death, including regulation of autophagy, intrinsic or extrinsic apoptotic signaling pathways (Fig. 2D). KEGG pathway enrichment analysis shows the PCD gene signatures to be significantly enriched in necroptosis, ferroptosis, autophagy, apoptosis, mTOR/ P53/ NF-κB and metabolic-associated pathways (Fig. 2E). Boxplot revealed that 16 PCD patterns differed between HCC and normal liver tissues (Fig. 2F). Furthermore, Pearson correlation analysis revealed a tightly correlation between several different PCD patterns and celluar metabolism. Moreover, the PCD patterns including ferroptosis and autophagy were associated with glycolysis and lipid metabolism (Fig. 2G). UpSet plot showed the intersection analysis among different PCD patterns, which suggested that they also share some target signatures that were involved in the PCD-associated biological process (Fig. 2H).

Fig. 2
figure 2

Distribution characteristics of different PCD patterns in HCC. (A,B) Heatmap and volcano plot of the PCD-related DEGs between HCC and normal liver tissues. (C) The t-SNE analysis of two sample groups using 1109 DEGs. Such normal and HCC tissues formed their own clusters, and appeared distant from each other on t-SNE. (D,E) GO enrichment analysis of molecular functions (MF), biological processes (BP), and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms. (F) Comparison of the relative levels of PCD patterns between two sample groups based on GSVA analysis. (G) Heatmap for the correlations between 16 PCD patterns and the 20 metabolism-related pathways by pearson’s analysis.

Full size image

Independent prognostic value of PCD-associated signatures

For each cohort, to evaluate associations of signature expression with overall survival, Cox proportional hazard regression models were fit upon the DEGs, using OS status as response, comparing high- to low-signature expression. 64 signatures were then combined with the intersection of 3 HCC cohorts through machine learning analysis to identify a series of gene signatures that can predict a patient to experience mortality. Using four feature selection algorithms (Elastic Net, RF, XGBoost, and Boruta) (Supplementary Fig. S2A–E), and a representative intersection (Supplementary Fig. S2F), we further reduce the number of signatures based on univariate Cox regression. Based on the signatures determined by each feature selection algorithm, we built multivariate Cox models. An intersection of 7 signatures remained with individual coefficients, HTRA2, SLC2A1, FTL, G6PD, and DLAT were identified as independent risk factors and integrated to build a PCD-related prognostic model based on step-wise multivariate Cox regression analysis. Forest plots showed the OS hazard ratios of cases in the TCGA cohort using a univariate Cox model (Supplementary Fig. S2G). Time-dependent ROC analysis confirmed that Elastic Net models outperformed the prediction of OS compared with other models at various time points (Fig. 3A,B) (P < 0.05). The intersection model achieved considerable performance with significantly fewer signatures, demonstrating the importance of the robustness of signatures confirmed by multiple ML algorithms.

Fig. 3
figure 3

Construction of 5-gene signature prognosis model (PCDI) for HCC risk stratification. (A,B) Performance of five feature selection strategies in prognostic models at different time points. An area greater than 0.75 under the ROC curve was considered good performance (C) Restricted cubic splines (RCS) determined the best fitting relationship and cutoff value between the risk score and hazard ratio of patient mortality. (D) Kaplan–Meier analysis demonstrated that patients with higher PCDI exhibited worse overall survival in the TCGA cohort. (E) Distribution of patients in different risk groups. (F,G) The two risk groups exhibited absolute dissimilarity in the PCA and t-SNE analysis (H,I) GO and KEGG enrichment analysis between high- and low-risk groups.

Full size image

Risk score (programmed cell death index, PCDI) was defined as the following formula: 0.021*SLC2A1 + 0.217*DLAT + 0.032*G6PD + 0.041*HTRA2 + 0.016*FTL. We further investigated the linearity assumption of the association between risk score and OS using restricted cubic splines with three knots and determined an optimal cutoff value (risk score = 0.84, P < 0.001) with a hazard ratio (HR) of 1 (Fig. 3C). Patients in each cohort were divided into a high- and low-risk group according to the redefined cutoff value. Kaplan–Meier (KM) survival curves showed that higher is associated with worse survival (HR = 2.42, P < 0.001) (Fig. 3D,E). PCA and t-SNE analyses identified two well-defined subgroups based on the risk model (Fig. 3F,G). For a more comprehensive assessment of the two risk groups, we performed GO and KEGG analyses where a series of pathways associated with cell proliferation were activated in the high-risk group, including chromosome segregation, mitotic nuclear division, cell cycle, ECM-receptor interaction, and hippo signaling pathway, etc. (Fig. 3H,I).

Validation of the PCDI in independent HCC cohorts

In three independent cohorts, patients were divided into two risk groups based on cutoff values in the TCGA cohort (Fig. 4A). Consistently, the groups with higher scores had worse outcomes compared to the group with lower risk scores (Fig. 4B, GEO cohort, HR = 1.80, P = 0.013; ICGC cohort, HR = 4.05, P < 0.001; Proteome cohort, HR = 3.81, P < 0.01, respectively). Figure 4C,D indicated the PCDI model had better discriminatory ability among high- and low-risk groups, based on dimensionality reduction analysis with both PCA and t-SNE. Importantly, the time-dependent AUC values demonstrated the robustness of the PCDI in three independent cohorts (Fig. 4E).

Fig. 4
figure 4

External validation for PCDI model in three independent HCC cohorts. Risk distribution, and survival scatter plot (A), Kaplan–Meier curve (B), PCA (C) and t-SNE (D) analysis, Time-dependent ROC curve analysis (E) in three validation cohorts.

Full size image

Inter-tumor TME heterogeneity in two risk groups

Next, we analyzed the heterogeneity of immune infiltration in the two subtypes. Inthe TCGA cohort, the low-risk group had more infiltrating CD8 + T cells (P < 0.05), while the high-risk group had more infiltrating M2 macrophages and Tregs (P < 0.05) (Fig. 5A). This finding remained significant within the ICGC cohort while diverging from the results observed in the GEO cohort. CD8 + T cells and Tregs were not significant in the GEO cohort (P > 0.05). Immune cell deconvolution analyses revealed that the risk score was positively correlated with M2 macrophages and Tregs infiltration in three HCC cohorts (Fig. 5B). Furthermore, the patients with higher M2 macrophage infiltration had shorter OS and RFS than those with lower infiltration (Fig. 5C). The findings strongly endorse the potential significance of M2 macrophages in HCC. In addition, our further analysis showed that higher CD8 + T levels predicted a more optimistic prognosis (TCGA, OS: P < 0.001, RFS: P = 0.022; ICGC, OS: P = 0.224, RFS: P = 0.359; ICGC, OS: P = 0.154) (Supplementary Fig. S3A). However, high Tregs often indicated a shorter OS or RFS (Supplementary Fig. S3B), but statistical differences were found only in the ICGC cohort.

Fig. 5
figure 5

Heterogeneity of immune infiltration between two risk groups. (A) Violin plot visualizing the abundance of 22 immune cells between high and low-risk groups in 3 HCC cohorts. (B) Analysis of the correlation between PCDI and infiltration levels of M2 macrophages, Tregs and CD8 + T cells (C) Higher M2 macrophages infiltration weas associated to worse prognosis in HCC.

Full size image

Unsupervised clustering redefines the three molecular subtypes associated with PCD

Unsupervised NMF consensus clustering was used to reveal a novel molecular classification of HCC patients based on the PCD-associated genes with median absolute deviation (MAD) value > 1. Next, k = 3 was regarded as the optimal number of clusters according to the consensus maps and silhouette width value (Fig. 6A,B). The t-SNE plot and transcriptome heatmap reveal three distinct blocks of genes whose expression patterns vary as a group across genotypes (Fig. 6C,D).The KM analysis showed a significant prognostic difference between three subgroups, with longer OS (P < 0.001) and RFS (P = 0.013) in PCDcluster 1 than another two subgroups (Fig. 6E). To determine the robustness of the molecular subtype, we performed NMF in three additional independent cohorts. The robustness of clustering assignment was further confirmed by consensus matrix heatmap, rank survey, and silhouette width analysis that in three independent cohorts (Supplementary Fig. S4, Supplementary Fig. S5A). The t-SNE plot generated by PCD-associated genes demonstrated the three HCC subtypes clearly separated from each other (Supplementary Fig. S5B). Further analysis revealed the consistent presence of a HCC subtype exhibiting shorter OS and RFS across four distinct cohorts (Supplementary Fig. S5C), irrespective of statistical significance.

Fig. 6
figure 6

Unsupervised hierarchical clustering based on genes related to PCD. (A-B) HCC samples were clustered in three robust subtypes by NMF algorithm in the TCGA cohort using1000 iterations when the optimal cluster number k = 3. The homogeneity of red-coloring and average silhouette seen in the graphs indicate the presence of 3 clusters of HCC patients. (C) Visualization of cluster results using t-SNE analysis. The three clusters were separated from each other (D) Heatmaps of differential DEGs amomg three subtypes. (E) Survival analysis revealed prognostic differences (right, OS; left, RFS) between three HCC subtypes.

Full size image

PCD clusters-specific clinicopathological characteristics of HCC

To reveal relevant clinical characteristics, we analyzed the distribution of clinical characteristics between the PCD clusters in the four HCC cohorts due to their clinical data integrity (Fig. 7). The results of the chi-square test revealed several noteworthy associations between clinicopathological characteristics and HCC subgroups. Specifically, PCD cluster 3 has a higher risk score and proportion of high-risk patients, in four HCC cohorts (P < 0.001). There was a significant correlation observed between PCD cluster 1 and lower pathological tumor grade (TCGA, P = 0.002; GEO, P = 0.022; ICGC, P = 0.002), TNM stage (TCGA, P = 0.027; GEO, P = 0.022; ICGC, P = 0.022), and AFP level (TCGA, P < 0.001; GEO, P < 0.001; Proteome, P < 0.001). Additionally, we identified the risk score was gradually increased corresponding to the tumor stage or pathological grade (Fig. 8A,B). Metabolic pathways play a crucial role in the regulation of PCD. Various metabolic pathways, such as glucose metabolism, lipid metabolism, and iron metabolism, can influence different types of PCD processes, including apoptosis, autophagy, ferroptosis, and necrosis. GSVA revealed elevated metabolic activity in pathways related to glucose, lipids, amino acids, and other metabolism-related processes in PCDcluster 1, as opposed to nucleotide metabolism. Additionally, PCDcluster 1 demonstrated an association with several HCC-associated subtypes that are correlated with a more favorable prognosis.. The oncogenic pathways associated with HCC progression, including hypoxia, TGF-β, Wnt and ferroptosis, were enriched in PCD cluster 2 and cluster 3 (Fig. 8C). Notably, the same conclusion applied to another three independent HCC cohorts (Supplementary Figs. S6, S7). In summary, metabolic pathways were important factors in regulating PCD.

Fig. 7
figure 7

Heatmap of differences in clinical and histological characteristics among subtypes and risk groups in four HCC cohorts. We compared categorical variables using the Chi-square (χ2) test for dichotomous variables and continuous variables using the t test.

Full size image
Fig. 8
figure 8

(A) Relative abundance of risk score between different histological grade and pathological stage (ANOVA test). The line and box represent median and upper and lower quartiles, respectively. (B) The levels of risk score were compared across different pathological stages in GEO (left), ICGC (middle), and Proteome cohort (right). (C) Heatmap of enriched pathways in each subtype using GSVA based on metabolism- and HCC-associated gene sets.

Full size image

Mutant landscape between two risk group

Previous studies have correlated intrinsic genomic and epigenomic alterations with regulating PCD biological process and metabolic features of tumor cells. As shown in Supplementary Fig. S8A, among the TCGA and Proteome cohorts, the top four significantly mutated genes were identified, including TP53 (28–59%), TTN (25–37%), CTNNB1 (19–24%), and MUC16 (16–23%). Consistently, TP53 mutation frequencies were relatively higher in the high-risk group than the low-risk. The KM analysis showed a significantly favorable long-term prognosis in TP53 wild-type group rather than in TP53-mutant group, in two HCC cohorts (TCGA, OS: P = 0.012, RFS: P = 0.011; Proteome, OS: P = 0.014), which was consistent with the conclusion of previous study (Supplementary Fig. S8B). We performed GSEA (Supplementary Fig. S8C) and GSVA (Supplementary Fig. S8D) analyses between subgroups, and the mutant group contained higher glycolytic activity while the wild-type group having higher lipid metabolic activity. Previous studies have shown that loss of p53 transcriptional activity by mutation (missense mutations) or decreased expression accelerates glycolysis and leads to glutamine and lipid accumulation in tumor cells, which is metabolized through TCA cycle to provide fuel for cancer progression. Taken together, these results strongly support a functional role for TP53 mutations in PCD.

Novel targeted therapies for HCC subclasses

The sensitivity prediction score using the PCDI model significantly correlated with several small molecule drug response in GDSC project. To evaluate the prediction performances depending on the characteristics of the cell lines and drugs, the drug IC50 values were used for drug sensitivity measurements between high-risk and low-risk groups in TCGA cohort. Supplementary Fig. S9 showed that the high-risk group was more sensitive and significantly correlated to commonly used kinase inhibitors such as Alisertib (R =  − 0.361, P < 0.05), AZD7762 (R =  − 0.32, P < 0.05), PD173074 (R =  − 0.363, P < 0.05) and Wee1 inhibitor (R =  − 0.489, P < 0.05), and molecular targeting agents such as Axitinib (R =  − 0.460, P < 0.05) and Nilotinib (R =  − 0.469, P < 0.05), rather than Doramapimod (R = 0.389, P < 0.05). In addition, the high-risk group was more sensitive to Docetaxel, a microtubule degregization inhibitor (R = -0.311, P < 0.05).

Comprehensive characteristics of programmed cell death at the single-cell level

After the implementation of stringent quality control and batch correction, 71,566 single cells were retained for the following analyses (Fig. 9A,B). Our scRNA-seq analysis identified 37 cell clusters using t-SNE visualization (Fig. 9C). Based on the expression of cell type-specific marker signatures, these 37 clusters were then classified into 8 distinct cell populations and these cell populations were present in each patient (Fig. 9D,E). The total counts of CD8 + T cells were decreased in the primary tumor (PT), portal vein tumor thrombus (PVTT), and metastatic lymph node (MLN) tissues compared with the non-tumor liver (NTL) tissues, while NTL tissues had lower infiltration of CD4 + T cells. Besides, NTL tissues showed significant enrichment of CD8 + T cells as compared to MLN and PVTT tissues (Fig. 9F,G). The functions and pathways of various major cell types were meticulously characterized using scGSVA, while the activity of pathways was evaluated through AUCell scoring. The assessment of PCD pathway activity across all cell populations was conducted by scoring all cell types within the comprehensive model. Given the interspersed presence of distinct PCD patterns throughout multiple different cell types, there was a high degree of cell-type-specific difference despite the latitudinal gradient (Fig. 10A). The results showed that cuproptosis, parthanatos, ferroptosis, oxeiptosis had higher distribution in malignant hepatocyte cells. Anoikis, and disulfidptosis had higher pathway activity in specific cell types, especially in endothelial, and fibroblast. Myeloid showed high pathway activity in pyroptosis, NETosis, Lysosome-dependent cell death, and immunogenic cell death (Fig. 10B,C).

Fig. 9
figure 9

scRNA-seq transcriptomic landscape of multicellular ecosystem in primary HCC and normal liver tissues. (A) Single cell transcriptome batch correction for different patient using the Harmony algorithm. (B) The distribution of specific cell populations across patients, samples, tissues, sites, viral infection status, and tumor stage. (C) The t-SNE plot of 37 cell clusters from the multicellular ecosystem of HCC patients. (D) Dot plot showing the canonical marker signatures of 8 major cell types in 37 cell clusters. (E) t-SNE plot of all cells with cell-type annotations. (F,G) t-SNE plot and stacked barplots showing the percentages of major cell types in each tumor site.

Full size image
Fig. 10
figure 10

Heterogeneity of PCD patterns at the single-cell level. (A) Heatmap of the average score of 16 PCD patterns in different cell subtypes based on scGSVA analysis. (B,C) T-SNE and violin plot of the relative PCD score between different cell types.

Full size image

Construction and validation of integrated Nomogram model for optimal risk stratification and survival prediction in HCC

To further translate PCDI in the clinic, a clinical prediction score was created to achieve further performance boost based on a nomogram model. Supplementary Fig. S10 listed the clinical and pathologic factors that were significant on univariate and multivariate stepwise Cox proportional hazards regression analysis. Higher PCDI, therefore, should be regarded as an independent risk factor for mortality that should be integrated into comprehensive clinical decision making in four HCC cohorts (TCGA cohort: HR: 1.270, 95%CI: 1.138–1.418, P < 0.001; GEO cohort: HR: 1.586, 95%CI: 1.171–2.148, P = 0.003; ICGC cohort: HR: 1.477, 95%CI: 1.135–1.921, P = 0.004; Proteome cohort: HR: 1.543, 95%CI: 1.301–1.830, P < 0.001). According to the results, we developed a nomogram based on PCDI (risk score), stage, and virus infection status in the TCGA cohort, corrected by age and gender (Fig. 11A). Subsequently, we performed a comparative analysis of the AUC values for the nomogram, PCDI, integrated clinical model, clinical and pathological characteristics at various time intervals to assess the predictive accuracy of different models. The nomogram achieved desirable predictive accuracies than other models (Fig. 11B). Additionally, a comparison was made between the overall survival probabilities of individuals with high and low nomogram scores. As depicted in Fig. 11C, it was observed that patients with high score experienced a significantly shorter OS compared to their counterparts in the TCGA cohort. The calibration curve demonstrated good agreement between the predictions made by the nomogram and the actual outcomes in the TCGA cohort (Fig. 11D). The graphical representation provided by DCA effectively demonstrated that Nomogram yielded greater net benefits in terms of overall survival compared to other parameters at four specific time intervals (Fig. 11E). The Nomogram performed well on internal and independent external validation cohorts further highlighted its potential usefulness (Fig. 11F).

Fig. 11
figure 11

Evaluation and validation of nomogram prediction model. (A) A hierarchical nomogram was developed to predict the prognosis of HCC patients. (B) The time-dependent AUC curves revealed that the Nomogram performed better than the other models. The level of the dotted line is 0.70 (C) KM analysis revealed prognostic differences in different subgroups based on nomogram score. (D) The calibration curve indicated the predictive probability of the nomogram at each time node. The calibration of the Nomogram predicted probabilities (solid line) were in close proximity compared to a perfectly calibrated model along the diagonal of the plot (dashed line). (E) Decision curve analysis (DCA) of the Nomogram predicts the net benefit in the mortality risk of HCC at 1, 2, 3, and 5 years. (F) The time-dependent AUC confirmed the robustness of the Nomogram in three externally independent validation cohorts.

Full size image

Verification of expression levels of key gene signatures related to PCDI

Among the five PCDI signatures incorporated in the risk model, HTRA2, SLC2A1, FTL, G6PD, and DLAT presented a continuous elevation in the high-risk group, compared to the low-risk group and normal liver tissue, which were consistent in three RNA-seq cohorts (Supplementary Fig. S11A). KM analysis also showed a significantly shorter OS (Supplementary Fig. S11B) and RFS (Supplementary Fig. S11C) of HCC patients with high expression of these four signatures (HTRA2, SLC2A1, G6PD, and DLAT) compared to the low expression level. It is well recognized that DNA methylation is one of the important epigenetic regulation mechanisms for controlling gene expression. We future analyzed the methylation status of a series of promoters corresponding to five signatures in different HCC and paracancer tissues. Our analysis revealed distinct differences in the levelof differential methylation observed for HCC and normal liver tissues (Supplementary Fig. S12A-B). Regression analysis was performed to test the association of gene expression with methylation levels in HTRA2, SLC2A1, FTL, G6PD, and DLAT genes (Supplementary Fig. S12C). These findings underscore the critical roles that specific methylation states play in oncogenesis. Further, we analyzed RNA expression of five gene signatures in 81 HCC cell lines including 22 in the CCLE dataset30. The heatmap illustrated the expression patterns in 81 HCC cell lines for five gene signatures (Supplementary Fig. S12D). Gene expression patterns that are specific to certain cell lines can help us to better understand multiple biological contexts in future studies. DLAT, also known as dihydrolipoamide s-acetyltransferase, played a crucial role in encoding component E2 of the pyruvate dehydrogenase complex (PDC), a multi-enzyme system involved in the metabolic pathway responsible for converting pyruvate to acetyl-CoA. Extensive research has demonstrated its oncogenic properties and association with copper metabolism disorders and cuproptosis, thus establishing DLAT as a novel gene signature in these contexts. Indeed, qRT-PCR and western blot showed consistent results that DLAT mRNA (Fig. 12A) and protein (Fig. 12B) levels were significantly increased in a panel of paired HCC tissues compared to those in normal liver tissues. An additional analysis was performed to verify the protein expression levels of DLAT available from the HPA (Human Protein Atlas) database. The findings of this study validated the cumulative expression of DLAT at the genetic levee within HCC tissues, in comparison to normal tissues (Fig. 12C). We next performed IHC analysis on DLAT protein levels and observed strong immunostaining intensity located in cytoplasmic and membranous in HCC tissues (Fig. 12D). These observations highlighted a potential role for DLAT in HCC pathogenesis. Furthermore, subcellular localization analysis indicated that the DLAT protein predominantly localizes to mitochondria, as evidenced by green immunofluorescence from the HPA database, indicating an important role in celluar metabolism (Supplementary Fig. S13A). Specifically, we have demonstrated an association between metabolism and PCD in previous studies, and further single-gene GSEA analysis revealed that high levels of DLAT, a key regulator in cuproptosis, are also associated with activation of other PCD patterns, including autophagy, ferroptosis, and lysosome-dependent cell death (LDCD). This further reveals that the crosstalk between different programmed cell death is a complex and dynamic process (Supplementary Fig. S13B). GSEA analysis showed that the high expression of DLAT was closely related to TCA cycle, glucose metabolism, lipid metabolism and pentose phosphate pathway (Supplementary Fig. 13C). Additionally, DLAT may be involved in the proliferation and metastasis of HCC (Supplementary Fig. S13C).

Fig. 12
figure 12

DLAT was highly expressed in HCC tissues. (A-B) qRT-PCR and WB confirmed high expression level of DLAT in HCC tissues. (C) Representative images of IHC analyses showing upregulation of DLAT in HCC compared with normal liver tissues in HPA database (C) and real-world cohort (D). (top: × 100 magnification, scale bar, 200 μm; bottom: × 400 magnification, scale bar, 50 μm).

Full size image

PCDI-associated signatures in Pan-cancer

Based on the expression data in 33 cancer types presented in the TCGA database, our analysis revealed a prevailing elevation in the expression levels of these five gene signatures across the majority of cancers, surpassing their corresponding normal tissue counterparts (Supplementary Fig. S14). Forest plots based on univariate Cox regression showed that, the dysregulation of DLAT exhibited a significant correlation with the prognostic outcomes of four distinct tumor types, including adrenocortical carcinoma (ACC, HR = 2.266, P = 0.039), breast invasive carcinoma (BRCA, HR = 1.457, P = 0.023), colon adenocarcinoma (COAD, HR = 0.622, P = 0.019), kidney renal clear cell carcinoma (KIRC, HR = 0.519, P < 0.001), brain lower grade glioma (LGG, HR = 1.696, P = 0.005), liver hepatocellular carcinoma (LIHC, HR = 1.721, P = 0.002), pancreatic adenocarcinoma(PAAD, HR = 1.760, P = 0.008), and rectum adenocarcinoma (READ, HR = 0.326, P = 0.012). In addition, other signatures have also shown prognostic value in a variety of tumors (Supplementary Fig. S15). These results indicated that these PCDI signatures might be involved in tumorigenesis and tumor development.

Discussion

The limitations of the current tumor-node-metastasis (TNM) classification and Barcelona Clinic Liver Cancer (BCLC) staging classification system of the International Union Against Cancer (7th Edition) impede the application to offer individual management to HCC patients. While the TNM/ BCLC stratified strategies described the extent of disease, staging in HCC is of limited prognostic value, as the clinical management of HCC is primarily determined by disease histology, while not molecular biological characteristics31,32. A robust risk stratification strategy is promising for personalized medicine or precision medicine, and further improves chances for favorable outcome benefited from the intrinsic and interpatient phenotypic or prognostic heterogeneity. Since several novel regulatory mechanisms have been reported as involved in the PCD biological process, PCD in mammals has been implicated in several disease states including cancer, autoimmune disease, neurodegenerative disease and aging. Considering the cross-talk between different PCD patterns or pathways can be activated, inhibited or altered, resulting in distinct clinical outcome and response to treatments. However, it is important to recognize that the role and the regulation of PCD in HCC or pan-cancer is very redundant and complex, which partly contributed to the disappointing outcome of clinical application.

Recent studies have highlighted several PCD patterns, such as cuproptosis, ferroptosis, pyroptosis and necroptosis, which are closely involved in progression of liver disease including HCC, alcoholic fatty liver disease, and non-alcoholic fatty liver disease, and represents a potential target for HCC diagnosis and treatment33,34,35. Numerous predictive gene expression-based models have been adopted based on various ML approaches, in response to the increasing advancements in high-throughput sequencing techniques and computational biology algorithms36,37,38. Benefiting from the integration of multiple feature selection algorithms, we developed and validated a PCDI model using the optimal subset (HTRA2, SLC2A1, FTL, G6PD, and DLAT), which was stably associated with prognosis in four independent HCC cohorts, and offered an available resource for developing disease biomarkers and therapy targets. Our study systematically revealed links between PCD index and HCC prognosis, and benefits from PCDI-based strategy classification. The index identified different subgroups of patients with low and high risk of mortality, providing an individualized risk assessment beyond the existing clinical and pathological parameters. The further findings suggested that the integration of tumor genomics and clinicopathologic features improved the C-statistic and metrics of risk reclassification. Whether such a molecule-based classifier, when combined with conventional clinicopathological characteristics and other previously reported index, will help to provide refinement for patient stratification strategies remains to be fully demonstrated in prospective clinical trials and independent cohorts. Moreover, through both large-scale survival analysis and a real-world cohort, we indicated that high DLAT expression in HCC tissues, which predicted an unfavorable prognosis and may therefore serve as a progression indicator in HCC.

Moreover, our study has provided evidence indicating a strong correlation between PCDI and a repressive TME as well as dysregulated metabolic processes. Specifically, higher PCDI indicated that increased infiltration of M2 macrophages and regulatory T cells (Tregs), which display diminished functionality in all metabolic pathways except those related to nucleotides. Notably, previous research has consistently reported a significant association between various PCD patterns and both TME and metabolic dysregulation. Metabolism and programmed cell death (PCD) are interrelated structures and processes that play critical roles in maintaining normal cellular functions and tissue integrity. Mitochondrial metabolism plays a crucial role in generating ATP, the primary energy source for the cell39. Metabolic pathways such as glycolysis, fatty acid oxidation, and the tricarboxylic acid (TCA) cycle also produce intermediate molecules that can be used for biosynthesis and cell signaling40. Recent studies have identified several links between PCD and metabolism. For instance, metabolic dysfunction or dysregulation can trigger PCD by different mechanisms, including the generation of reactive oxygen species (ROS), mitochondrial dysfunction, or the activation of death receptors41,42. Conversely, PCD can also induce metabolic changes by affecting key signaling pathways, such as mTOR, AMPK, or HIF-1α43,44. Several types of research have demonstrated that reprogramming metabolic pathways can control PCD. For example, inhibiting glycolysis and increasing mitochondrial metabolism attenuates cell death by apoptosis or necrosis45. Also, autophagy-triggering agents, such as rapamycin, activate cellular metabolism to enhance ATP generation and promote cell survival46. Indeed, all of the PCD patterns integrate metabolic processes to enable precise regulation of cell death during development, tissue remodeling, and defense against infections and diseases. A shared characteristic among these cell death pathways is their initiation by disruptions in cellular metabolism, which arise from either the depletion or excess of various metals or nutrients. This highlights the complex interaction between cellular metabolism and the regulation of diverse cell death modalities. Future research may elucidate further mechanistic similarities and distinctions among these pathways of cellular demise.

Apoptosis, the best-characterized form of PCD, involves a series of metabolic modifications in the mitochondria, and endoplasmic reticulum that lead to the activation of proteases, cleavage of structural proteins, and fragmentation of DNA40. During apoptosis, mitochondrial metabolism is altered, leading to decreased ATP production, enhanced vicious cycle of mitochondrial ROS (excessive reactive oxygen species) augmentation and autophagy flux, and a permeability increase of the inner mitochondrial membrane. Changes in metabolic fluxes and ROS levels modulate the expression or activity of pro-apoptotic and anti-apoptotic proteins, which transduce signals downstream of apoptotic receptors and activate caspases47,48.

Necroptosis, a type of programmed necrosis, entails the activation of receptor-interacting protein kinase 1 (RIPK1) and RIPK3 and the formation of necrosomes that initiate the assembly of membrane pores and the subsequent release of damage-associated molecular patterns (DAMPs)49. Metabolites such as glucose, amino acids, and fatty acids can regulate necroptosis by affecting the accumulation of reactive oxygen species (ROS), which can lead to oxidative damage and mitochondrial dysfunction50,51. Autophagy is another form of PCD that involves the lysosomal degradation of cytosolic components, proteins, and organelles. Autophagy is dependent on the cellular metabolic status, and nutrient-sensing pathways such as mTOR, AMPK, and SIRT1 modulate the initiation and progression of autophagic flux52. During nutrient deprivation or hypoxia, cells induce autophagy to recycle amino acids, fatty acids, and nucleotides as substrates for ATP generation and cell survival53. Finally, pyroptosis is an inflammatory form of PCD that results from the activation of caspase-1 and the formation of inflammasomes. Pyroptosis relies on the activation of glycolysis and mitochondrial respiration to fuel the production of interleukins and cytokines that amplify the immune response to infections and inflammation54. In conclusion, metabolism is closely related to PCD, and metabolic reprogramming can affect cell death decision and fate. A better understanding of the metabolic regulation of PCD can have significant impacts on pathological conditions, including cancer, neurodegenerative diseases, and autoimmune disorders, where dysregulated PCD contributes to disease progression.

Crosstalk between different programmed cell death (PCD) patterns or pathways is a complex mechanism that plays a central role in the pathogenesis and progression of many tumor types, including HCC47,55. The crosstalk between PCD patterns in HCC can be either synergistic or antagonistic and can impact different aspects of tumor growth and treatment response10,56. Autophagy can promote apoptosis by removing damaged proteins and organelles, thereby reducing the anti-apoptotic effects of these molecules. Moreover, therapeutic agents that inhibit autophagy can enhance the effect of apoptosis-inducing chemotherapy in HCC. Under stress sufficiency, autophagy activation can be processed as a protective mechanism under, such as hypoxia, and can prevent the progression of necrosis11,57. In contrast, necrosis can trigger the release of intracellular contents, including danger-associated molecular patterns (DAMPs), that stimulate autophagy and promote tumor growth58. Significantly, in certain circumstances, such as the loss, inhibition, or mutation of elements within the canonical pathway, a transition occurs where one cell death pathway is substituted by another. Overall, the crosstalk between different PCD pathways in HCC is a complex interplay between tumor cells and their microenvironment. A better understanding of the molecular mechanisms underlying PCD pattern crosstalk in HCC could lead to a better exploration of the pathogenesis of HCC and thereby the development of novel therapeutic approaches to the disease.

Further, we revealed that a higher frequency of TP53 mutations was associated with a higher PCDI, which was considered to be an important regulator of PCD in HCC. GSEA analysis showed that the mutant status of TP53 was associated with high glycolytic activity and disturbed lipid metabolism. TP53/p53 plays a transcription-dependent and -independent role in the regulation of apoptosis, autophagy, ferroptosis, metabolism, cellular oxidative-redox status, and other PCD patterns59,60,61,62. Furthermore, the deletion or mutation of TP53 has a significant impact on the recruitment and functionality of T cells, thereby creating an inhibitory immune microenvironment that ultimately facilitates the evasion of tumor cells from immune surveillance63,64. In conclusion, the interplay interactions between programmed cell death, mutation, and metabolism was complex and bidirectional. For example, mutations in metabolic genes can lead to metabolic reprogramming that promotes cell survival and suppresses PCD pathways. Conversely, the dysregulation of PCD pathways can lead to the accumulation of abnormal cells with mutations that promote tumor development and progression.

Understanding the complex interactions and regulatory mechanisms between PCD, mutation, and metabolism is critical to developing effective therapies for HCC and other diseases. By targeting the key regulatory pathways involved in PCD and metabolism, new and more effective treatments can be developed that can selectively promote cancer cell death while preserving normal cells.

Conclusion

In conclusion, our study recognized the comprehensive landscape of PCD-associated molecular subtypes. A PCD-based risk stratification strategy was presented, which was applied to predict the long-term outcome, profile the clinicopathological features distribution, characterize the TME and reveal the heterogeneity and signatures of multi-omics panoramic landscapes, to extend previous classification strategies for subtyping HCC patients. Understanding the complex regulation network rather than a linear cascade of PCD in HCC could lead to the improvements in the treatment of patients relying on the identification of novel predictive markers and treatment targets for support of individualized clinical decisions and targeted treatment strategies.