Relationships of Ferroptosis and Pyroptosis-Related Genes with Clinical Prognosis and Tumor Immune Microenvironment in Head and Neck Squamous Cell Carcinoma

Ferroptosis and pyroptosis are two new programmed cell death (PCD) modes discovered in recent years. However, the potential value of ferroptosis and pyroptosis-related genes (FPRGs) in prognosis prediction and the tumor immune microenvironment of head and neck squamous cell carcinoma (HNSCC) is still unclear. We obtained 21 significant FPRGs based on the training dataset (TCGA- HNSC) using the univariate Cox and differential expression analysis. The TCGA- HNSC (n = 502) dataset was clustered into two group (clusters A and B) based on the 21 significant FPRGs. 1467 differentially expressed genes (DEGs) between cluster A and B were put into univariate Cox and Least absolute shrinkage and selection operator (LASSO) analysis to build a risk model. The predictive capability of the risk model was successfully confirmed by internal validation, external validation, and clinical sample validation. To improve the clinical applicability, a nomogram model combined risk score and clinical information were constructed. Moreover, the patients with lower risk score were characterized by increased immune response and tumor mutation burden (TMB), while the patients with higher risk score were characterized by increased TP53 mutation rate. In conclusion, our comprehensive analysis of the FPRGs revealed their significant role in prognosis prediction and the tumor immune microenvironment. The risk model containing 9 FPRGs could be a potential prognostic markers and effective immunotherapy targets for HNSCC.


Background
Cell death is closely related to the basic processes of life and is an important way of growth and development, disease progression, and homeostasis of multicellular organisms [1]. In 2015, the Nomenclature Committee on Cell Death (NCCD) classified cell death into programmed cell death (PCD) and un-programmed cell death according to whether the process of death was regulated by procedures [2]. PCD is an active and orderly way of cell death to maintain the stability of internal environment. Specifically, it refers to the suicide protection measures initiated by gene regulation when cells are stimulated by internal and external environmental factors [3]. PCD contains apoptosis, necroptosis, autophagy, pyroptosis, ferroptosis, and other cell death modes, which play an important role in pathogen immunity and cancer cell clearance [4]. Among them, ferroptosis and pyroptosis are two new PCD modes discovered in recent years [5]. Cell pyroptosis is a novel pro-inflammatory programmed cell death mode, which depends on the activation of cysteinyl aspartate specific proteinase (caspase) and its mediated gasdermin D (GSDMD). Activated caspase mediates the hydrolysis of GSDMD into bioactive GSDMD-N, which embedded into the cytoplasmic membrane to form membrane perforation with a diameter of 10~15 nm. This causes increased cell permeability, imbalance of ion compensation, and water inflow into cells from the intercellular substance, resulting in cell swelling, and large release of lactate dehydrogenase and pro-inflammatory cytokines, such as IL-1β and IL-18 [6,7]. Ferroptosis is an irondependent regulatory form of cell death, including activation of reactive oxygen species (ROS), iron aggregation, activation of the mitogen-activated protein kinase (MAPK) system, reduced cysteine uptake, and glutathione depletion [8]. Ferroptosis is characterized by excessive accumulation of iron-dependent lipid peroxidation (LPO) on cell membranes, leading to cell necrosis, which can be inhibited by glutathione peroxidase 4 (Gpx4). Morphology of ferroptosis showed a loss of membrane integrity, accompanied by nuclear swelling and mitochondrial shrinkage, increased membrane density and mitochondrial outer membrane rupture [9].
Head and neck squamous cell carcinoma (HNSCC) is an immunosuppressive disease characterized by molecular heterogeneity and tumor-host interaction. Its morbidity and mortality are increasing year by year, and it is now the sixth most common cancer and the eighth leading cause of cancer death worldwide [10]. The treatment of early HNSCC is mainly surgery and radiotherapy, but the 5-year survival rate is less than 40% because most of patients have locally advanced disease at first diagnosis [11]. Platinum-based chemotherapy for advanced HNSCC has a poor prognosis, with a median survival less than 1 year [12]. Therefore, the treatment of HNSCC is in urgent need of new drug breakthrough. Exogenous activation of ferroptosis and pyroptosis has recently been shown to trigger powerful antitumor effects [13]. Some chemotherapy drugs can switch caspase-3-mediated apoptosis to pyroptosis by cleaving GSDME into GSDME-N in GSDME-expressing tumor cells [14,15]. GSDME, as a tumor suppressor, can improve the antitumor immunity by activation of pyroptosis, while inflammasome activation induced by pyroptosis further enhances the therapeutic efficacy of some immune checkpoint blockers [16,17]. Ferroptosis can improve the cytotoxicity of cisplatin of resistant HNC cells [18] and the efficacy of radiotherapy [19]. Ferroptosis can be induced by sorafenib, a kinase inhibitor, which has been reported that it can increase the radiosensitivity and antiproliferative effect of cisplatin in HNSCC cells [13,20]. Therefore, induction of ferroptosis and pyroptosis may provide an effective treatment strategy of HNSCC.
In this study, we systematically investigated the role of ferroptosis and pyroptosis-related genes in prognosis prediction and the tumor immune landscape of HNSCC. We first constructed and validated a risk model based on the ferroptosis and pyroptosis-related genes (FPRGs). The HNSCC patients were clustered into high-and low-risk group based on the median cut-off of risk score. Then, we assessed the clinical features, tumor mutation burden (TMB), cancer stem cell (CSC) characteristics, and immune infiltration in the two groups. This study paves a novel road for prognosis prediction and treatment strategy of HNSCC.

Materials and Methods
2.1. Data Acquisition. The workflow of this study was shown in Figure 1. The mRNAs-seq data, somatic mutation data, copy number variation (CNV) data, and corresponding clinical information of TCGA-HNSC dataset including 44 normal samples and 502 HNSCC samples were downloaded from The Cancer Genome Atlas (TCGA) database. The mRNAs-seq data and clinical information of GSE65858 dataset (270 HNSCC samples) based on the platform GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip, Illumina Inc., San Diego, CA, USA) were obtained from the Gene Expression Omnibus (GEO) database. It was generated from samples of peripheral blood mononuclear cells (PBMCs) of patients. The "limma" package was used to normalized the expression profiles data. The baseline information is shown in table 1.

Unsupervised Clustering for FPRGs. We downloaded 313
FPRGs from the predecessors' study [21] and extracted the mRNAs-seq data of the 313 genes from the TCGA-HNSC dataset. The differentially expressed (DE) FPRGs between normal and HNSCC samples were screened using the "limma" package based on the selection criteria of log jFCj ≥ 1 and p < 0:05 [5]. The FPRGs with prognostic value (p < 0:05) was selected by univariate Cox regression analysis. Then, the DE FPRGs with prognostic value were subjected to consensus clustering algorithm. The "ConsensuCluster-Plus" package was performed with 1000 times repetitions to guarantee the stability of classification [22].

GSVA and GO Functional Enrichment Analysis.
Gene set variation analysis (GSVA) was used to investigate the differentially activity of molecular pathways between different subtypes using the "GSVA" packages in R software [23]. The gene file of "c2.cp.kegg.v7.4.symbols.gmt" was downloaded from MSigDB database for GSVA analysis, and p < 0:05 was considered as statistically significance. The differentially expressed genes (DEGs) between different subtypes were screened using the "limma" package based on the selection criteria of log jFCj ≥ 1 and p < 0:05. Then, the DEGs were subjected to univariate Cox regression analysis. The DEGs meet the screening criteria p < 0:05 and were considered as significant DEGs for subsequent analysis. Gene Ontology (GO) functional enrichment analysis was performed using the "clusterProfiler" package to explore the potential molecular function of significant DEGs [24].

Construction and Validation of a Prognostic Risk Model.
The significant DEGs were included in Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis using the "glmnet" package in R software, and a 10-fold cross-validation/leave-one-out was performed to avoid model overfitting [25]. The significant DEGs with nonzero regression coefficients obtained by LASSO regression analysis were subjected to multivariate Cox regression analysis to further narrow down the genes and build a risk model. The risk score of each patient is calculated using the following formula: risk score = Σ ðexpression value of each gene × and its coefficientÞ. The HNSCC patients were clustered into high-and low-risk group based on the median cut-off of risk score. The Kaplan-Meier (KM) 2 Oxidative Medicine and Cellular Longevity curve was plotted to evaluate the prognosis of the risk model using the "survminer" R package. The receiver operating characteristic (ROC) curves at 1, 3, and 5 years were drawn to assess the prognostic predictive performance of the risk model using the "survival ROC" R package. Univariate and multivariate Cox regression analyses was used to identify whether the risk model is an independent prognostic factor for HNSCC. The TCGA-HNSC dataset was randomly divided into TCGA-training dataset (n = 250) and TCGA-testing dataset (n = 249) to confirm the performance of the risk model by internal validation. External validation was performed in GSE65858 dataset (n = 270). The prognostic predictive performance of the risk model was validated in internal and external validation using the same methods mentioned above.    60  256  112  30  None  1  0  0  Sex  Female  134  47  7  Male  368  223  61  None  0  0  0  OS_Event  Alive  283  176  51  Dead  218  94  17  None  1  0  0  PFS_Event  FALSE  137  TRUE  133  None  0  Grade  1  6 2  2  300  3  119  4  2  None  19  Stage  I  2 5  1 8  6  II  81  37  12  III  90  37  6  IV  306  178  44  None  0  0  0  Alcohol_history  Yes  333  239  59  No  158  31  9  None  11  0  0  Hpv16_status  Positive  31  60  11  Negative  72  209  57  None  399  1  0  Perineural_invasion  Yes  165  No  188  None  149  New tumor event after initiative treatment  Tumor free  275  With tumor  107  None  120  Smoking_history  Yes  222  53  No  48  15  None  0  0   4 Oxidative Medicine and Cellular Longevity Real time quantitative polymerase chain reaction (RT-qPCR) was performed based on SYBR Premix Exaq (Takara, Japan). GAPDH was used as an internal reference to calculate the relative expression levels of genes in healthy samples and HNSCC samples according to the 2-ΔΔCt method. Supplementary Table 1 presents the primer sequences of the genes. We then compared the differential expression level of genes between healthy samples and HNSCC samples. Finally, we validated the prognostic predictive performance of the risk model based on the 68 clinical specimens using the same methods mentioned above.

Construction of a Nomogram Model.
To improve the clinical applicability, we constructed a nomogram model combined risk score and clinical information to predict the survival of HNSCC patients at 1, 3, and 5 years using the "rms" R package [26]. Calibration curve was used to assess the differential predicted OS probability between the actual and the nomogram model. Decision curve analysis (DCA) curve and ROC curve were used to compare the differential performance of the nomogram to risk score and clinical information.

Exploration of the Clinicopathological Features and
Stemness Characteristics of the Prognostic Risk Model. The "compare" R package was used to compare the risk score in different cluster and clinicopathological features including age, sex, stage, grade, hpv16 status, and alcohol history. Gene mutation rate and tumor mutation burden (TMB) between high-and low-risk groups was compared by Wilcox test.
The "maftools" R package was used to visualize the differential gene mutation in high-and low-risk groups [21]. The correlation between TMB and risk score was identified by Spearman correlation analysis. The statistical significance was set at P < 0:05.
2.9. Tumor Immune Characteristics Analysis. The abundance of tumor immune cell infiltration in HNSCC samples was calculated by single sample gene set enrichment analysis (ssGSEA) algorithm [21]. The "estimate" R package was used to calculate the stromalscore, immunescore, and ESTI-MATEScore of the HNSCC samples. Wilcox test was used to compare the differential immune cell infiltration, immune checkpoint genes expression, stromalscore, immunescore, and ESTIMATEScore in different groups (high-vs low-risk group). Spearman correlation analysis was used to analyze the correlation between immune cell infiltration abundance and genes and risk score.
2.10. Statistical Analysis. One-way ANOVA and Kruskal-Wallis tests were used to compare differences between groups. Kaplan-Meier (K-M) curve was plotted for prognostic analysis in high-and low-risk groups. The "RCircos" R package was used to present the CNV of the DE FPRGs with prognostic value in chromosomes [27]. The "forestplot" R package was performed to calculate and visualize the hazard ratios (HR) of the DE FPRGs in TCGA-HNSC dataset [22]. All parametric analyses were based on two-tailed tests, the statistical significance of which was set at P < 0:05. All statistical analyses were performed using R 4.0.0.

Landscape of 21 FPRGs in TCGA-HNSC Dataset.
We obtained 91 DE FPRGs between normal and HNSCC samples through differential expression analysis (Figure 2(a), Supplementary table 2) and 48 prognostic related ferroptosis and pyroptosis by univariate Cox regression analysis (Supplementary table 3). Then, 21 integrated FPRGs were required and visualized by Venn diagram (Figure 2 3.2. Unsupervised Clustering Based on 21 FPRGs. The "Con-sensusClusterPlus" R package was performed to cluster the HNSCC patients in TCGA-HNSC dataset into two different subtypes based on the expression of the 21 FPRGs (cluster A and cluster B, Figure 3(a)). PCA revealed that we can completely distinguished cluster A and cluster B based on the expression level of the 21 FPRGs (Figure 3(b)). K-M analysis for the two different subtypes revealed that the patients in cluster B group have poor outcome than the patients in cluster A group (Figure 3(c)). In addition to the gene GZMB, the other 20 FPRGs were higher expressed in cluster B group than in cluster A group (Figures 3(d) and 3(e)). To compare the different biological behaviors between the two subtypes, GSVA analysis was performed. As shown in Figure 3(f), we found that cluster B mainly enriched in pathways associated with malignant progression of cancer such as _MAPK_SIGNALING_PATHWAY, P53_SIG-NALING_PATHWAY, and CHEMOKINE_SIGNALING_ PATHWAY, which verifies the poor prognosis of cluster B patients (Figure 3(f)). To further explore the biological function of each subtype, we identified 165 prognostic DEGs through differential expression analysis (Supplementary Table 4) and univariate Cox regression analysis (Supplementary Table 5). GO enrichment analysis was 5 Oxidative Medicine and Cellular Longevity     Oxidative Medicine and Cellular Longevity       14 Oxidative Medicine and Cellular Longevity  0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Time (years)   0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Time    16 Oxidative Medicine and Cellular Longevity The HNSCC patients in TCGA-HNSC dataset were randomly divided into TCGA-training dataset and TCGA-testing dataset to confirm the performance of the risk model by internal validation. In the TCGA-training dataset, the patients in low-risk group have longer OS than in highrisk group (Figure 6(a)). The AUC values at 1-, 3-, and 5year OS predicted by the risk model were 0.741, 0.768, and 0.811, respectively ( Figure 6(b)). PCA analysis indicated that the expression level of the 9 prognostic DEGs can separate the patients in high-risk group from the low-risk group completely ( Figure 6(c)). The risk score, survival state, and heat map of the 9 prognostic DEGs in the TCGA-training dataset are presented in Figure 6  21 Oxidative Medicine and Cellular Longevity multivariate Cox regression analysis indicated that the risk score was an independent prognostic predictor for OS ( Figure 6(e)). In the TCGA-testing dataset, the OS in the high-risk group was shorter than the low-risk group (Figure 6(f)). The AUC values at 1-, 3-, and 5-year OS predicted by the risk model were 0.562, 0.638, and 0.707, respectively ( Figure 6(g)). The high-risk group can be separated from the low-risk group based on the 9 prognostic DEGs (Figure 6(h)). Figure 6(i) showed the risk score, sur-vival state, and heat map of the 9 prognostic DEGs in the TCGA-testing dataset. The risk score was an independent prognostic predictor for OS, as revealed by univariate and multivariate Cox regression analysis (Figure 6(j)).
The HNSCC patients in GSE65858 dataset was used to confirm the performance of the risk model by external validation. The OS and relapse-free survival (RFS) in the lowrisk group was longer than in the high-risk group (Figures 7(a) and 7(f)). The AUC values at 1-, 3-, and 5-     29 Oxidative Medicine and Cellular Longevity HNSCC samples. We then built a risk model based on the 9 prognostic DEGs using the same methods mentioned in TCGA-HNSC dataset. We successfully verified the favorable prognostic predictive performance of the risk model according to the Figures 8(a)-8(e). We also compared the differential expression of the 9 prognostic DEGs between healthy samples and HNSCC samples. The results suggested that AC117422.1, AC117422.1, AC128687.2, AL161431.1, and FCRL1 were elevated and LRADT1, PDCL2, PLA2G3, and SPRR3 were declined in the HNSCC samples (Figure 8(f)). Finally, univariate and multivariate Cox analysis indicated that the risk score was the independent factor for the prognosis of the HNSCC patients (Table 2).

Construction of a Nomogram Model.
To improve the clinical applicability, a nomogram model combined risk score and clinical information was constructed (Figure 9(a)). Calibration curve at 1-, 3-, and 5-year indicated that the predicted OS probability of the nomogram model was close to the actual (Figure 9(b)). DCA curve suggested that the nomogram model has the highest net benefit compared with the individual features (Figure 9(c)). ROC

32
Oxidative Medicine and Cellular Longevity curve indicated that the nomogram model has the optimum sensitivity and specificity in prognostic prediction than the individual features (Figure 9(d)).

Exploration of the Clinicopathological Features and
Stemness Characteristics of the Prognostic Risk Model. The mutation frequency of the TP53 associated with adverse outcome of cancer was found higher in high-risk group (67%, Figure 10(a)), compared with the low-risk group (58%, Figure 10(b)). To investigate the relationship between risk score and clinicopathological features, the "compare" R package was performed. We observed that the patients with age ≤60 had higher risk score than the patients with age>60, whereas the risk score has no statistic difference in other stratified clinicopathological features (Figures 10(c)-10(h)). The patients in cluster A group corresponds to higher risk score (Figure 10(i)) and the relationship of cluster, risk, and fustat was showed in Sankey diagram (Figure 10(l)). Increasing evidence revealed that the patients with higher TMB can more benefit from immunotherapy [28]. Figures 10(j) and 10(k) showed that the risk score has negative correlation with TMB, suggesting that the patients in low-risk group were more sensitive to immunotherapy.
3.6. Tumor Immune Characteristics Analysis. The abundance of immune cells was calculated using the ssGSEA algorithm, and the different levels of immune cell infiltration between high-and low-risk group was compared, finding that the HNSCC samples in low-risk group has increased immune response (Figure 11(a)). We also investigated the differential expression levels of immune checkpoints, finding that there were 19 immune checkpoints overexpressed in the low-risk group than that in the high-risk group (Figure 11(b)). In addi-tion, the HNSCC samples in low-risk group were related to higher immune score and ESTIMATEScore (Figure 11(c)). The network presented the interactions, regulator connection, and prognostic value of the 23 types of immune cells (Figure 11(d)). We found that there were strong positive correlation and mutual regulation between the 22 types of immune cells (except CD56dim natural killer cells, Supplementary table 7). Combined with the network graph and K-M curve (Supplement Figure 1), 17 types of immune cells were associated with the prognosis of the HNSCC patients. Among them, 16 types of immune cells were protective factors, and neutrophilia was a risk factor. Spearman correlation analysis was performed to evaluate the correlation between the 9 prognostic DEGs in risk model and the 23 types of immune cells. We observed that FCRL1 has the strongest positive correlation with the 23 types of immune cells, while the AC128687.2 has the strongest negative correlation with the 23 types of immune cells (Figure 12(a)). Figure 12(b) revealed that the infiltration abundance of 17 immune cells was reduced as the risk score.

Discussion
Head and neck tumors commonly occur in the oral cavity, nasopharynx, oropharynx, hypopharynx, and larynx [29].
As the most common pathologic type, HNSCC ranks 6th in the incidence rate of malignant tumors worldwide, and more than 800,000 new cases are diagnosed every year [30]. At present, surgery is the main treatment, supplemented by radiotherapy and chemotherapy, but the 5-year survival rate is still not ideal [31]. Especially, local recurrence and distal organ metastasis often occur in advanced HNSCC after treatment, with higher mortality [32]. Targeted

35
Oxidative Medicine and Cellular Longevity and immunotherapy are promising for patients with advanced HNSCC. Currently, cetuximab is the only molecular-targeted drug approved for clinical treatment of HNSCC. Combined with platinum-based chemotherapy, cetuximab is the standard treatment for patients with advanced HNSCC [33,34]. Nivolumab and pembrolizumab are two immunocheckpoint inhibitors permitted for immu-notherapy in HNSCC patients [33,35]. However, drug resistance and immune escape are major problems faced by targeted and immunotherapy. Induction of apoptosis to halt tumor growth is the aim of many HNSCC treatment strategies [13]. Ferroptosis and pyroptosis are two new PCD modes discovered in recent years and play important roles in the malignant processes and immune microenvironment

38
Oxidative Medicine and Cellular Longevity of HNSCC. Therefore, target ferroptosis and pyroptosisrelated genes may be improved the prognosis of HNSCC.
In this study, we successfully constructed and verified a risk model based on the FPRGs. Firstly, we obtained 21 significant FPRGs differential expression analysis and univariate Cox regression analysis. The patients in TCGA-HNSC dataset were divided into two different subtypes based on the expression of the 21 FPRGs using the "ConsensusClusterPlus" R package. Further, we acquired the 165 prognostic DEGs between the two subtypes to explore the molecular differences of the two subtypes. GO enrichment analysis indicated that the 165 prognostic DEGs mainly enriched in GO:0005198~structural molecule activity, GO:0005882~intermediate filament, GO:0005615~extracellular space, etc. These extracellular components were an important part of the immune microenvironment [36]. Ferroptosis and pyroptosis have been widely reported to have extremely complicated crosstalk with tumor immune microenvironment [5]. The 165 DEGs with prognostic significance was subjected into the LASSO-multivariate Cox regression analysis to a built a risk model with 9 prognostic    CD244  CD27  CD28  CD40LG  CD44  CD48  CTLA4  ICOS  LAG3  LAIR1  LGAS9  PDCD1  TIGIT  TMIGD2  TNFRSF14  TNFRSF4  TNFSF4

40
Oxidative Medicine and Cellular Longevity DEGs. K-M curve revealed that the patients with higher risk score have poor outcome compared the patients with lower risk score. ROC curve indicated that the sensitivity and specificity of the risk score for prognostic prediction in HNSCC patients was favorable. More importantly, the risk model was successfully verified with a stable prognostic value through internal validation, external validation, and clinical sample validation. In addition, the risk score was found to be independent of other clinical information in predicting the prognosis of HNSCC. Regardless of the clinical characteristics, stratified    43 Oxidative Medicine and Cellular Longevity prognostic analysis showed that the HNSCC patients in the high-risk group continued to have poor outcomes except for the HPV+ group. This result may be due to the small sample size of HPV+ group (n = 31). Finally, a nomogram model combined risk score and clinical information was constructed to improve the clinical applicability. Reviewing previous studies, some of these 9 prognostic DEGs have been found to be involved in the occurrence and progression of solid tumors. For example, Lu Yu et al. [37] reported that the expression level of SPRR3 was reduced as the malignant progression of oral squamous cell carcinoma (OSCC), which was consistent with our analysis results. AL161431.1 was found upregulated in endometrial carcinoma, lung cancer, and pancreatic cancer and associated with the immune microenvironment, proliferation, migration, epithelial-mesenchymal transformation, and poor prognosis [38,39]. Randall S. Davis [40] identified that FCRL1 overexpressed in breast, melanoma, and lung cancer may be a potential biomarker and therapeutic target. LRATD1 also named FAM84A has been revealed to be related to the occurrence and development of papillary thyroid cancer, liver tumor, and colon cancer [41][42][43]. PLA2G3 was upregulated in ovarian cancer, melanoma, and colorectal cancer and improved the poor prognosis and malignant progression of cancer [44][45][46].
TP53 is one of the most frequently altered genes in human cancers, which is present in about 50% of invasive tumors [47]. Genomic data showed that TP53 was the most common mutant gene in HNSCC and associated with shorter survival outcome of HNSCC patients [48]. Our research found that the patients in high-risk group with poor survival outcome have higher TP53 mutation frequency (67%) than the patients in low-risk group (58%), which was consistent with the results of previous studies. TMB refers to the number of nonsynonymous mutations in somatic cells per mega base pair (Mb) in a specific genomic region, which can indirectly reflect the ability of tumor to produce neoantigens [49]. Patients with higher TMB are more likely to benefit from immunotherapy [49]. In this study, we found that the risk score has negative correlation with TMB, suggesting that the patients in low-risk group were more sensitive to immunotherapy.
To reveal the mechanism of risk model in tumor immune microenvironment, we firstly calculated the abundance of immune cells using the ssGSEA algorithm and compared the differential immune cell infiltration among the high-and low-risk group. The results identified that Activated.B.cellna, Activated.CD8.T.cellna, Eosinophilna, Immature.B.cellna, MDSCna, Macrophagena, Mast.cellna, Monocytena, Natural.killer.cellna, T.follicular.helper.cellna, Type.1.T.helper.cellna, and Type.17.T.helper.cellna had decreased infiltration as the risk score increased. K-M curve revealed that the above 12 type of immune cells except T.follicular.helper.cellna are the protective prognostic factor for HNSCC. B cells are the main effector cells of humoral immunity, which can directly kill tumor cells and inhibit tumor development by secreting immunoglobulin [50]. Xin Feng et al. showed that the B cells act a favorable role in the prognosis of HNSCC. The higher infiltration of B cell and their subtypes may improve the prognosis of HPV+ HNSCC patients [51]. Activated CD8 T cells as the most important antitumor effector cells can recognize tumor associated antigens by expressing T cell receptors and kill tumor cells [52]. Many studies have shown that the HNSCC patients can benefit from the increased infiltration of the activated CD8 T cells [53]. Eosinophils can kill tumors directly or indirectly by releasing cytotoxic proteins or chemoattractants, which may extend the prognosis of HNSCC patients [54]. As the first line of defense against tumor, natural killer cells have been reported to play an important role in antitumor immunity of HNSCC [54]. T helper cells, as the most important helper cells in tumor immunity, can promote the recruitment of natural killer cells to the tumor and activate death receptors on the surface of tumor cells and the CD8 T cells by releasing cytokines [21]. Contrary to our analysis, monocytes, macrophages, mast cells, and myeloid derived suppressor cells (MDSCs) were considered to relate to the malignant progression and poor prognosis of HNSCC [55][56][57][58].