A novel prognostic model based on immunogenic cell death-related genes for improved risk stratification in hepatocellular carcinoma patients

Abstract

Purpose: Hepatocellular carcinoma (HCC) is a prevalent primary malignant tumor with increasing incidence and mortality rates in recent years. The treatment options for advanced HCC are very limited. Immunogenic cell death (ICD) plays an important role in cancer, in particular immunotherapy. However, the specific ICD genes and their prognostic values in HCC remain to be investigated.

Methods: The TCGA-LIHC datasets were obtained from TCGA database, LIRI-JP datasets were obtained from ICGC database, and immunogenic cell death (ICD) genes datasets were obtained from previous literature. WGCNA analysis identifies ICD-related genes. Functional analysis was used to investigate the biological characteristics of ICD-related genes. Univariate Cox analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis was used to select prognostic ICD-related genes and construct a prognostic risk score. Prognostic independence of ICD risk scores was determined by univariate and multivariate Cox regression analyses. A nomogram was then constructed and the diagnostic value was assessed using decision curve analysis. Immune infiltration analysis and drug sensitivity analysis were used to investigate immune cell enrichment and drug response in HCC patients classified as low or high risk based on their risk score.

Results: Most of the ICD genes were differentially expressed in normal and HCC patients, and some ICD genes were differentially expressed in different clinical groups. A total of 185 ICD-related genes were identified by WGCNA. Prognostic ICD-related genes were selected using a univariate Cox analysis. A model comprising nine prognosis ICD-related gene biomarkers was developed. Patients was divided into high-risk and low-risk groups, and patients in high-risk groups had poorer outcomes. Meanwhile, the reliability of the model was verified by external independent data. The Independent prognostic value of the risk score in HCC was investigated by univariate and multivariate Cox analyses. Diagnostic nomogram was constructed to predict prognosis. Through immune infiltration analysis, we found that some innate and adaptive immune cells were significantly different between low- and high-risk groups.

Conclusion: We developed and validated a novel prognostic predictive classification system for HCC based on nine ICD-related genes. In addition, immune-related predictions and model could help predict the outcomes of HCC and could provide a reference for clinical practice.

Keywords: Hepatocellular carcinoma; Immunogenic cell death; LASSO; Prognosis; WGCNA.

Fig. 1

Altered ICD genes in HCC patients. a The expression of ICD genes in HCC tumor and normal tissues. b–d The expression of ICD genes in different HCC clinical characters. e Waterfall plots show somatic mutations of ICD genes in HCC patients. f Circos plot shows the distribution of ICD genes in chromosome. *, **, *** and **** indicate P value less than 0.05, 0.01, 0.001 and 0.0001, respectively

Fig. 2

WGCNA on ICD ssGSEA score. a Clustering dendrogram of samples with trait heatmap. b Analysis of network topology for various soft-thresholding powers. The left panel shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis). c Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned module colors. d Module–trait associations. Each row corresponds to a module, and each column corresponds to a trait. Each cell contains the corresponding correlation and p value. The table is color coded by correlation according to the color legend. e A scatterplot of gene significance (GS) for ICD ssGSEA score vs. module membership (MM) in the saddlebrown (left) and darkorange (right) module. There is a highly significant correlation between GS and MM in these two modules

Fig. 3

Construction and validation of the ICD risk signature model. a Lasso Cox analysis identified nine genes most associated with OS in the TCGA dataset. b Risk score distribution, survival status of each patient and heatmaps of prognostic nine-gene signature in the TCGA database. c Kaplan–Meier analyses demonstrate the prognostic significance of the risk model in the TCGA cohort. d Time-dependent receiver operating characteristic (ROC) curve for 1-, 3-, and 5-year survival of HCC patients. e qPCR showed the difference expression of PNLIPRP2, FKBP6, FGF9 and RAMP3 in HCC tumor tissues (T) and normal tissues adjacent to the tumor (N)

Fig. 4

Clinical characteristics and prognostic analysis of ICD risk signature in the training cohort. a Box plots present differential risk score of multiple clinical characters in the training cohort. b, c Univariate and multivariate Cox analyses evaluate the independent prognostic value of ICD risk signature in the training cohort. *, **, *** and **** indicate p value less than 0.05, 0.01, 0.001 and 0.0001, respectively

Fig. 5

An ICD-related gene signature-based prediction model. a The nomogram prediction model for HCC patients’ OS. b Calibration graphs of the training set’s 1-year, 2-year, 3-year and 5-year survival probabilities. c 1-year, 3-year and 5-year survival benefit in the training cohort

Fig. 6

The association of ICD risk signature with tumor microenvironment. a Relative proportion of EstimateScore and ImmuneScore in ICD-high and ICD-low groups. b Scatter plots show the correlation of risk score with the infiltration of EstimateScore and ImmuneScore. c The differential estimated proportion of 36 XCELL immune cell types in ICD-related subtypes. The central line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). *, **, *** and **** indicate p value less than 0.05, 0.01, 0.001 and 0.0001, respectively

Fig. 7

Drug sensitivity analysis of HCC patients. Lollipop chart shows the correlation of risk score with drug IC50 from GDSC (a) and CCLE (b, c) database. Box plots present differential estimated AUC value of multiple drug from GDSC (d) and CCLE (e, f) database. *, **, *** and **** indicated p value was less than 0.05, 0.01, 0.001 and 0.0001, respectively