Establishing a signature based on immunogenic cell death-related gene pairs to predict immunotherapy and survival outcomes of patients with hepatocellular carcinoma

Sample data collection and identification of immunogenic cell death (ICD)-related genes

RNA sequencing results, clinical features of HCC patients, and normal liver specimens were available through The Cancer Genome Atlas (TCGA) project and Genotype-Tissue Expression (GTEx) database. After excluding patients with no clinical data, repeated data, incomplete clinical data, and a follow-up time of 0 days, 374 HCC patients and 160 normal samples were finally obtained. For the validation set, the GSE14520 dataset was downloaded from Gene Expression Omnibus (GEO) database. Additionally, 34 ICD-related genes were obtained from previous articles [12], and are listed in Supplementary Table 1. The distribution of clinical features of HCC cohorts were detailed in Table 1.

Identification of ICD-related gene pairs

Pairs of ICD-related genes were screened by cyclically singly paired. If the expression level of ICD-related gene A is greater than that of ICD-related gene B, it is recorded as 1; otherwise, it is recorded as 0, and the 0-or-1 matrixes are established. We consider that there was no relationship between pairs and prognoses if ICD-related gene pair = 0 or 1 because the patient’s survival outcome cannot be correctly predicted due to pairs without a certain rank. When the number of pairs equal to 1 or 0 was >20% and <80% of the total number of pairs, the match was deemed valid.

Development and verification of ICD-related prognostic signature

Univariate cox regression was employed to retrieve OS-associated gene pairs. These gene pairs were then used to develop an ICD-related gene pair signature capable of predicting HCC patient prognosis through a LASSO regression method. The formula developed based on this analysis was then established as follows:

Risk score = regression coefficient (genei) × expression value (genei).

By setting the median value of the risk score calculated in the training set as the threshold, patients in the training and validation sets (GSE14520) were stratified into two risk subgroups. Survival analysis was conducted to analyze differences in survival between two risk HCC patients. The subsequent receiver operator characteristic (ROC) analysis and C-index were employed to assess the prognostic accuracy of the developed signature.

Construction of a nomogram

To further confirm the applicability of the signature, we investigated the association between the risk score and clinical data, including age, sex, grade, surgical procedure, recurrence status, vascular tumor cell type, pathologic stage, and tumor status. The prognostic utility of the ICD-related signature was also examined through univariate analyses in both sets, and its ability to independently predict patient prognosis was examined via a multivariate approach. A nomogram including the risk score and various clinical traits (surgery procedure and stage) was constructed using the “regplot” package. The predictive capability of the nomogram was evaluated by AUC values, calibration plots, and decision-curve analysis (DCA). The area under the curve (AUC) of the ROC was used to assess the accuracy of the nomogram in 3- and 5-year survival predictions and the predictive prognostic performance between risk score and TNM stage were compared.

Assessment of immune landscape

The interplay between the risk scores and the immune microenvironment in HCC patient tumors was assessed via the ESTIMATE algorithm, which examined the predicted stromal and immune cell content in these tumor samples. Calculations of immune infiltration statuses among the HCC patients were conducted using by CIBERSORT algorithm. A Spearman correlation analysis was conducted to investigate the correlations between risk score and the infiltration of immune cells by the CIBERSORT algorithm.

Immunotherapeutic response analyses and drug efficacy assessment

To investigate the predictive ability of risk scores in the benefit of immunotherapeutic treatment, we first the differences in the expression of immune checkpoints between the two risk subgroups using “ggpubr” R package. Secondly, the TIDE method was utilized to predict cancer immunotherapy response and tumor immune escape in low- and high-risk populations [13]. Moreover, we analyzed the immunophenoscore (IPS) of HCC samples in the TCIA database. The correlations between tumor mutation burden (TMB) and risk score were analyzed using “ggExtra”, “ggplot2” and “ggpubr” packages in R software. The optimal cutoff value of TMB in survival data was identified through the function surv_cutpoint of the “survminer” R package, and then patients were stratified into low- and high-TMB subgroups. The OS of CRC samples between two subgroups was compared by the Kaplan-Meier method. Next, a combined survival analysis was performed for TMB and risk score. To evaluate the response of the two risk groups to different drugs, the R package “pRRophetic” was employed to evaluate the half inhibitory concentration (IC50) of some common chemotherapeutic drugs in every HCC specimen.

Functional enrichment analysis

To understand the biological functions of risk-related differential genes and the potential signaling enrichment pathways, the functional enrichment analyses were conducted with the package “clusterProfiler”, and presented as bar graphs and bubble plots, respectively.

Statistical analysis

R (v4.0.3; http://www.Rproject.org) was used for statistical comparisons. P-value < 0.05 would be considered statistically significant.

Availability of data and materials

Publicly available datasets were analyzed in this study. This data can be found at TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/).

Identification of ICD-related gene pairs and establishment of the risk model in HCC

This study progressed according to the flow chart (Figure 1). A total of 89 ICD-related gene pairs were screened with the algorithm described in “Methods”. Of the 89 ICD-related gene pairs identified above, 15 were evidently associated with the OS of HCC patients in the training set (Figure 2A). These survival-associated ICD-related gene pairs were then used to generate a prognostic signature via a LASSO Cox regression approach (Figure 2B, 2C), ultimately identifying 8 gene pairs for inclusion in the developed signature. The coefficients of 8 ICD-related gene pairs were employed to calculate the risk score (Table 2). Its prognostic value was next assessed in the training set, with HCC patients being stratified into low- and high-risk groups. HCC patients with the high-risk score had a worse OS (P < 0.05; Figure 2D). The same formula and cut-off threshold were then applied to the validation set, which similarly revealed worse OS (Figure 2E). The AUCs of the training and validation sets were 0.720 and 0.710 at 3 years, 0.717 and 0.707 at 5 years, respectively, indicating that the signature has good predictive efficacy (Figure 2F, 2G). We then calculated the C-index of risk score and several clinical parameters. The highest C-index of the risk score affirmed the predictive utility of our signature in both training and validation sets (Figure 2H, 2I). Additionally, we compared the risk model with existing signatures in HCC. ROC curves indicated that our signature achieved significantly favorable predictive power compared with previously published prognostic models (Supplementary Figure 1). Risk score and OS distributions were shown in Figure 2J and it was noticed that the risk of death gradually increased with the increment of the risk score. The same analyses were then performed in the validation set and similar distribution was found (Figure 2K).

Figure 1. The flow chart of this study.

Figure 2. Construction and verification of ICD-related prognostic gene signature. (A) Forest plots showing the results of the Cox univariate regression of 15 ICD-related prognostic gene pairs. (B) The LASSO coefficient profile of 8 ICD-related prognostic gene pairs. (C) Selection of optimal LASSO model parameters for HCC patients (λ). (D, E) Kaplan-Meier survival curves of patients in low- and high-risk groups in training (D) and validation (E) sets. (F, G) The AUC values of the signature in the training (F) and validation sets (G). (H, I) C-index of the risk score and other clinical traits in the training (H) and validation sets (I). (J, K) Risk score and OS distributions were assessed in the training (J) and validation sets (K).

Development and validation of a nomogram

We explored the relationship between the risk score and clinicopathological traits, including age, gender, surgical procedure, grade, vascular tumor cell type, recurrence status, tumor status, and stage. Significant correlations were observed between risk score and surgical procedure, vascular tumor cell type, recurrence status, tumor status, and stages (Figure 3A–3E). HCC patients who underwent segmentectomy had evidently higher risk scores than those who underwent lobectomy (Figure 3A). Significantly higher risk scores were observed in the Macro and Micro groups (Figure 3B). High-risk scores were more common in patients with recurrence status (Figure 3C) and with tumor status (Figure 3D). In addition, patients in stage pathological III-IV had higher risk scores than those in stage I-II (Figure 3E). Univariate and multivariable Cox regression analyses demonstrated the independent predictive roles of surgery, stage, and risk score in the training cohort (Figure 3F, 3G). The independent predictive roles of stage and risk score were further confirmed in the validation set (Figure 3H, 3I). To enhance the clinical use of the signature, risk score, surgery procedure, and stage were combined to build a prognostic nomogram (Figure 3J). The combined nomogram showed an excellent predictive ability and performed better than any TNM stage in predicting 3- and 5-year survival (Figure 3K, 3L). Calibration curves demonstrated that the nomogram performed well (Figure 3M). Meanwhile, the DCA curve demonstrated the contribution of risk and nomogram in clinical decision-making (Figure 3N).

Figure 3. Generation and validation of a nomogram scoring system. (A–E) The relationship between the risk score, surgical procedure, vascular tumor cell type, recurrence status, tumor status, and stages. (F, G) Forest plots showing the results of the univariate (F) and multivariate (G) Cox analysis in the training set. (H, I) Forest plots showing the results of the univariate (H) and multivariate (I) Cox analysis in the validation set. (J) Nomogram predicting the 3-year and 5-year overall survival of HCC patients. (K, L) Comparison of the 3- (K) and 5-year AUC (L) of nomogram and TNM staging system. (M) The calibration plots of the nomogram at 3 and 5 years. (N) The DCA curves of the nomogram at 3 and 5 years.

Assessment of immune landscape

The ESTIMATE algorithm was further used to process HCC patient sample data to compare the relative contribution of immune and stromal cells to the obtained patient samples. Compared with low-risk patients, patients with high-risk scores had significantly lower immune and stromal scores and higher tumor purity (P < 0.001; Figure 4A–4D). Based on the CIBERSORT algorithm, we observed increased CD8+ T cell, plasma cells, activated memory CD4+ T cell, M1 macrophage, and T follicular helper cell infiltration in low-risk patients, whereas high-risk patients exhibited enhanced resting memory CD4+ T cell, M0 macrophage, M2 macrophage, and neutrophils infiltration (P < 0.05; Figure 4E). The scatter plot (Figure 4F–4L) further showed that high-risk score was positively associated with several infiltrations of immune cells, such as M0 macrophage, M2 macrophage, and neutrophils (Figure 4F–4H), while the infiltration level of CD8+ T cell, activated memory CD4+ T cell, M1 macrophage, and Tfh cells was associated with the high-risk group (Figure 4I–4L).

Figure 4. Assessment of immune cell infiltration and the immune microenvironment in different subgroups. (A–D) Differences in immune scores, stromal scores, and tumor purity between high-risk and low-risk groups. (E) Differences in immune cells infiltration between high-risk and low-risk groups. (F–L) Correlations between the risk score and immune cells infiltration. ^*, ^**, and ^*** represent p < 0.05, 0.01, and 0.001, respectively.

Immunotherapeutic response analyses and drug efficacy assessment

We compared the expression of major immune checkpoints including HAVCR2, PD-1, and PD-L1, between the two subgroups. Overall, most of the immune checkpoints were substantially elevated in the low-risk group (Figure 5A), indicating that patients with low-risk scores may achieve better ICIs therapy results. Moreover, we use the TIDE algorithm to evaluate the likelihood of HCC benefiting from ICI therapy. The results demonstrated that the TIDE score was higher in high-risk patients (Figure 5B), suggesting tumors in high-risk patients could acquire immune escape more easily. Furthermore, we compared the immunotherapy efficacy of the two risk subgroups by IPS scores and found that in the anti-CTLA4 single-drug group, the anti-PD1 single-drug group, and the anti-CTLA4 and anti-PD1 combination group, the IPS scores were higher in the low-risk group (Figure 5C–5F), suggesting that patients with low-risk HCC may have greater sensitivity to immunotherapy. Growing evidence suggests that TMB may determine the individual response to cancer immunotherapy. Correlation analysis demonstrated that the risk scores were negatively associated with TMB (Figure 5G), indicating low-risk patients may respond better to immunotherapy. In addition, K-M survival analysis demonstrated that low-TMB patients enjoyed a much longer survival time than their counterparts (Figure 5H). Combining TMB and risk score allowed us to classify patients into four groups. The combined survival analysis showed that the high TMB and high- risk groups had the worst prognoses, and, conversely, the low TMB and low-risk groups had the best prognoses (Figure 5I).

Figure 5. Comparison of the immunotherapeutic response of low- and high-risk HCC patients. (A) Comparison of immune checkpoint expression between the two subgroups. (B) Differences in TIDE scores in low- and high-risk individuals. (C–F) Differences in IPS scores in low- and high-risk individuals. (G) Differences in TMB in low- and high-risk individuals. (H) Kaplan-Meier survival curves of HCC patients in different TMB groups. (I) Kaplan-Meier survival curves in different risk scores and TMB subgroups.

To identify a drug therapy target, a crucial way is to clarify the correlation between ICD-related signatures and antitumor drugs. As shown in Figure 6, the IC50 values of Temsirolimus, Bortezomib, Metformin, Paclitaxel, and Sunitinib were higher in the high-risk subgroup (Figure 6A–6E), while the IC50 value of AKT inhibitor VIII was higher in the low-risk subgroup (Figure 6F).

Figure 6. Chemotherapy sensitivity and functional enrichment analyses. (A–F) Boxplots of the IC50 values of the Temsirolimus, Bortezomib, Metformin, Paclitaxel, Sunitinib, and AKT inhibitor VIII between different risk subgroups. (G, H) Heatmap (G) and volcano plots (H) demonstrate DEGs between two risk subgroups. (I) GO terms of DEGs among different risk subgroups. (J) KEGG enrichment analyses of DEGs among different risk subgroups.

GO and KEGG enrichment analyses

To characterize the functional enrichment in risk groups, we compared the mRNA expression difference between two risk subgroups to identify DEGs (|log2 (fold change, FC) | > 1 and adjusted P < 0.05), and a total of 314 DEGs were identified (Figure 6G, 6H). As shown in Figure 6I, GO terms indicated that these DEGs were mainly correlated with antigen binding, MHC protein binding, and immunoglobulin receptor binding. Additionally, KEGG pathway analysis demonstrated that these DEGs were mainly associated with immune-related pathways, including the T cell receptor signaling pathway, Cytokine-cytokine receptor interaction, Cell adhesion molecules, and Th17 cell differentiation (Figure 6J).