Pyroptosis patterns and immune infiltrates characterization in head and neck squamous cell carcinoma.

BACKGROUND: Pyroptosis plays an essential role in tumor immune responses and inflammation related to chemotherapy. Herein, we studied the characteristic patterns of pyroptosis in head and neck squamous cell carcinoma (HNSCC) to determine their prognostic and therapeutic effects.
METHODS: Consensus clustering analysis was performed to classify patients into pyroptosis or gene clusters. A novel pyroptosis score was constructed by principal component analysis. Kaplan-Meier survival curves were used to show the prognostic value. We also assessed the functional enrichment, tumor mutation burden, immune cell infiltration, and the sensitivity to chemotherapy and immunotherapy between high and low pyroptosis score group.
RESULTS: Two distinct pyroptosis clusters were defined based on the mRNA expression profiles of PRGs, which were related to immune activation in HNSCC. Notably, a pyroptosis score was constructed according to different expression gene signatures, and then, each HNSCC patient was classified into a low or high pyroptosis score group. Patients with low pyroptosis scores had better immunotherapeutic responses and higher sensitivities to chemotherapeutic agents (paclitaxel, docetaxel, and gemcitabine). Kaplan-Meier survival curves showed that the pyroptosis patterns were independent prognostic indicators regardless of the level of tumor mutation burden.
CONCLUSIONS: Pyroptosis plays an essential role in immune infiltration in HNSCC. Quantifying the pyroptosis score of individual patients might suggest prognostic, immunotherapeutic, and chemotherapeutic strategies for HNSCC.

INTRODUCTION

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide; approximately 600 000 new cases are diagnosed each year. The primary treatments of HNSCC depend on the clinical stage and generally include surgery, radiotherapy, and chemotherapy. There are limited treatments for advanced HNSCC patients, and the outcome is poor; the 5‐year overall survival rate is only 63%. The introduction of immunotherapy and targeted therapy substantially improved HNSCC outcomes. Immune checkpoint inhibitor (ICI) antibodies targeting programmed death ligand‐1 (PD‐L1) or programmed cell death protein 1 (PD1) were approved by the US Food and Drug Administration for use in combination with chemotherapy as first‐line treatment. Nevertheless, only 10%–20% of HNSCC patients benefit from PD‐L1/PD1 antibody therapy, , , suggesting that treatment strategies for HNSCC need to be improved.

Pyroptosis is a type of programmed cell death that is highly proinflammatory. Pyroptosis is mediated by the gasdermin (GSDM) family, including GSDMA, GSDMB, GSDMC, GSDMD, GSDME and pejvakin (PJVK). As shown in Figure 1A, GSDMs are usually cleaved to liberate an N‐terminal fragment by the activated inflammatory caspase family. There are two common mechanisms for pyroptosis. One is caspase‐1‐dependent typical inflammasomes, in which pro‐caspase‐1 was activated to caspase‐1 by the multiprotein complex; this phenomenon is associated with apoptosis‐associated speck‐like protein specks, NOD‐like receptors (NLRs), and non‐NLR (such as AIM2) that are switched on by immune stimulation such as pathogen‐associated molecular patterns, damage‐associated molecular patterns, or others. Another is caspase‐4/5‐dependent nontypical inflammasome, in which pro‐caspase‐4/5 was activated by cytosolic bacterial lipopolysaccharide released by gram‐negative bacteria. The pyroptosis pathway mediated by caspase‐3/GSDME is also associated with tumor chemotherapy. , In the pathway, caspase‐3 is spliced by Asp270, which converts cells that express GSDME noninflammatory apoptotic death into inflammatory pyroptotic death. , GSDMB or GSDMC might be spliced by caspase‐6 or caspase‐8 after tumor necrosis factor (TNF)‐mediated death receptor signaling to trigger pyroptosis, while the process of activated GSDMA is currently unknown. Activated caspase‐1/4/5 cleaves GSDMD into an N‐terminal GSDMD fragment, which binds to membrane phosphatidylinositol in eukaryotic cells and forms molecular pores similar to the pore‐forming toxins of bacteria. , Subsequently, cells swell and bubble, cell membranes dissolve, cellular contents and interleukins (including IL1β, IL18, caspase‐1/4/5, and other cytokines) are released, and an inflammatory reaction ensues. This is a typical feature of pyroptosis. Research had reported that pyroptosis was related to carcinogenesis and antitumor immunity. Abnormal expression of GSDMD positively correlates with the number and activity of CD8+ T lymphocytes. A recent study reported that 15% of cancer pyroptosis induced by GSDMA is sufficient to prompt antitumor immunity and clear 4T1 mammary tumors.

FIGURE 1

(A) Schematic diagram of pyroptosis mechanisms. Pyroptosis is mediated by the GSDM family, which usually cleaved to liberate an N‐terminal fragment by the activated inflammatory caspase family. (B) The workflow of this study. Abbreviations: TCGA: The Cancer Genome Atlas; CNV: copy number variation; DEGs: differentially expressed genes; PCA: principal component analysis; GSVA: gene set variation analysis; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; TMB: tumor mutation burden

In brief, pyroptosis plays an essential regulatory role in carcinogenesis and sensitivity to antitumor therapies. Recently, a study reported that reduced expression of calcium ion regulator CD38 prevents inflammasome‐induced pyroptosis in HNSCC. Cai et al. found that triptolide inhibited HNSCC cells by inducing GSDME‐mediated pyroptosis. Herein, we evaluated the characterization of pyroptosis patterns based on pyroptosis‐related genes (PRGs) in HNSCC, and we found that pyroptosis patterns were related to immune cell infiltration. Subsequently, we constructed a novel pyroptosis score algorithm using principal component analysis. We also assessed the functional enrichment, tumor mutation burden (TMB), immune cell infiltration, and the sensitivity to chemotherapy and immunotherapy between high and low pyroptosis score group.

MATERIALS AND METHODS

HNSCC datasets

The workflow chart of this study is shown in Figure 1B. The transcriptome data and corresponding clinicopathological characteristics of three HNSCC cohorts (TCGA‐HNSC, GSE41613, and GSE65858) from The Cancer Genome Atlas (TCGA) and the Gene‐Expression Omnibus (GEO) database were used. For the TCGA‐HNSC datasets, the fragments per kilobase of transcript per million (FPKM) values of gene expression, somatic mutation, and copy number variation (CNV) data were downloaded from the National Cancer Institute's Genomic Data Commons (https://portal.gdc.cancer.gov/). FPKM values were then transformed into transcripts per kilobase million values identical to those resulting from microarrays. Batch influences from nonbiological technical biases were rectified using the “ComBat” algorithm of the “sva” package. Finally, after removing patients without survival information, a total of 866 HNSCC patients were used for further analysis.

Multiomics landscape analysis of PRGs in HNSCC

Based on literature mining, , , , a total of 33 genes that were related to pyroptosis were identified (Table S1). The expression of 33 PRGs was extracted and compared between 502 HNSCC and 44 normal samples using the “limma” and “pheatmap” package from TCGA‐HNSC cohort. We identified 26 expressed PRGs in TCGA‐HNSC cohort. Because of the lack of expression data in two GEO cohorts, four genes (GSDME, NLRP6, PJVK, and SCAF11) were removed; the remaining 22 differentially expressed PRGs that included the CASP family (CASP1, CASP3, CASP4, CASP5, CASP6, and CASP8), the GSDM family (GSDMB, GSDMC, and GSDMD), the NLR family (NLRC4, NLRP1, NLRP3, and NLRP7), and other genes (AIM2, GPX4, IL1B, NOD1, NOD2, PLCG1, PYCARD, TNF, and ELANE) were selected for further evaluation. Mutation data of PRGs and TMB value of HNSCC samples were evaluated, summarized, and illustrated using the “maftools” package. CNV frequency of PRGs was also estimated. Univariate Cox regression analysis was performed to explore the prognostic value of PRGs using the “survival” and “survminer” packages. The interaction network for PRGs was constructed using |Pearson R| > 0.2 and p < 0.05 and the “igraph” and “reshape2” packages.

Consensus clustering analysis

Based on expression profiles of PRGs, consensus clustering analysis was performed to identify distinct pyroptosis signatures and classify HNSCC patients for subsequent evaluation. We increased the clustering variable (k) from 2 to 9 and repeated the operation 1000 times using the “ConsesusClusterPlus” packages to ensure the best stability of classification. Principal component analysis (PCA) was used to estimate two pyroptosis clusters in all HNSCC patients using the “Rtsne” and “ggplot2” R packages. Overall survival was compared between two pyroptosis clusters using Kaplan–Meier analysis. Using the Pearson's chi‐squared test, a heatmap was used to visualize the relationship between clinical factors and pyroptosis clusters.

Functional enrichment analysis and immune cell infiltration feature of pyroptosis clusters

To compare biological functions between two pyroptosis clusters, gene set variation analysis (GSVA) enrichment was performed using the “GSVA” R package. Characteristic gene data were obtained from the Molecular Signature Database. A single sample gene set enrichment analysis (ssGSEA) algorithm was used to evaluate the relative abundance of immune cell infiltration based on the 23 immune cell‐associated gene sets using the “GSVA” and “GSEABase” package.

Differentially expressed genes (DEGs) and pathway enrichment analysis of pyroptosis clusters

The “DESeq2” package filtered DEGs between two pyroptosis clusters according to specific criteria (|log2FC| ≥ 0.5 and p < 0.001) using the “limma” package. Based on these DEGs, Gene Ontology (GO) analysis was performed to identify enriched GO terms (BP: biological process; CC: cellular component; and MF: molecular function) using the “clusterProfiler” package with an adjusted P‐value of < 0.05. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was also performed to evaluate the related pathways between pyroptosis clusters.

Construction of pyroptosis‐related DEGs clusters and pyroptosis scores

First, consensus clustering analysis was employed to categorize the patients as the DEGs (‘gene clusters’). Kaplan–Meier survival curves were drawn, and the log‐rank tests were used to calculate differences among gene clusters. A heatmap was used to visualize the correlation patterns between gene clusters and the clinical variables using the “heatmap” package. The PCA method was then used to establish a pyroptosis score. PC1 and PC2 were brought in to act as signature scores. This algorithm focused the score on the largest module of well‐related (or nonrelated) genes in the set while downweighting contributions from genes that do not track with other set members. , , We calculated the pyroptosis score for each sample as follows:

pyroptosis score =, where i represents the level of DEGs.

Correlation between pyroptosis score and clinical subtypes

Cut‐point survival analysis was used to determine the optimal cut‐offs for dividing patients into low and high pyroptosis score groups. Sankey plots were used to display group changes of each patient, including pyroptosis clusters, gene clusters, pyroptosis score, and survival status. The Wilcoxon signed‐rank test was performed to compare the differences in pyroptosis score between the distinct two pyroptosis clusters or three gene clusters. Overall survival time was compared between high and low pyroptosis score groups using Kaplan–Meier analysis and the log‐rank tests. Subgroup survival analysis was used to assess the effectiveness of the pyroptosis score system in the age (<60 or ≥60 years), gender (female or male), and clinical‐stage subgroups (I‐II or III‐IV).

Correlation between pyroptosis score and TMB in HNSCC

Somatic mutation data were annotated based on the hg19 reference genome and were visualized using the “GenVisR” package. The “maftool” package was used to analyze the mutated genes in HNSCC, and the top 20 frequently mutated genes belonging to both high and low pyroptosis score groups were shown in waterfall plots. The total number of mutations counted was divided by the exome size (38 megabases) to calculate the TMB of each patient. Subsequently, “survival” and “survminer” packages were used to compare the overall survival based on TMB and pyroptosis score, visualizing using the Kaplan–Meier curves and log‐rank tests.

Immunotherapy and chemotherapy

The expression of ICI‐related genes (PD1, CTLA4, and PD‐L1) was compared between the high and low pyroptosis score groups. The immunophenoscore (IPS) of each HNSCC patient was collected from the TCIA database (https://tcia.at/home) to evaluate the response of immunotherapy, which was calculated using the expression of the immune‐related genes to denote four types of immune cells (immunosuppressive cells, MHC molecules, effector cells, and selected immunomodulators). The half‐maximum inhibitory concentration (IC50) of four conventional chemotherapeutic agents (paclitaxel, gemcitabine, docetaxel, and cisplatin) for HNSCC patients was counted using the “pRRophetic” package.

Statistical analysis

All statistical analyses and graphs were generated using R software (version 4.0.3). The Pearson chi‐squared and Wilcoxon signed‐rank tests were used for the categorical and continuous variables analyses. Pearson correlation analysis was used for correlation evaluations. p < 0.05 indicated statistical significance unless specified otherwise.

RESULTS

Multiomics landscape of PRGs in HNSCC

As shown in schematic diagram (Figure 1A), pyroptosis is inflammatory process involving multiple genes, which mediated the GSDM family cleaving by various activated caspases. We explored the characterization of these genes in HNSCC and designed the workflow for this study (Figure 1B). According to the flow chart, the expression of 33 PRGs was compared between 502 HNSCC and 44 normal tissues from TCGA. As shown in Figure 2A, a total of 25 genes were found to be upregulated, and ELANE was downregulated in tumor tissues (p < 0.05). Considering the lack of expression data in the two GEO cohorts, four genes (GSDME, NLRP6, PJVK, and SCAF11) were removed, and the remaining 22 differentially expressed PRGs (AIM2, CASP1, CASP3, CASP4, CASP5, CASP6, CASP8, GPX4, GSDMB, GSDMC, GSDMD, IL1B, NLRC4, NLRP1, NLRP3, NLRP7, NOD1, NOD2, PLCG1, PYCARD, TNF, and ELANE) were selected for further evaluation. To reveal the genetic variation landscape of PRGs, somatic mutations and CNV analysis were performed. The overall mutation rate of PRGs was relatively low in HNSCC patients, approximately 21.54%. CASP8 exhibited the highest mutation rate of 10%, while others showed lower mutations frequencies (0%–2%) in HNSCC samples (Figure 2B). CNV analysis revealed that GSDMC and GSDMD had relatively high amplification rates, whereas ELANE and GPX4 were mainly copy number deletions (Figure 2C). The chromosomal location of PRGs with CNV alteration is shown in Figure 2D. The comprehensive landscape of these PRGs interactions and prognostic significance for HNSCC patients were shown in the network chart (Figure 2E).

FIGURE 2

Multiomics landscape of PRGs in HNSCC. (A) The heatmap of 33 PRGs between HNSCC samples and normal tissues. *p < 0.05, **p < 0.01, and ***p < 0.001. (B) The mutation frequency of PRGs. Each column indicates individual patients. The upper bar plot shows TMB. (C) The CNV frequency of PRGs in TCGA‐HNSC cohort. (D) The location of PRGs with CNV mutation on chromosomes. (E) Interaction and outcome among PRGs in HNSCC patients. The size of each circle represents the prognostic value of the corresponding gene, p < 1e‐04, p < 0.001, p < 0.01, p < 0.05, and p < 1. The left half of each circle colored with red, gray‐blue, yellow, or gray represents the CASP family, GSDM family, NLR family, or other genes. The right half of each circle shows prognostic value, purple shows risk factors, and green shows favorable factors. The pink line represents a positive correlation with p < 0.0001, and the blue line means a negative correlation with p < 0.0001

Classification of HNSCC patients based on the PRGs

To explore the relationship of HNSCC characteristics and PRGs, consensus clustering analysis was performed in TCGA‐HNSC, GSE41613, and GSE65858. By adding clustering variables (k) from 2 to 9, we observed the highest intragroup and low intergroup correlations when k = 2, suggesting that these HNSCC samples could be classified into two clusters based on 22 differentially expressed PRGs; we called them pyroptosis clusters A (n = 420) and B (n = 446) (Figure 3A). PCA charts confirmed that two pyroptosis clusters were well distinguished (Figure 3B). Kaplan–Meier curves for overall survival of HNSCC patients in the two clusters were drawn; however, there was no statistical significance (p = 0.237, Figure 3C). The relationship between two pyroptosis clusters and clinicopathologic features including clinical stage (I‐II or III‐IV), gender (female or male), age (<60 or ≥60 years), and survival status (alive or dead) was illustrated in a heatmap. There were few differences in clinicopathologic features between the two pyroptosis clusters (Figure 3D).

FIGURE 3

Classification of HNSCC patients based on the PRGs. (A) HNSCC patients were grouped into two clusters according to the consensus clustering matrix for k = 2. (B) The PCA scatter plot indicating that two pyroptosis clusters were well distinguished. (C) Kaplan–Meier curves showing the overall survival of HNSCC patients between two pyroptosis clusters. (D) Heatmap showing the PRGs expression and clinicopathologic features of HNSCC patients between two pyroptosis clusters

Functional enrichment analysis and immune cell infiltration feature of pyroptosis clusters in HNSCC

To explore the biological functions of pyroptosis clusters, GSVA enrichment was performed. We observed that pyroptosis cluster A was markedly enriched in immune fully activation pathways, including the activation of the chemokine signaling pathway, the JAK‐STAT signaling pathway, cytokine‐cytokine receptor interactions, toll‐like receptor signaling pathways, and natural killer cell‐mediated cytotoxicity (Figure 4A). Relative abundance of immune cell analysis showed that pyroptosis cluster A was substantially richer in immune activation cells such as activated dendritic cells, NK cells, macrophages, eosinophils, mast cells, and plasmacytoid dendritic cells (Figure 4B). To further explore the heterogeneity of two pyroptosis clusters, we recognized the DEGs between the groups. 784 DEGs were identified, including 385 upregulated and 399 downregulated genes in pyroptosis cluster A (Table S2). GO enrichment (Figure 4C) and KEGG pathway analysis (Figure 4D) indicated that these DEGs were mainly related to immune responses, chemokine activity, and TNF signaling pathways. This finding suggests that pyroptosis might play a non‐negligible role in the immune regulation of HNSCC progression.

FIGURE 4

Functional enrichment analysis and immune cell infiltration feature of pyroptosis clusters in HNSCC. (A) Heatmap showing the biological pathways of two pyroptosis clusters by GSVA enrichment analysis. (B) The relative abundance of each immune cell of pyroptosis cluster A and cluster B. *p < 0.05, ***p < 0.001, ns: not significant. (C) GO functional annotation of the DEGs between two pyroptosis clusters. (D) KEGG pathways analysis of the DEGs between two pyroptosis clusters

Clinicopathologic features and overall survival among three gene clusters

To comprehensively annotate the molecular characteristics of pyroptosis patterns, consensus clustering analysis was presented based on 784 DEGs between two pyroptosis clusters and when k = 3 was identified with optimal clustering stability. All samples were regrouped into three subtypes (gene cluster A [n = 252], gene cluster B [n = 249], and gene cluster C [n = 365]) (Figure 5A). A PCA chart confirmed this finding (Figure 5B). We found that the overall survival of gene cluster A was significantly worse than the other two gene clusters (p < 0.001, Figure 5C). We then compared the expression of 22 PRGs among three gene clusters and found that the expression of all genes except NLRP7 were significantly different among three gene clusters (p < 0.001, Figure 5D). The clinicopathologic features, including clinical stage (I‐II or III‐IV), gender (female or male), age (<60 or ≥60 years), and survival status (alive or dead), were shown in a heatmap. There were few differences between the three gene clusters and the two pyroptosis clusters (Figure 5E).

FIGURE 5

Clinicopathologic features and overall survival among three gene clusters. (A) Consensus clustering matrix (k = 3) was used to classify HNSCC patients into three different genomic subtypes, gene cluster A, gene cluster B, and gene cluster C. (B) The PCA scatter plot indicating that three pyroptosis clusters were well distinguished. (C) Kaplan–Meier curves showing three gene clusters that are significantly related to overall survival in HNSCC patients (log‐rank, p < 0.001). (D) The differential expression of PRGs among three gene clusters. **p < 0.01, ***p < 0.001, ns: not significant. (E) Heatmap revealing the clinicopathologic features of HNSCC patients among three gene clusters and two pyroptosis clusters

Construction of pyroptosis score and overall survival for individual HNSCC patients

A novel pyroptosis score for each HNSCC sample was constructed using a PCA algorithm to evaluate the pyroptosis patterns in individual patients accurately. According to optimal cut‐offs of survival scores, all patients divided into high (n = 602) or low (n = 264) pyroptosis score groups. Sankey plots were used to visualize the group changes of individual patients, including pyroptosis clusters, gene clusters, pyroptosis scores, and survival status (Figure 6A). We compared the pyroptosis scores in two pyroptosis clusters and three gene clusters and found that pyroptosis cluster A showed lower pyroptosis scores than pyroptosis cluster B (Figure 6B). Gene cluster A represented a significantly lower pyroptosis score (Figure 6C). Kaplan–Meier survival plots and bar plots revealed that low pyroptosis score patients had shorter survival times (log‐rank, p < 0.001, Figure 6D) and higher death rates (Figure 6E, 54% vs. 36%) than those with high pyroptosis scores. It is reasonable that low pyroptosis scores significantly decreased in the deceased patients than survivors (p = 1.7e‐06, Figure 6F). Gene cluster A showed the worst overall survival and the lowest pyroptosis score, confirming the previous result (Figures 5C and 6C).

FIGURE 6

Construction of pyroptosis score and overall survival for individual HNSCC patients. (A) Sankey plot showing the relationship of pyroptosis clusters, gene clusters, pyroptosis score, and survival status. (B) Comparison of pyroptosis scores between two pyroptosis clusters. (C) Boxplot of pyroptosis scores for three gene clusters. (D) Kaplan–Meier curves indicating that the overall survival of the low pyroptosis score group was worse than the high score group (log‐rank, p < 0.001). (E) The survival proportion of patients in low and high pyroptosis score groups. (F) The correlation between pyroptosis score and survival status

Relationship between pyroptosis score and clinical subtypes in HNSCC patients

Kaplan–Meier curves were drawn to compare the survival time of HNSCC patients with different clinical subtypes in high and low pyroptosis scores. We found that patients with high pyroptosis scores presented significant clinical benefits and remarkably prolonged survival times in the age (<60 years, p < 0.001, Figure 7A; ≥60 years, p < 0.001, Figure 7B), male (p < 0.001, Figure 7D), and stage III‐IV(p < 0.001, Figure 7F) subgroups, while there was no significant difference in female (p = 0.096, Figure 7C) and stage I‐II (p = 0.968, Figure 7E) subgroups.

FIGURE 7

Relationship between pyroptosis score and clinical subtypes of HNSCC patients. Kaplan–Meier survival curve of clinical subtypes for age (<60 [A] and ≥60 [B]), gender (female [C] and male [D]), and clinical stage (I‐II [E] and III‐IV [F])

Correlation between the pyroptosis score and TMB in HNSCC

A waterfall plot showed the top 20 mutation genes in the low (Figure 8A) and high pyroptosis score groups (Figure 8B). Although there was no visible difference between low and high pyroptosis scores, TP53 had a high mutation frequency in the two groups, especially the low group (73%). Because of the critical role of TMB in clinical prognostic value, we calculated the TMB for each HNSCC patient and divided them into high and low TMB groups. Patients with low TMB had longer survival times significantly than high TMB (p = 0.004, Figure 8C). Considering the prognostic value of TMB and pyroptosis score in HNSCC, we performed a subgroup survival analysis based on TMB. In both high and low TMB score subgroups, the low pyroptosis group showed worse overall survival than the low pyroptosis group (p < 0.001, Figure 8D).

FIGURE 8

Correlation between pyroptosis score and TMB in HNSCC. Waterfall plot showing the top 20 mutation genes in low pyroptosis score group (A) and high pyroptosis score group (B). (C) Overall survival of the high TMB group was worse than the low TMB group (Log‐rank, p < 0.001). (D) Overall survival curves for HNSCC patients classified by pyroptosis score and TMB. “H” means high, and “L” means low level (log‐rank, p < 0.001)

The influence of pyroptosis score on immunotherapy and chemotherapeutic agents for HNSCC patients

We further analyzed the correlation between pyroptosis score and immune cells. As shown in Figure 9A, pyroptosis scores were significantly positively related to activated B cells, activated CD8T cells, and type 17T helper cells, while negatively related to CD56 bright NK cells, NK T cells, and regulatory T cells. Therefore, to explore the role of pyroptosis score in antitumor immunity in HNSCC patients, the correlation between pyroptosis score and immunotherapy, as well as chemotherapeutics, were evaluated. We found that the expression of PD‐L1 (p = 6.6e‐11, Figure 9D) was significantly lower in the high pyroptosis score than the low score group, suggesting that patients with low pyroptosis score might have a better response to ICIs. However, there is no significant difference in the expression of PD1 (Figure 9B) and CTLA4 (Figure 9C). The violin charts based on IPS confirmed that patients with low pyroptosis scores have better responses to PD1 inhibitor therapy alone (p = 0.00019, Figure 9E) or a combination with CTLA4 inhibitor and PD1 inhibitor (p = 0.0017, Figure 9G), while not for CTLA4 inhibitor alone (p = 0.38, Figure 9F). The IC50 of each HNSCC patient was measured based on the pRRophetic algorithm. There was a significantly lower IC50 in the low pyroptosis score group for paclitaxel (p = 0.0061, Figure 9H), docetaxel (p = 6.2e‐13, Figure 9J), and gemcitabine (p = 2.8e‐11, Figure 9K), whereas there was no difference between two pyroptosis score groups for cisplatin (p = 0.57, Figure 9I). These findings suggest that HNSCC patients with low pyroptosis scores are more sensitive to paclitaxel, docetaxel, and gemcitabine therapy.

FIGURE 9

Evaluation of the influence of pyroptosis score on immunotherapy and chemotherapeutics for HNSCC patients. (A) Correlation matrix of the pyroptosis score and immune cell proportions. Blue: negative correlation, red: positive correlation, *p < 0.05. The relationship between pyroptosis score and expression of PD1 (B), CTLA4 (C), and PD‐L1 (D). The correlation between pyroptosis score and immunophenoscore (IPS) for anti‐PD1 monotherapy (E), combination anti‐PD1 with anti‐CTLA4 immunotherapy (F), and anti‐CTLA4 monotherapy (G). Drug sensitivity analysis of paclitaxel (H), cisplatin (I), docetaxel (J), and gemcitabine (K) for HNSCC patients in high and low pyroptosis score groups

DISCUSSION

Several lines of evidence indicated that pyroptosis plays a crucial role in regulating antitumor immunity and inflammation. Latest two studies had reported the prognostic prediction of pyroptosis signature in HNSCC. , They identified both few PRGs with a 4‐gene and a 5‐gene signature, respectively, and analyzed the immune cell infiltration of the risk model simply. The influence of pyrolysis patterns on immunotherapy and chemotherapy is still unclear, Thus, a comprehensive analysis of the pyroptosis characteristics and immune infiltrates in HNSCC patients is needed to identify effective immunotherapeutic and chemotherapeutic strategies.

In this study, we constructed a novel pyroptosis score system based on DEGs between pyroptosis clusters instead of a few different expressed PRGs in HNSCC. Importantly, we focus on the correlation between pyroptosis score and the prospect of immunotherapy and chemotherapy. Firstly, we identified 22 differentially expressed PRGs from 33 PRGs and showed a comprehensive landscape of the somatic mutation, CNV frequency, chromosome location, and interactions in HNSCC. We then recognized two distinct pyroptosis clusters for HNSCC patients by consensus clustering analysis based on the 22 PRGs. Because of the highly proinflammatory features and the critical role of pyroptosis in antitumor immunity, , we performed GSVA analysis (which is usually used for predicting the variation in the pathway and biological functions activity in samples ). We found that pyroptosis cluster A was markedly enriched in immune full‐activation pathways such as the chemokine signaling pathway, cytokine‐cytokine receptor interactions, the toll‐like receptor signaling pathway, NK cell‐mediated cytotoxicity, and the NOD‐like receptor signaling pathway. These results were confirmed in subsequent comparison of immune cell infiltration in two pyroptosis clusters. We found that pyroptosis cluster A was significantly more infiltrated with activated dendritic cells, NK cells, macrophages, eosinophils, mast cells, and plasmacytoid dendritic cells than pyroptosis cluster B. Consistent with this finding, GO and KEGG enrichment showed that the DEGs between two pyroptosis clusters were related to inflammatory responses and activation of the immune system. However, there was no significant difference in clinic factors and survival time between the two clusters. To further evaluate the prognostic value of these PRGs, three gene clusters were identified as pyroptosis patterns for all HNSCC patients based on 784 DEGs. The overall survival result suggested that the three gene clusters had significantly distinct clinical outcomes.

To accurately evaluate the pyroptosis features in each HNSCC patient, we constructed a pyroptosis score system using the PCA algorithm. We found that patients with low pyroptosis scores usually presented shorter survival times and higher death rates, suggesting that low pyroptosis scores might predict poor outcomes in HNSCC. Consistently, the pyroptosis cluster A and gene cluster A had worse overall survival and lower pyroptosis scores than the other pyroptosis cluster and gene clusters. Further clinical subgroup features analysis based on age, gender, and clinical stage showed that the prognostic patterns were effective in all subgroups except the female and stage I‐II subgroups (possibly because of the low percent weight of case numbers), suggesting the universality of the prognostic patterns for pyroptosis score.

Study had reported a correlation between gene mutations and immunotherapy response or tolerance. Our findings suggest that a different mutation frequency existed in several genes between the high and low pyroptosis score patients, and some of these genes were related to sensitivity or resistance to immunotherapy. We also calculated overall survival in the TMB and pyroptosis score groups. We found that the low pyroptosis score group had worse overall survival regardless of the level of TMB, suggesting that the pyroptosis score system gave accurate individual predictive values for HNSCC patients. Several studies demonstrated that HNSCC is highly infiltrated with immune cells. , , Intriguingly, we found that pyroptosis score positively correlated with activated B cells, activated CD8T cells, and type 17T helper cells, suggesting that the high pyroptosis score patient have an activated immune system. Therefore, it is not difficult to understand why the high pyroptosis score group had a relatively long survival time.

Immunotherapy is becoming the fourth pillar in tumor treatment juxtaposed with surgery, radiotherapy, and chemotherapy; it aims to activate the immune system to kill tumor cells using the endogenous immune system. Immunotherapy can treat heterogeneous or disseminated tumors where targeted therapy is ineffective. Immunotherapeutic strategies, especially ICIs, have achieved breakthroughs in many tumor treatments. With the recent approval of pembrolizumab and nivolumab, ICIs for anti‐PD1/PD‐L1 for first‐line treatment in recurrent and metastatic patients in the United States, and the cure rates of HNSCC have advanced substantially. , Nevertheless, tumors can display inherent or acquired resistance to immunotherapy. Only 10%–20% of HNSCC patients treated with anti‐PD1/PD‐L1 antibodies showed favorable clinical outcomes. Therefore, exploring prognostic markers and assessing treatment sensitivity are essential for applying tumor immunotherapy. Understanding their interactions might improve immunotherapy because of the inflammatory characteristics of tumors and the role of antitumor immunity via pyroptosis. We found that the expression of PD‐L1 was significantly higher in low pyroptosis score group than high score group, suggesting that low pyroptosis score patients might obtain clinical effect after treatment with anti‐PD1/PD‐L1 antibodies. We then found that for immunotherapy with a PD1 inhibitor alone or combined with a CTLA4 inhibitor, HNSCC patients of low pyroptosis scores were more likely to benefit than high score patients. Therefore, we speculate that the pyroptosis score system might help develop an individualized and precise immunotherapy strategy.

The relationship between pyroptosis and cancer chemotherapy is complex. Several lines of evidence suggest that therapeutic strategies (chemotherapy, immunotherapy, targeted therapy, and radiotherapy) induce pyroptosis in tumor cells and then trigger a local or systemic antitumor response. However, the inflammatory reaction triggered by pyroptosis is a common side effect for normal cells. Zhang et al. reported that paclitaxel and cisplatin‐induced pyroptosis mediated by caspase‐3/GSDME in lung cancer. A polo‐like kinase1 (PLK1) inhibitor‐induced pyroptosis via a caspase‐3/GSDME‐mediated mechanism in esophageal squamous cell carcinoma significantly enhanced cisplatin sensitivity in vitro and in vivo. Lu et al. observed that lung cancer cells exhibited pyrolysis after treatment with tyrosine kinase inhibitors (TKIs), including trametinib, erlotinib, and ceritinib. These studies found that pyroptosis promoted the antitumor effect of TKIs. In the present study, we analyzed the role of pyroptosis score in predicting chemosensitivity of HNSCC. The low pyroptosis score group had a markedly higher sensitivity than the high score group in treatment of paclitaxel, docetaxel, and gemcitabine. This finding suggests that the pyroptosis caused by chemotherapeutic drugs was more highly activated in the low score group, and the pyroptosis score might serve as a potential predictive marker for individual HNSCC patients before chemotherapy.

This study also has some limitations, including the lack of mechanism of the effects of pyroptosis patterns on immune infiltration and chemotherapy in HNSCC. Our conclusions would be more reliable if there were an additional validation set. Therefore, large‐scale, multicenter studies are needed to confirm our findings.

In conclusion, we successfully constructed a novel pyroptosis score algorithm by PCA analysis, and we found that HNSCC patients with low pyroptosis scores had better immunotherapeutic responses and higher sensitivities to chemotherapeutic agents. We speculate that the pyroptosis score might predict immunotherapeutic responses and chemosensitivity in individual HNSCC patient, and this would assist oncologists in generating accurate and individualized strategies.

CONFLICT OF INTEREST

The authors declare that they have no competing conflicts of interest.

AUTHOR CONTRIBUTIONS

Chongchang Zhou and Hongxia Deng conceived and designed this study. Hongxia Deng and Yangli Jin drafted the manuscript. Chongchang Zhou analyzed and expressed the data. Zhengyu Wei, Shijie Qiu, Dong Ye, Shanshan Gu, Yi Shen, and Zhisen Shen took part in data analysis. All authors reviewed and approved the final manuscript.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.