Integrated multiomics analysis identified comprehensive crosstalk between diverse programmed cell death patterns and novel molecular subtypes in Hepatocellular Carcinoma

Abstract

Hepatocellular carcinoma (HCC) is a highly aggressive malignancy with increasing global prevalence and is one of the leading causes of cancer-related mortality in the human population. Developing robust clinical prediction models and prognostic stratification strategies is crucial for developing individualized treatment plans. A range of novel forms of programmed cell death (PCD) plays a role in the pathological progression and advancement of HCC, and in-depth study of PCD is expected to further improve the prognosis of HCC patients. Sixteen patterns (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) related to PCD were collected from the literatures and used for subsequent analysis. Supervised (Elastic net, Random Forest, XgBoost, and Boruta) and unsupervised (Nonnegative Matrix Factorization, NMF) clustering algorithms were applied to develop and validate a novel classifier for the individualized management of HCC patients at the transcriptomic, proteomic and single-cell levels. Multiple machine learning algorithms developed a programmed cell death index (PCDI) comprising five robust signatures (FTL, G6PD, SLC2A1, HTRA2, and DLAT) in four independent HCC cohorts, and a higher PCDI was predictive of higher pathological grades and worse prognoses. Furthermore, a higher PCDI was found to be correlated with the presence of a repressive tumor immune microenvironment (TME), as determined through an integrated examination of bulk and single-cell transcriptome data. In addition, patients with TP53 mutation had higher PCDI in comparison with TP53 WT patients. Three HCC subtypes were identified through unsupervised clustering (NMF), exhibiting distinct prognoses and significant biological processes, among the three subtypes, PCDcluster 3 was of particular interest as it contained a large proportion of patients with high risk and low metabolic activity. Construction and evaluation of the Nomogram model was drawn based on the multivariate logistic regression analysis, and highlighted the robustness of the Nomogram model in other independent HCC cohorts. Finally, to explore the prognostic value, we also validated the frequent upregulation of DLAT in a real-world cohort of human HCC specimens by qPCR, western blot, and immunohistochemical staining (IHC). Together, our work herein comprehensively emphasized PCD-related patterns and key regulators, such as DLAT, contributed to the evolution and prognosis of tumor foci in HCC patients, and strengthened our understanding of PCD characteristics and promoted more effective risk stratification strategies.

Keywords: Hepatocellular carcinoma; Integrated multi-omics analysis; Machine learning model; Molecular subtype; Programmed cell death.

Fig. 1

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

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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.

Fig. 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).