Prognosis and Characterization of Immune Microenvironment in Head and Neck Squamous Cell Carcinoma through a Pyroptosis-Related Signature

Pyroptosis, as a novel identified programmed cell death, is closely correlated with tumor immunity and shows potential roles in cancer treatment. Discerning a pyroptosis-related gene signature and its correlations with tumor immune microenvironment is critical in head and neck squamous cell carcinoma (HNSCC). Transcriptome data and corresponding clinical data were downloaded from TCGA and GEO databases. Tumor mutation burden (TMB) data were obtained from TCGA database. Firstly, univariate and least absolute shrinkage and selection operator (LASSO) regression analyses were used to construct a six pyroptosis-related gene signature. Kaplan–Meier analysis, receiver operating characteristic (ROC) curves, and principal component analysis (PCA) results verified that the risk model has good performance in predicting the survival. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses revealed that the pyroptosis-related gene signature was immune related. Finally, the immune landscape and immunotherapy sensitivity prediction capabilities of the risk model were further explored. There were close correlations between the overall survival (OS) and various immune cells and immune functions. Single-sample gene set enrichment analysis (ssGSEA) showed that high risk group had decreased expression of various immune cells and lower activities of immune functions. Meanwhile, tumor mutation burden (TMB) data combining risk score could well predict the OS of HNSCC patients. However, tumor immune dysfunction and exclusion (TIDE) analysis revealed that there was no significant difference in the sensitivity to immunotherapies between high and low risk groups. Finally, a nomogram based on risk score and clinicopathological parameters was constructed. And, the risk model demonstrated better sensitivity and specificity than TIDE scores and T-cell-inflamed signature (TIS). In conclusion, although the risk model could not well predict the immune escape and response to immunotherapies, the signature established by pyroptosis-related genes, with better sensitivity and specificity than TIDE scores and TIS signature, could be used for predicting prognosis and immune status of HNSCC patients.


Introduction
Head and neck cancer, including tumors originated from oral cavity, nasopharynx, oropharynx, hypopharynx, larynx, and tongue, accounts for more than 830,000 newly diagnosed cases worldwide [1]. Approximately 90% of head and neck cancers are head and neck squamous cell carcinoma (HNSCC). e primary treatment strategies for HNSCC include surgery, radiotherapy, and chemotherapy. e 5-year survival rate is approximately 50-60% due to the high heterogeneity of HNSCC [2,3]. Recently, immune checkpoint inhibitors (ICIs), such as anti-PD1/PD-L1 and anticytotoxic T lymphocyte antigen 4 (CTLA-4), have been applied in tumor treatment [4]. However, only a small population of patients could respond to PD-1/PD-L1 antibodies [5,6]. us, identifying biomarkers that could be used as prognostic factors and treatment targets is an urgent need in HNSCC. Pyroptosis, also known as caspase 1-dependent programmed cell death, has gained a lot of attentions recently. Pyroptosis could be triggered by microbial infection and noninfectious stimuli, such as cancers [7]. Pyroptosis is morphologically and mechanistically distinct from other forms of cell death, which is characterized by cell swelling, rapid plasma membrane rupture, and release of proinflammatory cytokines [8,9]. e biological functions of pyroptosis in the tumor process are still controversial. On the one hand, as a lytic form of programmed cell death, pyroptosis could mediate tumor cell killing [10,11]. On the other hand, as a form of proinflammatory death, pyroptosis could promote tumor growth by forming a suitable tumor growth microenvironment [12]. Multiple pyroptosis-related genes have been identified in various malignancies, including head and neck cancer. GSDME-mediated pyroptosis promoted the cytotoxicity of triptolide against head and neck cancer cells [13]. NLRP3 inflammasome was involved in pyroptosis activation of HNSCC cells [14]. However, the roles of pyroptosis related genes in predicting the survival of HNSCC patients are still largely unknown.
Chemotherapy and targeted therapy could eliminate tumor cells by inducing the pyroptosis of tumor cells [15,16]. And induction of pyroptosis has been viewed as a potential cancer treatment strategy [13]. Moreover, pyroptosis could mediate tumor regression by initiating antitumor immunity [17]. However, the therapeutic roles of pyroptosis in HSCC still need to be elucidated.
In this study, we firstly establish a pyroptosis gene signature for the prediction of OS in HSNCC. Secondly, the performance and biological functions of the risk model were further investigated. Moreover, the associations with immune microenvironment and immunotherapy responses were explored. Finally, a nomogram was established to predict the OS of patients with HNSCC.

Data Download and
Processing. RNA-sequencing data and corresponding clinical data were downloaded from the TCGA database (https://tcga-data.nci.nih.gov/tcga/) and the NCBI Gene Expression Omnibus (GEO) database. e accession number of GEO data was GSE65858, and the data platform number was GPL10558. TCGA-HNSCC cohort included 44 normal tissues and 502 HNSCC tumor tissues. GSE65858 contained RNA-sequencing data and complete clinical information of 270 HNSCC patients. All data were obtained from online public sources. us, ethical approval was not required.

Identification of Pyroptosis-Related Genes.
A list of 52 pyroptosis genes was prepared by screening published literatures [17][18][19][20]. e expression of pyroptosis genes was retrieved using "limma" package in R software. Differentially expressed pyroptosis-related genes with significant differences were compared between normal tissues and HNSCC tumor tissues with P < 0.05. Differentially expressed pyroptosis genes were visualized using a heatmap by utilizing "pheatmap" package in R software.
Protein-protein interaction (PPI) network of these differentially expressed pyroptosis-related genes was constructed using Search Tool for Recurring Instances of Neighbouring Genes (STRING) database (https://stringpreview.org/). e correlations between these differentially expressed pyroptosis genes were visualized using "igraph" and "reshape2" packages in R software.

Consensus Clustering Analysis.
Patients were classified into different groups based on the best k-means clustering by using the ConsensusClusterPlus R package. Cluster consistency clustering analysis was based on the expression of 6 pyroptosis-associated genes by setting the cluster variables (k) set from 2 to 10.
Kaplan-Meier survival analysis between different clusters was performed by using "survival" and "survminer" packages in R software. e correlations of different clusters with clinical parameters (Age, Gender, Grade, Stage, and T, N, and M stages) were assessed and demonstrated by a heatmap using "pheatmap" package.

Univariate Cox Regression Analysis.
Pyroptosis genes with prognostic values were identified using univariate cox regression analysis with P < 0.05. Least absolute shrinkage and selection operator (LASSO) regression analysis was used to construct a risk model by using "glmnet" and "survival" packages.

Performance of the Risk Model.
Firstly, risk score was calculated by the following formula: expression value of gene1 * coefficient of gene1+. . .+Expression of gen-eN * coefficient of geneN. Coefficient of each gene was generated from multivariate Cox regression. HNSCC patients were divided into high risk and low risk group by using the median risk score as a cutoff. e performance of the risk model was evaluated using Kaplan-Meier survival analysis and the area under curve (AUC) of the receiver operating characteristic (ROC) curves. Kaplan-Meier survival analysis was used to compare the survival rate between high risk and low risk groups by using "survival" and "survminer" packages in R software. And the ROC curves were drawn by using "survival," "survminer" and "timeROC" packages. Moreover, the distributions of risk scores and survival statuses were plotted by "pheatmap" package.

Principal Component Analysis (PCA)
. PCA with nonlinear methods was used to detect sample-to-sample heterogeneity by using "Rtsne" and "ggplot2" packages. e performances of the risk model were verified in both TCGA and GEO databases.

Independent Prognostic Value of the Risk Model.
Univariate and multivariate Cox regression analyses were performed to determine the independent prognostic values of risk score and clinical factors (Age, Gender, Grade, Stage, and T, N, and M stages) using "survival" package. Clinical data were retrieved from both TCGA and GEO databases. Moreover, a heatmap was used to evaluate the differences if prognostic pyroptosis genes and risk scores in clinicopathological parameters by using "pheatmap" R package.

Gene Enrichment and Functional Annotation Analysis.
Firstly, differentially expressed genes between high risk and low risk groups were identified by using "limma" package. Gene functional enrichment analyses were explored using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. GO enrichment analysis was carried out using "clusterProfiler", "org.Hs.eg.db," "enrichplot," and "ggplot2" R packages. KEGG analysis was developed using "clusterProfiler", "org.Hs.eg.db," "enrichplot," and "ggplot2" R packages.

Immune Cells' Infiltration and Immune Function.
CIBERSORT algorithm was used to estimate the abundances of various immune cell subsets in each sample. ESTIMATE method was used to calculate the immune scores, stromal scores, and estimate scores in TCGA-HNSCC cohort. Kaplan-Meier survival analysis was used to predict the prognostic values of various immune cells and immune functions. Immune cells' infiltration in different risk groups was explored using the single-sample Gene Set Enrichment Analysis (ssGSEA) based on the gene expression in TCGA-HNSCC cohort and GEO database. ere were 16 types of immune cells and 13 types of immune functions being compared between high risk and low risk groups using "reshape" and "ggpubr" packages.

Calculation of the Tumor Mutation Burden.
e somatic mutation data of HNSCC patients was downloaded and retrieved from the TCGA database. HNSCC patients were divided into high TMB and low TMB groups based on the median counts of TMB. Kaplan-Meier survival analysis was used to compare the survival rate between different groups by using "survival" and "survminer" packages in R software.

Tumor Immune Dysfunction and Exclusion (TIDE)
Analysis. Data of tumor immune dysfunction and exclusion (TIDE) was downloaded from http://tide.dfci.harvard.edu/. TIDE could be used to validate the performance on predicting anti-PD1 and anti-CTLA4 responses. e correlations of the risk model with the dysfunction of tumor infiltrating cytotoxic T cells, exclusion of cytotoxic T cells by immunosuppressive factors, microsatellite instability (MSI), and TIDE prediction scores were assessed using "limma" and "ggpubr" packages in R.

Construction and Validation of a Nomogram Model.
A nomogram model was constructed to comprehensively assess the survival probability of HNSCC patients, incorporating clinical factors (Age, Gender, Grade, Stage, and T, N, and M stages) and risk score. e predictive probabilities for 1-, 3-, and 5-year clinical outcomes were depicted by calibration curves. e sensitivity and specificity of the nomogram model in predicting 1-, 3-, and 5-year survival were determined using ROC curves. Moreover, the performance of the risk model in predicting OS was compared with TIDE scores and TIS signature using ROC curves.

Statistical Analysis.
All the data procession and statistical analyses were performed by using Perl software and R programming language (version 4.1.0). In all data analyses, a two-tailed P < 0.05 was considered statistically significant.

Identification of Pyroptosis-Related Genes.
Firstly, the expression of 52 pyroptosis-related genes were retrieved from the TCGA-HNSCC cohort. e expressions of these 52 pyroptosis-related genes were compared between 44 normal tissues and 502 HNSCC tumor tissues. As shown in Figure 1(a), there were 38 pyroptosis-related genes with significant differences between normal tissues and tumor tissues (P < 0.05). Seven pyroptosis-related genes (ELANE, IL18, CHMP2B, CHMP4C, CASP9, CHMP6, and CHMP2A) were downregulated in tumor tissues, while 31 pyroptosis-related genes were upregulated in tumor tissues (Figure 1(a)). e interactions of these pyroptosis-related genes were analyzed by protein-protein interaction (PPI) network according to STRING database (Figure 1(b)). Moreover, the coexpressed network of these pyroptosis genes are visualized in Figure 1(c).

Patients Classification and Correlations with Clinical
Parameters. HNSCC patients were classified using the consensus clustering analysis based on the differentially expressed pyroptosis-related genes. As shown in Supplementary Figure 1(a), when the k-value was 3, the HNSCC patients could be well classified. Kaplan-Meier survival analysis revealed that there were significant differences in the overall survival of HNSCC patients (P � 0.002, Supplementary Figure 1(b)). Cluster 1 demonstrated better OS than cluster 2 and cluster 3, while cluster 3 exhibited poorer OS than cluster 1 and cluster 2. Moreover, the correlations between the gene expression profiles and variable clinical parameters (Age, Gender, Grade, Stage, and T, N, and M stages) were further studied. e results were depicted in a heatmap, where there were significant differences in Grade (P < 0.05) and Gender (P < 0.05) between different clusters (Supplementary Figure 1(c)).

Validation of the Pyroptosis-Related Gene Prognostic
Signature. HNSCC patients were divided into high risk and low risk groups based on the median value of the risk scores. Survival data were retrieved from both TCGA and GEO databases. TCGA-HNSCC data showed that patients in high risk group demonstrated poorer overall survival than low risk group (Figure 3

Prognostic
Significance of the Risk Model. Moreover, the prognostic significance of the risk score was assessed using univariate and multivariate Cox regression analyses. Data retrieved from the TCGA database showed that Age (P � 0.007), T (P � 0.010), N (P � 0.001), and risk score (P < 0.001) were independent prognostic factors for HNSCC patients (Figures 4(a) and 4(b)). Data retrieved from GEO database also verified that Age (P � 0.012), T (P � 0.022), N (P � 0.034), and risk score (P � 0.015) were independent prognostic factors (Figures 4(c) and 4(d)).
Moreover, a heatmap depicted the correlations between the pyroptosis-related genes and risk scores. As shown in Figure 4(e), BAK1, GSDME, and IL1A were high-risk genes, while CHMP7, GZMB, and NLRP1 were low-risk genes. Moreover, there were significant differences in Grade between high risk and low risk groups.

Functional Enrichment Analysis of the Pyroptosis Gene
Signature. We next explored the biological functions and signaling pathways underlying the risk model by using Gene Ontology (GO) enrichment analysis and Kyoto.
Encyclopedia of Genes and Genomes (KEGG) Pathway Analyses. To begin with these studies, we firstly extracted 20 differentially expressed genes between high risk and low risk groups based on the TCGA database. GO enrichment analysis consisted of biological processes (BPs), molecular functions (MFs), and cellular components (CCs). GO results revealed that the functions of risk model were immune related ( Figure 5(a)). KEGG results showed that the underlying signaling pathways included cytokine-cytokine receptor interaction and chemokine signaling pathway, and so on ( Figure 5(b)).

Immune Landscape of HNSCC.
e immune landscape of HNSCC was analyzed using CIBERSORT algorithm. As shown in Figure 6(a), a barplot showed the infiltration of various immune cells. Next, the correlations between infiltrating immune cells and overall survival of HNSCC patients were analyzed by Kaplan-Meier survival curve. e high infiltration of aDC ( Figure 6

Immune Significance of the Pyroptosis Gene Signature.
Single-sample gene set enrichment analysis (ssGSEA) was employed to assess tumor immune-cells' infiltration and immune functions between high risk and low risk groups in both TCGA and GEO databases. In the TCGA-HNSCC cohort, high risk group demonstrated generally lower levels of 16 types of immune cells (Figure 8(a)). Meanwhile, high risk group also demonstrated lower activity of 10 kinds of immune functions, such as APC coinhibition, APC costimulation, CCR (chemokine and chemokine receptor), checkpoint, and so on (Figure 8(b)). Similar results were generated from the GEO database. e infiltration of 16 types of immune cells and 11 kinds of immune functions was lower in high risk group than that in low risk group (Figures 8(c) and 8(d)).
Moreover, the correlations between the prognostic pyroptosis-related genes and various immune cells are summarized in Figure 8(e). Here, red color represents positive correlation, while blue color represents negative correlation.

Significance of Pyroptosis-Related Gene Signature in
Immunotherapy. We next studied whether pyroptosis-related gene signature could contribute to clinical benefit of immune checkpoint inhibitor treatment. TIDE analysis showed that there were no significant differences in dysfunction of tumorinfiltrating cytotoxic Tcells (A), exclusion of cytotoxic Tcells by immunosuppressive factors (B), microsatellite instability (MSI) (C), and TIDE prediction scores (D) between high and low risk groups (Supplementary Figure 2).  the median value of TMB. Kaplan-Meier (K-M) survival analysis with log-rank tests indicated that the low-TMB group demonstrated better OS than high-TMB group (Figure 9(a), P � 0.004). We next combined the risk score and TMB to predict the survival of HNSCC patients. As shown in Figure 9(b), high TMB and high risk predicted the worst OS, while low-TMB and low risk group demonstrated the best prognosis. ese results indicated that risk score and TMB could be utilized simultaneously for predicting patient prognoses.  Journal of Oncology 13

Construction of a Prediction
( Figure 10(a)).
e performance of the nomogram was assessed using AUC index of ROC analysis and calibration curve. e AUCs of nomogram model for predicting 1-, 3-, and 5-year overall survival rates were 0.677, 0.738, and 0.744, respectively (Figure 10(b)). Calibration curve for the probability of 1-, 3-, and 5-year OS demonstrated good consistency between the actual observed survival and the nomogrambased prediction (Figure 10(c)). What is more, the predictive capabilities of risk score, TIDE, and TIS were 0.719, 0.509, and 0.481, respectively, indicating that the pyroptosis-related gene signature has better sensitivity and specificity in predicting the OS of HNSCC patients (Figure 10(d)).

Discussion
Pyroptosis was reported to be involved in the tumor progression and demonstrated great potential in cancer treatment [23]. As a form of programmed necrotic cell death,   pyroptosis was like a double-edged sword in malignancies [10][11][12]. On the one hand, pyroptosis could kill the tumor cells and mediate tumor regression [10,11]. On the other hand, pyroptosis could modify the tumor microenvironment to promote tumor growth [12]. e effects of pyroptosis on the biological functions of tumor cells depended on various pyroptosis-related genes. Moreover, the expression of pyroptosis-related genes were correlated with the survival of some cancer patients, including ovarian cancer [24], lung adenocarcinoma [25], skin cutaneous melanoma [20], gastric cancer [26], and laryngeal squamous cell carcinoma [27]. However, the expression and prognostic values of pyroptosis-related genes in HNSCC still need to be comprehensively elucidated.

T cells regulatory (Tregs) T cells gamma delta T cells follicular helper T cells CD8 T cells CD4 naive T cells CD4 memory resting T cells CD4 memory activated
Here, a list of pyroptosis-related genes was retrieved by literature searching, which includes 52 pyroptosis-related genes. And thirty-eight out of 52 pyroptosis-related genes were differentially expressed between normal tissues and HNSCC tumor tissues. ree clusters, divided by consensus clustering analysis, demonstrated significant differences in the OS of HNSCC patients, as well as other clinical factors including gender and grade. What is more, univariate Cox regression analysis and LASSO Cox analysis were performed to construct a pyroptosis-related gene prognostic signature. e risk model was consisted of 6 pyroptosis-related genes, including BAK1, GSDME, IL1A, CHMP7, GZMB, and NLRP1. Most of these six prognostic pyroptosis-related genes were once reported to be correlated with the progression and prognosis of cancer patients. BAX1, also known as autophagy-related gene and apoptosis signaling molecule, could predict the survival of head and neck cancer [28][29][30]. Pan-cancer analysis revealed that GSDME was highly expressed in most malignancies and significantly correlated with patients' survival [31]. Ibrahim et al. identified GSDME as a promising biomarker in the detection of colorectal cancer [32]. IL1A, as a proinflammatory cytokine, was upregulated in many types of cancers, such as breast cancer and lung cancer [33,34]. Induction of IL1A could promote proliferation, angiogenesis, and metastasis of tumor cells [33]. You et al. found that abnormal IL1A expression was correlated with poor prognosis of triple negative breast cancer [35]. CHMP7, as a ESCRT-III-related protein, was related to nuclear envelope reformation [36]. CHMP7 was identified as a prognostic biomarker in gastric cancer [37]. Granzyme B (GZMB), as an inflammatory gene and cytotoxic gene, participated in immune response and tumor cell killing. GZMB was identified as a progression biomarker in basal-like breast cancer [38]. NLR family, pyrin domain containing 1 (NLRPL) could form inflammasome and was involved in immune responses and cell death [39]. Low expression of NLRP1 was correlated with poor overall survival in lung adenocarcinoma [25,40]. e expression and prognostic significance of most of these prognostic pyroptosis-related genes have never been reported in HNSCC. In this study, all these genes were significantly correlated with the OS of HNSCC patients. Multiple genes risk model usually performed better sensitivity and specificity than single gene in predicting OS, as single gene usually showed lower sensitivity and specificity based on ROC analysis. Here, our risk model demonstrated good performance according to data generated from TCGA and GEO database. High risk group had a significantly poorer overall survival than the low risk group. Univariate and multivariate Cox regression analyses further verified that the risk score could be used as an independent prognostic factor for HNSCC patients.
We next explored the potential functional enrichment of the pyroptosis-related gene signature. GO and KEGG analyses indicated that the risk model was mainly involved in immune responses, cytokine activity, and chemokine and chemokine receptor activity. ese results strongly indicated that pyroptosis exhibited a significant correlation with immune microenvironment. CIBERSORT estimation results showed the abundance of immune cells infiltration and immune functions in the tumor microenvironment of HNSCC. Indeed, previous studies have indicated that there were close correlations between pyroptosis and immune response and immune cells infiltration. Pyroptosis could stimulate effective immune responses. Wang et al. found that nanoparticle-conjugated gasdermin could sensitize 4T1 tumor cells to anti-PD-1 therapy strategy [16]. e expression of Gasdermin E (GSDME) was linked with the cancer-associated fibroblast infiltration in multiple malignancies [31]. Increased expression of GZMB was connected with high numbers of mutations in non-small-cell lung cancer [41]. Diminished expression of GZMB was associated with higher PD-1, PD-L1, LAG-3, and TIM-3 expression in invasive bladder cancer [42]. e expression of NLRP1 was positively correlated with the degree of tumor-infiltrating immune cells in lung adenocarcinoma. Accordingly, both TCGA and GEO data revealed that there were significant differences in the immune cells infiltration and immune functions between the high risk and low risk groups. High risk group demonstrated decreased expression of immune cells and repressed activities of immune functions. All the prognostic pyroptosis-related genes were significantly correlated with the infiltration of both positive and negative immune cells, such as Tregs, macrophage M0, macrophage M1, macrophage M2, CD8+T cells, and so on. All these data strongly indicated that our pyroptosis related risk model was associated with the tumor immune microenvironment. However, there were no significant differences in tumor immune dysfunction and exclusion between high and low risk groups. Even ICIs targeting PD-1, such as Nivolumab and Pembrolizumab, have been approved in the treatment of HNSCC [4,43]. About 70% of HNSCC patients failed to benefit from ICIs [43]. Our data partially explained the low response rate to ICI in HNSCC. What is more, the immune-related mechanisms underlying pyroptosis might not be PD-1/PD-L1 or CTLA-4 related.
ere must be other immune-related mechanisms that needed to be further elucidated. Finally, a novel HNSCC nomogram was constructed according to the risk score and clinical parameters. AUC values and calibration curves demonstrated good predictive ability. Meanwhile, the sensitivity and specificity of the risk model were much better than TIDE score and TIS signature.
Intensive studies have shown that HNSCC exhibited enormous tumor heterogeneity, leading to different clinicopathological characters and therapeutic responses based on different tumor regions, distinct mutational profiles, and variable molecular characters [44]. e main limitation of this study was that we constructed the pyroptosis gene signature to predict OS and immune response incorporated all HNSCC patients without considering tumor heterogeneity. It is worthwhile to further explore the relationships between pyroptosis genes and various phenotypes of HNSCC, exploring their potential clinical applications to target different types of HNSCC patients.

Conclusion
In conclusion, our study provided a novel model for predicting the survival of HNSCC. is model showed good correlations with tumor immune microenvironment.