An Immunogenic Cell Death-Related Classification Predicts Prognosis and Response to Immunotherapy in Head and Neck Squamous Cell Carcinoma

Abstract

Immunogenic cell death (ICD) has been classified as a form of regulated cell death (RCD) that is sufficient to activate an adaptive immune response. Accumulating evidence has demonstrated the ability of ICD to reshape the tumor immune microenvironment through the emission of danger signals or DAMPs, which may contribute to the immunotherapy. Currently, identification of ICD-associated biomarkers that stratify patients according to their benefit from ICD immunotherapy would be of great advantage. Here, we identified two ICD-associated subtypes by consensus clustering. ICD-high subtype was associated with the favorable clinical outcomes, abundant immune cell infiltration, and high activity of immune response signaling. Besides, we established and validated an ICD-related prognostic model that predicted the survival of HNSCC and was associated with tumor immune microenvironment. In conclusion, we established a new classification system of HNSCC based on ICD signatures. This stratification had significant clinical outcomes for estimating prognosis, as well as the immunotherapy of HNSCC patients.

Keywords: head and neck squamous cell carcinoma; immunogenic cell death; immunotherapy; prognosis; tumor microenvionment.

Figure 1

Identification of ICD-associated subtypes by consensus clustering. (A) Protein–protein interactions among the ICD-associated genes; (B) Heatmap shows 34 ICD gene expression profiles among normal and HNSCC samples in TCGA database; (C) Immunohistochemistry validated the expression pattern of CALR in normal and tumor sample; (D) Heatmap depicts consensus clustering solution (k = 2) for 36 genes in 502 HNSCC samples; (E) Delta area curve of consensus clustering indicates the relative change in area under the cumulative distribution function (CDF) curve for k = 2 to 10; (F) Heatmap of 34 ICD-related gene expressions in different subtypes. Red represents high expression and blue represents low expression; (G) Kaplan–Meier curves of OS in ICD-high and ICD-low subtypes. *P < 0.05, **P < 0.01, ***P < 0.001, & ****P < 0.0001.

Figure 2

Identification of differentially expressed genes (DEGs) and underlying signal pathways in different subtypes. (A) Volcano plot presents the distribution of DEGs quantified between ICD-high and ICD-low subtypes with threshold of |log2 Fold change| > 1 and P < 0.05 in TCGA cohort; (B) Heatmap shows the DEG expression in different subtypes; (C) Dots plot presents the KEGG and GO signaling pathway enrichment analysis. The size of the dot represents gene count, and the color of the dot represents – log₁₀ (p. adjust-value); (D) GSEA analysis determines the underlying signal pathway between ICD-high and ICD-low subtypes.

Figure 3

Comparison of somatic mutations between different ICD subtypes. (A, B) Oncoprint visualization of the top ten most frequently mutated genes in ICD-high subtype (A), and ICD-low subtype (B).

Figure 4

Immune landscape of ICD-high and ICD-low subtypes. (A) Violin plots show the median, and quartile estimations for each immune score, and tumor purity score; (B) Relative proportion of immune infiltration in ICD-high and ICD-low subtypes; (C) Violin plot visualizes significantly different immune cells between different subtypes; (D, E) Box plots present differential expression of multiple immune checkpoints (D), and HLA genes (E) between ICD-high and ICD-low subtypes. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.

Figure 5

Construction and validation of the ICD risk signature. (A) Univariate Cox analysis evaluates the prognostic value of the ICD genes in terms of OS; (B) Lasso Cox analysis identified 12 genes most associated with OS in TCGA dataset; (C) Risk scores distribution, survival status of each patient, and heatmaps of prognostic 12-gene signature in TCGA database; (D, E) Kaplan–Meier analyses demonstrate the prognostic significance of the risk model in TCGA and GSE41613 cohort.

Figure 6

The association of ICD risk signature with tumor microenvironment. (A, B) Scatter plots show the correlation of risk score with the infiltration of CD8, activated NK cell, and activated CD4 memory cell (A), which was further validated by GSE41613 cohort (B); (C) Box plot presents the association of ICD risk score with immunotherapy response; (D, E) Univariate and multivariate Cox analyses evaluate the independent prognostic value of ICD risk signature in HNSCC patients.