A Ferroptosis-Related Prognostic Signature Based on Antitumor Immunity and Tumor Protein p53 Mutation Exploration for Guiding Treatment in Patients With Head and Neck Squamous Cell Carcinoma

Background: Due to the lack of accurate guidance of biomarkers, the treatment of head and neck squamous cell carcinoma (HNSCC) has not been ideal. Ferroptosis plays an important role in tumor suppression and treatment of patients. However, tumor protein p53 (TP53) mutation may promote tumor progression through ferroptosis. Therefore, it is particularly important to mine prognostic-related differentially expressed ferroptosis-related genes (PR-DE-FRGs) in HNSCC to construct a prognostic model for accurately guiding clinical treatment. Methods: First, the HNSCC data obtained from The Cancer Genome Atlas (TCGA) was used to identify PR-DE-FRGs for screening candidate genes to construct a prognostic model. We not only used a variety of methods to verify the accuracy of the model for predicting prognosis but also explored the role of ferroptosis in the development of HNSCC from the perspective of the immune microenvironment and mutation. Finally, we explored the correlation between the prognostic model and clinical treatment and drew a high-precision nomogram to predict the prognosis. Results: Seventeen of the 29 PR-DE-FRGs were selected to construct a prognostic model with good predictive performance. Patients in the low-risk group were found to have a greater number of CD8 + T cells, follicular helper T cells, regulatory T cells, mast cells, T-cell costimulations, and type II interferon responses. A higher tumor mutation burden (TMB) was observed in the low-risk group and was associated with a better prognosis. A higher risk score was found in the TP53 mutation group and was associated with a worse prognosis. The risk score is closely related to the expression of immune checkpoint inhibitors (ICIs)-related genes such as PD-L1 and the IC50 of six chemotherapeutic drugs. The nomogram we constructed performs well in predicting prognosis. Conclusion: Ferroptosis may participate in the progression of HNSCC through the immune microenvironment and TP53 mutation. The model we built can be used as an effective predictor of immunotherapy and chemotherapy effects and prognosis of HNSCC patients.

A Ferroptosis-Related Prognostic Signature Based on Antitumor Immunity and Tumor Protein p53 Mutation Exploration for Guiding Treatment in Patients With Head and Neck Squamous Cell Carcinoma Xin Fan 1,2 , YangShaobo Ou 1,2 , Huijie Liu 1,2 , Liangzhen Zhan 3 , Xingrong Zhu 1,2 , Mingyang Cheng 1,2 , Qun Li 1,2 , Dongmei Yin 1,2 and Lan Liao 1,2 * Background: Due to the lack of accurate guidance of biomarkers, the treatment of head and neck squamous cell carcinoma (HNSCC) has not been ideal. Ferroptosis plays an important role in tumor suppression and treatment of patients. However, tumor protein p53 (TP53) mutation may promote tumor progression through ferroptosis. Therefore, it is particularly important to mine prognostic-related differentially expressed ferroptosis-related genes (PR-DE-FRGs) in HNSCC to construct a prognostic model for accurately guiding clinical treatment.
Methods: First, the HNSCC data obtained from The Cancer Genome Atlas (TCGA) was used to identify PR-DE-FRGs for screening candidate genes to construct a prognostic model. We not only used a variety of methods to verify the accuracy of the model for predicting prognosis but also explored the role of ferroptosis in the development of HNSCC from the perspective of the immune microenvironment and mutation. Finally, we explored the correlation between the prognostic model and clinical treatment and drew a high-precision nomogram to predict the prognosis.
Results: Seventeen of the 29 PR-DE-FRGs were selected to construct a prognostic model with good predictive performance. Patients in the low-risk group were found to have a greater number of CD8 + T cells, follicular helper T cells, regulatory T cells, mast cells, T-cell costimulations, and type II interferon responses. A higher tumor mutation burden (TMB) was observed in the low-risk group and was associated with a better prognosis. A higher risk score was found in the TP53 mutation group and was associated with a worse prognosis. The risk score is closely related to the expression of immune checkpoint inhibitors (ICIs)-related genes such as PD-L1 and the IC50 of six chemotherapeutic drugs. The nomogram we constructed performs well in predicting prognosis.

INTRODUCTION
Head and neck squamous cell carcinomas (HNSCCs) are neoplasms affecting different tissues and organs in the head and neck region including the tongue, mouth, nasopharynx, larynx, and throat (Fang et al., 2018). As the sixth most common cancer in the world, head and HNSCC affects more than 700,000 new patients each year, with a mortality rate of approximately 40-50% (Bray et al., 2018;Cohen et al., 2019). Due to the lack of effective early monitoring and screening factors, HNSCC is difficult to be detected in the early stages (Cui et al., 2021). At present, surgery, radiotherapy and chemotherapy based on the location and clinical stage of HNSCC are still identified as commonly used clinical methods for the treatment of HNSCC (Yi L. et al., 2020). However, as one of the core treatment options, recurrence induced during radiotherapy has become the main reason for the failure of HNSCC treatment . The treatment that can be used is very limited in view of the poor prognosis of advanced HNSCC. In recent years, the prognosis of patients with advanced malignant tumors has been significantly improved with the clinical application of more and more cancer-corresponding immune checkpoint inhibitors (ICIs). After the anti-PD-1 monoclonal antibodies pembrolizumab and nivolumab were approved for the treatment of recurrent or metastatic HNSCC in 2016, first, ICIs have shown increasing potential in the treatment of HNSCC (Quan et al., 2019). Unfortunately, during HNSCC treatment, ICIs may cause severe morbidity or, in rare cases, mortality, and the objective response rate is only 20% (Wang et al., 2019a;Tada et al., 2020). Therefore, it is urgent to find highly sensitive markers that can predict the benefits of various clinical treatments to significantly reduce the side effects of treatment.
Ferroptosis is a kind of regulatory necrosis driven by irondependent phospholipid peroxidation, which is regulated by cell metabolism, redox homeostasis, and various signaling pathways related to cancer (Yi J. et al., 2020). As a non-apoptotic form, the role of ferroptosis in the physiological processes of many diseases, especially its important role in tumor suppression, has been confirmed by many studies (Jiang et al., 2015;Zhang et al., 2018). Recent results indicated that ferroptosis mediates the tumor suppressor activity of interferon gamma secreted by CD8 + T cells in response to immune checkpoint blockade, indicating that the immune system may partly prevent tumorigenesis through ferroptosis (Wang et al., 2019c). Luo et al. (2021) also concluded through meta-analysis that ferroptosis could inhibit tumor growth and contribute to chemotherapy sensitivity. In addition, the activation of ferroptosis has also been proven to help cancer treatments, such as immune checkpoint blockade (Wang et al., 2019c) and radiotherapy (Lei et al., 2020). Ferroptosis induction to achieve cancer cell death may provide a new perspective for the treatment of metastatic cancer and malignant tumors resistant to traditional therapies (Liang et al., 2019;Xia et al., 2020).
As a reflection of the number of mutations in tumor cells, tumor mutation burden (TMB) has been observed to be closely related to the efficacy of ICIS and prognosis in many studies (Devarakonda et al., 2018;Cao et al., 2019;Samstein et al., 2019;Wang et al., 2019c;Wu et al., 2019). As a tumor suppressor gene, missense mutations in the TP53 gene are common in human cancers (Mantovani et al., 2019). In many sporadic cancers, TP53 mutations are associated with poor prognosis (Kandoth et al., 2013). Mutated p53 not only no longer plays the role of suppressing cancer but also assists in the development of cancer by depriving the cell of the tumor suppressor response (Mantovani et al., 2019). In addition to the loss of tumor suppressor function, TP53 mutations also inhibit ferroptosis by altering cellular iron acquisition and metabolism, thereby promoting cancer progression . Therefore, it is very likely to find some potential roles of ferroptosis in the development of HNSCC from the perspective of TP53 mutation.
In view of the unsatisfactory treatment methods and lack of ideal models to guide treatment, it is particularly important to explore the ferroptosis-related genes that play a role in HNSCC and construct a new effective prognosis model to predict the prognosis of HNSCC and guide clinical treatment precisely. Moreover, the role of ferroptosis in HNSCC is unclear, and it is urgent to explore the mechanism of ferroptosis from the perspective of immunity and mutation. For this purpose, we used the HNSCC cohort data from The Cancer Genome Atlas (TCGA) to screen the prognostic-related differentially expressed ferroptosis-related genes (PR-DE-FRGs) to construct a prognostic model. We not only used a variety of methods to verify the performance in predicting prognosis of the model but also explored the potential mechanism of the ferroptosisrelated model in the development and prognosis of HNSCC from the perspective of immune microenvironment and mutation. Taking into account the close correlation between ferroptosis and clinical treatment, we have further explored the clinical application value of prognostic models, including the guiding value of ICIS treatment, chemotherapy, and the combined model constructed by the prognostic model to predict the prognosis.
In general, our research will provide a model of ferroptosisrelated genes with good predictive performance. In addition, this model performs well in clinical applications. It can be used to guide immunotherapy and chemotherapy and predict the prognosis of HNSCC patients. In addition, our research will also explore the possible role of the immune microenvironment and TP53 mutation in the development of HNSCC, providing a new perspective for subsequent research.

Data Acquisition and Identification of Differentially Expressed Ferroptosis-Related Genes
The RNA-sequencing data in level 3 and corresponding clinical data of 504 HNSCC samples and 44 adjacent normal tissues were downloaded from TCGA database 1 for the analysis on February 3, 2021. In addition, somatic gene mutation and the copy number variations (CNVs) data on TCGA also have been downloaded. The lists of FRGs used in our analysis were derived from the FerrDb Database, 2 which is the first database known for the manual collection and management of ferroptosis-related markers and regulatory factors and ferroptosis-related diseases in the world. A total of 258 FRGs from 784 articles on ferroptosis in PubMed database were included in FerrDb. The annotation of these genes revealed 108 driver genes, 69 suppressor genes, and 111 gene markers. In view of the fact that these data were both publicly available in TCGA, the approval of the local ethics committee was exempted for this research. The gene expression quantification of the same 246 FRGs between the lists of FRGs and RNA-sequencing data of HNSCC samples from TCGA was extracted. R package limma was used for differential expression analysis between 504 HNSCC samples and 44 adjacent normal tissues to identify the DE-FRGs. Since there were fewer FRGs, our analysis set the filter conditions to | log2 fold change| (| log2FC| > 0) along with false discovery rate <0.01.

Construction of Regulatory Networks and Prognostic Model Based on Prognostic Differentially Expressed Ferroptosis-Related Genes
We integrated data from TCGA to obtain the same samples that both held complete overall survival (OS) and messenger RNA (mRNA) expression data. Univariate Cox analysis with the cutoff value of p < 0.05 was performed to select FRGs with prognostic values. After the minimum required interaction score was set at medium confidence (0.400), protein-protein interaction (PPI) network consisting of 29 prognostic DE-FRGs was generated in the STRING database, which was widely known for searching known protein interaction relationships online. In addition, R package reshape2, psych, RColorBrewer, and igraph were used for correlation analysis of 29 genes and visualization of related networks. A total of 499 samples with complete OS data and mRNA expression value of prognostic DE-FRGs were randomly matched to the training set (n = 300) and the test set (n = 199) at a ratio of 6:4. The expression values of 29 prognostic DE-FRGs from the training group were used in the least absolute shrinkage and selection operator (LASSO) regression analysis, which could screen out highly related genes and minimize the risk of overfitting for screening signatures to accurately forecast the clinical efficacy in HNSCC cases. Finally, a prognostic model based on 17 DE-FRGs was constructed by selecting the optimal penalty parameter, which was determined by the minimum 10fold cross-validation. Next, we used the coefficients obtained by the LASSO regression algorithm in the following risk score equation: risk score = sum of corresponding coefficients × its gene's expression value (Gao et al., 2010). In addition, Kaplan-Meier survival curve was used to analyze the relationship between expression of each gene in the model and OS.

Validation of the Prognostic Model
The training set, test set, and whole set were used to evaluate the predictive power of the prognostic model simultaneously. However, three sets obtained from the TCGA data were from a single cohort. The risk score of each case that would be divided into high-and low-risk groups according to the median value of their risk scores was calculated separately in the training set, test set, and whole set. As soon as risk score was obtained, we used the R tool to visualize the specific risk score value and survival status of each sample in each set in the model. R package survival and survminer were used to perform Kaplan-Meier to verify the performance of the model to distinguish the survival of HNSCC patients. The multifactor receiver operating characteristic (ROC) curves were used not only to verify the accuracy of the model but also to optimality of the model compared to other factors in predicting survival. We also performed principal component analysis (PCA) with prcomp function of the R package stats based on the expression of genes in each set. In addition, univariate and multivariate cox regression analysis were performed on the variables (age, gender, grade, stage, T, N, and risk score) from three sets to determine whether the risk score could be used as an independent prognostic predictor of OS.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis
The R package limma was utilized to determine the differential genes (DEGs) between the high-and low-risk group among the three sets based on the filter condition (| log2FC | ≥ 1, FDR < 0.05). We used R package clusterProfiler to run Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) to elucidate which functions and pathways of these DEGs are enriched.

The Relationship Between Prognostic Model Immune Cells/Immune Functions
To further reveal the correlation between the prognostic model and immune tumor immune microenvironment, the singlesample gene set enrichment analysis based on R packages GSEABase and gsva was used to quantify the scores of immune cells and functions. The difference analysis of the scores of immune cells and functions between different risk groups based on the prognostic model were performed to draw a box plot.
Since ferroptosis is closely related to cancer immunity, we also downloaded the data of the content of T cells CD8 + and regulatory T cells (Tregs) in each HNSCC sample calculated by six advanced algorithms from the TIMER 2.0 database, 3 including TIMER, MCPCOUNTER, CIBERSORT, XCELL, QUANTISEQ, and EPIC. In further analysis, we separately analyzed the correlation between the content of these cells and the PR-DE-FRGs used in constructing the model.

Mutation and Copy Number Variation Analysis
We downloaded the somatic gene mutation data and corresponding clinical data of HNSCC samples in the TCGA dataset. These mutation data were used to evaluate the relationship between the mutation and the prognostic model of HNSCC patients who owned available mutation data. After using VarScan to detect the mutation annotation format (MAF) files of somatic mutations in HNSCC samples, the R package GenVisR was used to visualize the 30 most frequently mutated genes in the high-and low-risk groups, respectively, (Koboldt et al., 2009). The R package maftools was used to calculate the TMB who was defined as the number of somatic, coding, indels mutations, and base substitutions in a million bases of the genome. To explore the influences of TMB on survival, Kaplan-Meier (KM) analysis was used to compare the differences in survival between highand low-TMB groups. According to the TP53 mutation status, the TCGA samples were assigned to the wild group and the mutant group. To explore the association between ferroptosis and TP53 mutation, we compared the difference in risk scores between the TP53 mutation group and TP53 wild group. Not only that, KM analysis was used to compare the difference in OS. We extracted the CNV data of 16 genes used to build the model in 526 HNSCC samples. After statistics on the CNV frequency of these genes, the corresponding results were visualized. The locations of CNV changes in these 16 genes on the chromosome were also visualized in the circle graph.

Correlation Between Prognostic Model and Clinical Treatment
The immunophenoscore (IPS) data of HNSCC patients were obtained from The Cancer Immunome Atlas (TCIA) database. 4 The patient's IPS was obtained by evaluating the gene expression of the four cell types (effector cells, immunosuppressive cells, MHC molecules, and immunomodulators) that determine immunogenicity . As the IPS score increases, the immunogenicity increases . To analyze the potential correlation between prognostic models and immunotherapy, Spearman correlation analysis was run to evaluate the correlation between risk score and four types of IPS and the expression of eight ICIs-related genes (CD274, CD40,CTLA4, PDCD1, CD96, TIGIT, LAG3, and HAVCR2), 3 timer.comp-genomics.org 4 tcia.at/home respectively. The results of the correlation were visualized in the form of a lollipop graph through the R package ggplot2. To predict the sensitivity of HNSCC patients to chemotherapy drugs recommended by National Comprehensive Cancer Network (NCCN) guidelines in the clinic, the R package pRRophetic was used to predict the half maximal inhibitory concentration (IC50) of different risk sample to chemotherapy drugs. The IC50 difference between high-and low-risk populations was compared. The cell line expression data in Cancer Drug Sensitivity Genomics (GDSC) database and RNA sequencing transcriptome data in TCGA database were used to construct the ridge regression model to predict the IC50 of the drug in this R package (Geeleher et al., 2014).

Construction and Evaluation of Nomogram
To create a more clinically applicable quantitative tool for predicting 1-, 2-, and 3-year OS of HNSCC patients, the factors that may affect survival containing risk scores and clinical risk factors were used to construct a multivariate cox regression model by the R package rms to draw nomogram. The calibration curve was then used to evaluate the accuracy of survival prediction in the nomogram. In addition, we not only drew the ROC curves in 1, 2, and 3 years of nomogram but also compared ROC curves of nomogram with risk score and other clinical factors to verify the optimality of nomograph in predicting prognosis.

Immunohistochemical Staining
In order to further verify the differential expression of the 17 PR-DE-FRGs used to construct the model between HNSCC tissues and normal tissues of the head and neck, we used the result of immunohistochemistry (IHC) staining from the human protein atlas database. 5 Finally, the immunohistochemical images of the protein expression of 15 PR-DE-FRGs in HNSCC tissues and normal tissues of the head and neck were obtained.

Statistical Analysis
According to the distribution characteristics, Student's t-test or Mann-Whitney test was used to compare the difference between continuous variables, while the chi-square test or Fisher exact test was used to compare the difference between categorical variables. Univariate Cox regression analysis was used to confirm DE-FRGs related to prognosis, which would be screened to construct the optimal prognostic model by LASSO regression. The nomogram was constructed using multivariate Cox regression. Kaplan-Meier curve with log-rank test was utilized to compare survival between different groups. The prognostic ability of each factor for prognosis was evaluated by ROC curve. The independent predictive value of risk scores of prognostic model for survival was evaluated with univariate and multivariate Cox regression analyses. All analyses in our study were performed in R programming language (version 4.0.3) and SPSS Statistics software 22.

Identification of Differentially Expressed Ferroptosis-Related Genes in Head and Neck Squamous Cell Carcinoma
We showed the workflow of this study in Figure 1. In the first, gene transfer format files from ensembl were used to annotate the data of 504 HNSCC and 44 adjacent normal tissues. The gene expression quantification of 246 FRGs were extracted to identify 157 DE-FRGs (Figure 2A). Among 157 FRGs, 36 genes were identified as downregulated and 121 genes were identified as upregulated ( Figure 2B).

Validation of Prognostic Model
The distribution of risk scores, OS, and OS status of 17 signatures in training, test, and whole sets, which all indicate that as the risk score increases, there is an increasing trend of deaths, is shown in Figures 5A-C. After Kaplan-Meier analysis, we all found that patients in the high-risk group dramatically had lower survival probability in three sets (Figures 5D-F). Figures 5G-I showed the PCA results of the three sets respectively. The univariate COX regression analysis showed that the risk score of the training, test, and whole sets were all significantly correlated with OS (HR = 1.323, 95% CI = 1.193-1.466, p < 0.001; HR = 1.217, 95% CI = 1.130-1.311, p < 0.001; HR = 1.253, 95% CI = 1.183-1.327, p < 0.001) (Figures 5J-L). Risk score was still identified as an independent predictor of OS (HR = 1.287, 95% CI = 1.157-1.431, p < 0.001; HR = 1.196, 95% CI = 1.105-1.295, p < 0.001; HR = 1.242, 95% CI = 1.169-1.320, p < 0.001) in each set by the multivariate Cox regression analysis after adjusting for other clinical confounding factors (Figures 5M-O). To verify the optimality of the model, we drew ROC curve based on the training, test, and whole sets. Among the 3 years, the results in each set showed that all the area under the curve (AUC) values are above 0.7 (Figures 6A-I), which means that this model has better predictive value. It was worth mentioning that consistent results were obtained in the three sets: the risk scores have the largest AUCs of 1, 2, and 3 years in all factors, which indicates that the risk model has the best predictive effect (training set,

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis
The cellular components (CCs), molecular functions (MFs), and biological processes (BPs) enriched in training, test, and whole sets are shown in Figures 7A-C, respectively. Surprisingly, BPs, MFs, and CCs, which were enriched in all three sets, were almost all immune related, such as immune response-activating cell surface receptor signaling pathway, humoral immune response, adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, lymphocyte-mediated immunity, immunoglobulinmediated immune response, B-cell-mediated immunity, humoral immune response mediated by circulating immunoglobulin, complement activation classical pathway, complement activation, antigen binding, immunoglobulin receptor binding, external side of plasma membrane, immunoglobulin complex, immunoglobulin complex circulating, and blood microparticle. After KEGG pathway analysis, it was found that the cytokinecytokine receptor interaction pathway was enriched in all sets. Figures 7D-F, respectively, showed all results of KEGG pathway for training, test, and whole sets. The blue rectangles in Figures 7D,F also marked the 10 same pathways enriched in training and whole sets.

The Immune Infiltrating Cells/Immune Functions Content Associated With Prognostic Model
After difference analysis of the scores of 16 immune cells between different risk groups, we found that the score of B cells and mast cells were all significantly higher in the lowrisk group in training, test, and whole sets (adjusted p < 0.05, Figures 8A-C). It was worth mentioning that the scores of B cells statistically differed between the high-and low-risk groups (adjusted p < 0.001, Figures 8A-C), which was closely related with B-cell-mediated immunity enriched significantly in the GO analysis.
In addition, the scores of CD8 + T cells, immature dendritic cells (iDCs), T-helper cells, follicular helper T cells (Tfh), tumorinfiltrating lymphocyte (TIL), and regulatory T cells (Tregs) in the low-risk group all were significantly higher than those in the high-risk group in training and whole set (Figures 8A,C). In terms of 13 immune functions, in both the training and whole sets, it was worth mentioning that the scores of T-cell costimulation and type II IFN response were both observed to be significantly higher in the low-risk groups (Figures 8D,F). Figure  8E shows the difference analysis results of immune functions based on the test set data.
Finally, in addition to ZFP69B and LINC00336, the correlation results of the remaining 15 PR-DE-FRGs with T cells CD8 + /Tregs content were obtained and shown in Figure 8G. Through a comprehensive analysis of the results obtained by multiple calculation methods, it is shown that the expression of ASNS, BNIP3, FTH1, CISD2, SLC7A5, and ATG5 is negatively correlated with the content of CD8 + T cells, whereas the expression of MAP3K5, FBXW7, and SOCS1 is positively correlated with the content of CD8 + T cells. In addition, from the correlation graph, the content of Tregs was not only found to be significantly and negatively correlated with the expression of AURKA, BNIP3, and SLC7A5 but also was found to be significantly and positively correlated with the expression of MAP3K5, DRD4, SLC2A3, BAP1, FBXW7, MAP1LC3A, and SOCS1. Figures 9A,B, respectively, show the mutations of the top 30 most frequent genes of 240 samples in high-risk group and 226 samples in the low-risk group. In the box plot (Figure 9C), the TMB of patients in the low-risk group was shown to be significantly lower. It was observed that OS of samples with low TMB were better in further survival analysis ( Figure 9D). Patients in the TP53 mutation group were found to have a higher risk score ( Figure 9E). In addition, patients in the TP53 mutation group were observed to have a worse prognosis (Figure 9F). From Figure 9G, we observed that CNVs were common in these 16 genes. Except for BNIP3, SLC7A5, ATG5, FBXW7, AURKA, and BAP1, other genes have a higher frequency of acquiring  Figure 9G). The corresponding locations of these 16 genes on the chromosome and the synthesis of CNV are shown in Figure 9H.

Clinical Treatment Application of Prognostic Model
At present, ICIs such as PD-1 inhibitor pembrolizumab was being used clinically as the first-line treatment of relapsed HNSCC. Unfortunately, none of the four combinations of IPS showed a significant correlation with the risk score ( Figure 10A). However, we found that the expression of CD274 (PD-L1) was significantly positively correlated with the risk score, while that of CD40, PDCD1 (PD-1), CTLA4, CD96, and TIGIT were significantly negatively correlated with risk score (all p < 0.05; Figure 10B). This result suggests that our model may be used to predict the expression of PD-L1 in patients to predict the effect of immunotherapy. We predicted IC50 of each sample of six chemotherapy drugs (cisplatin, paclitaxel, docetaxel, gemcitabine, methotrexate, and BIBW2992) recommended by NCCN guidelines for HNSCC. It was found that except for methotrexate, the IC50 of the other five chemotherapy drugs in the low-risk group was higher (all p < 0.05), which implied that high-risk patients were more sensitive to the other five drugs, while low-risk patients were more sensitive to methotrexate ( Figure 10C).

Construction and Evaluation of Nomogram
Age, gender, T stage, N stage, and risk score were observed to be significantly related to OS by univariate cox regression analysis based on the training set (all p < 0.05, Figure 5J). From the nomogram constructed by these factors, we observed that risk score was the most important factor affecting patients' survival, followed by N stage, age, gender, and T stage ( Figure 11A). Through the calibration curve, we confirmed that the nomogram all had a good consistency between the predicted and the actual 1-, 2-, and 3-OS in training, test, and whole sets (training set, Figures 11B-D; test set, Figures 11E-G

Immunohistochemistry Staining
We obtained immunohistochemical staining images of the protein expression of 15 PR-DE-FRGs (ASNS, AURKA, FTH1, SLC2A3, SLC7A5, CISD2, PRDX6, ATG5, BAP1, MAP1LC3A, SOCS1, BNIP3, MAP3K5, ZFP69B, and FBXW7) in HNSCC and normal tissues of the head and neck (Figure 13). By comparing the images of immunohistochemical staining in HNSCC tissue and normal tissues of the head and neck, we found that 11 DE-FRGs (ASNS, AURKA, FTH1, SLC2A3, SLC7A5, CISD2, PRDX6, ATG5, BAP1, MAP1LC3A, and SOCS1) showed differences in proteins expression (Figures 13A-K). The protein expression of ASNS, AURKA, FTH1, SLC2A3, SLC7A5, CISD2, ATG5, BAP1, and SOCS1 in the HNSCC tissue is higher than that in normal tissues of the head and neck (Figures 13A-F,H,I,K). On the contrary, the protein expression of PRDX6 and MAP1LC3A was observed to be lower in the HNSCC tissue (Figures 13G,J). All these support the results of the differential expression of these genes in our analysis.

DISCUSSION
In the research, we obtained 29 PR-DE-FRGs through difference and univariate cox regression analysis and screened 17 PR-DE-FRGs through LASSO regression to construct a prognostic model. The prognostic model showed good performance in all data randomly divided into three sets. Through the correlation analysis between the immune infiltrating cells and immune functions and the model, many meaningful results have been observed. Also, from the perspective of mutation, a series of associations between mutation and ferroptosis, even the development and prognosis of HNSCC have been initially established. In terms of clinical application, the prognostic model shows good distinguishing performance in the expression of ICIS gene and the sensitivity of commonly used chemotherapy drugs. Not only that, the combined model constructed by the prognostic model is also excellent in predicting the prognosis of HNSCC patients. As we all know, different ferroptosis-related genes play different roles in ferroptosis. Among the 17 genes used to construct the prognostic model, 4 were identified as genes that inhibit ferroptosis (FTH1, CISD2, LINC00336, and PRDX6), 5 were identified as genes that drive ferroptosis (ATG5, BAP1, FBXW7, MAP1LC3A, and SOCS1), and 8 (ASNS, AURKA, BNIP3, DRD4 MAP3K5, SLC2A3, SLC7A5, and ZFP69B) were identified as genes that suggest ferroptosis in previous research.
Past studies have found that ferroptosis can inhibit tumor growth and contribute to chemotherapy, immune checkpoint blockade, and radiotherapy, which have a positive effect on the prognosis of patients (Wang et al., 2019c;Lei et al., 2020;Luo et al., 2021). The role of these ferroptosis-related genes in the prognosis of HNSCC samples is also supported in our study. First, in our analysis, FTH1, CISD2, and PRDX6 were identified as risk factors, while BAP1, FBXW7, MAP1LC3A, and SOCS1 were identified as protective factors. In further survival analysis, it was also found that patients with high expression of BAP1, FBXW7, MAP1LC3A, and SOCS1 had a better prognosis, while patients with high expression of FTH1 and PRDX had a worse prognosis.
Except for PRDX6, FBXW7, and MAP1LC3A, the other 14 genes were upregulated in HNSCC. AURKA is abundantly expressed in a variety of cancer types (Chow et al., 2017) and plays an important role in the proliferation of ovarian cancer cells (Yuan and Huang, 2019). AURKA overexpression can also inhibit ferroptosis of cancer cells by inhibiting GPX4 (Gomaa et al., 2019). Not only that, the high expression of AURKA may also be a poor prognostic marker for adrenocortical carcinoma, renal clear cell carcinoma, hepatocellular carcinoma (HCC), lung adenocarcinoma (ADC), and mesothelioma (Du et al., 2021). BNIP3 can promote tumor growth by promoting necrosis and autophagy and is associated with poor prognosis of patients (Vakkila and Lotze, 2004;Rosenfeldt and Ryan, 2009;Petrova et al., 2015). BNIP3 has also been upregulated in lung cancer and is associated with poor prognosis (Petrova et al., 2015). It is generally believed that SLC7A5 is related  to tumor development, angiogenesis, and poor prognosis of cancer patients as an ferroptosis regulator involved in energy metabolism (Giglia et al., 2014). Overexpression of CISD2 in head and neck cancer can enhance the resistance of sulfasalazine to ferroptosis (Kim et al., 2018). CISD2, which was found to be significantly upregulated in lung ADC samples, stimulates cancer cell proliferation and survival through elevated reactive oxygen species levels and ultimately leads to poor prognosis . LINC00336 was found to be upregulated in lung cancer and acts as a competitive endogenous RNA to inhibit ferroptosis in lung cancer to induce tumor formation (Wang et al., 2019b).
PRDX6 has the ability to remove lipid reactive oxygen species (LOOH) and inhibit ferroptosis (Lu et al., 2019). PRDX6 is overexpressed in a variety of cancers and is involved in the tumor progression of different tumors, such as lung cancer (Ho et al., 2010), thyroid (Nicolussi et al., 2014), and colorectal cancer (Huang et al., 2018). BAP1 induces ferroptosis by inhibiting cystine uptake and glutathione synthesis based on SLC7A11 inhibition, thereby inhibiting tumor development (Zhang et al., 2018. Renal cell carcinoma lacking BAP1 has also been found to be associated with a poor prognosis (Kapur et al., 2013). As a tumor suppressor in the process of human carcinogenesis, Fbxw7 is often observed to be downregulated in a variety of human cancers and is associated with a poor prognosis (Welcker and Clurman, 2008). SOCS1 can induce cells to be sensitive to ferroptosis by regulating the expression of p53 target genes (Saint-Germain et al., 2017). SOCS1, which was identified as a tumor suppressor gene for hepatocellular carcinoma (HCC), was found to be upregulated in hepatocytes and was associated with a good prognosis (Khan et al., 2020). The conclusions of these studies are unexpectedly consistent with our prognostic results. Currently, the role of ferroptosis in HNSCC is unclear. Although the function of many genes in cancer is still unclear, our research provides a new perspective on their potential role in HNSCC. Not only that, it can be found in the PPI and correlation network of these genes that they closely affect each other, which means that they may affect the progress of HNSCC through a complex interaction network. Wang et al. (2019d) found that CD8 + T cells release cytokines to drive ferroptosis to kill tumor cells, including tumor necrosis factor and interferon γ; blocking them will eliminate the ferroptosis induced by T cells. Yee et al. (2020) also found that ferroptosis induced by neutrophils promoted tumor necrosis in the progression of glioblastoma. These studies have shown that the immune system can inhibit tumorigenesis in part by promoting ferroptosis of cancer cells. In the enrichment analysis of GSEA, it is found that our ferroptosis-related signatures may be involved in various immune-related BPs, CCs, and MFs, suggesting that ferroptosis is closely related to tumor immunity in HNSCC. For this reason, we further explored the interaction between tumor immunity and ferroptosis in HNSCC from the perspective of immune microenvironment including immune infiltrating cells and immune function. Costimulation was found to promote the proliferation and survival of CD8 + T cells and TRegs (Chen and Flies, 2013). Combining our results, we infer that the increased T cell costimulation in patients in the lowrisk group may promote the proliferation of CD8 + T cells and TRegs, resulting in an increase in their content in HNSCC. The increased CD8 + T cells release interferon γ to increase the response of interferon γ to drive more ferroptosis to kill tumor cells, inhibit the progression of HNSCC, and ultimately lead to a better prognosis. TFH has been proven to play an antitumor immune cell effect by assisting CD8 effector T cells or directly eliminating tumor cells as cytotoxic T cells (Tran et al., 2014;Shi et al., 2018). Therefore, the higher TFH in HNSCC patients in the low-risk group may also induce ferroptosis by assisting CD8 effector T cells exerting antitumor effects. In breast cancer, mast cells are almost only found around the tumor and seem to have cytolytic activity on tumor cells (Rovere et al., 2007). Therefore, it is regarded as an antitumor factor and is associated with a good prognosis (Dabiri et al., 2004), which seems to explain the presence of more mast cells in HNSCC.
Gene mutations, especially TP53 mutations, which account for more than 50% of tumor mutations, play an important  role in the occurrence and development of tumors. In our study, TP53 mutations were found to be widely present in 72% of HNSCC patients from CGA database. As a reflection of the number of mutations in tumor cells, the impact of TMB on the prognosis of tumor patients is still controversial. TMB is considered a poor prognostic marker in neuroblastoma, resectable pancreatic cancer, and diffuse glioma (Elisa et al., 2018;Hwang et al., 2018;Luo et al., 2020), while it is a good prognostic marker in non-small cell lung cancer (Devarakonda et al., 2018). In HNSCC, we found that the TMB of the high-risk group was significantly higher, which may explain that patients in the high-TMB group may have fewer ferroptosis and lead to a worse prognosis. In our study, it was found that the TP53 mutation group had a worse prognosis, which was also supported by many related studies. In many sporadic cancers, TP53 mutations have also been observed to be associated with poor prognosis (Kandoth et al., 2013). In addition to tumor suppressor function, TP53 can also mediate the process of ferroptosis and inhibit tumor growth by regulating cell cystine metabolism and ROS response (Le et al., 2017). In contrast, TP53 mutations can inhibit ferroptosis by altering cellular iron acquisition and metabolism, thereby promoting cancer progression . We observed in HNSCC that patients in the TP53 mutation group had a higher risk score and were associated with a worse prognosis. This suggests that TP53 mutations in HNSCC may also inhibit ferroptosis, promote tumor progression, and lead to poor prognosis.
Keynote042, 024, and 021 based on pembrolizumab, Checkmate057 based on nivolumab, and IMpower131 and 150 based on atezulizumab all show that patients with high expression of PD-L1 can benefit more in immunotherapy. At present, the efficacy evaluation of PD-L1 detection for immunotherapy treatment of NSCLC has obtained the latest NCCN level 1 recommendation. In view of the conclusion that ferroptosis contributes to chemotherapy and immune checkpoint blocking therapy, we evaluated the possibility of prognostic models as markers for chemotherapy and immune checkpoint blocking treatment effect. From the corresponding results, it is found that our prognostic model can be used as expression value of PD-L1, PD-1, CTLA4, CD96, and TIGIT and six commonly used chemotherapeutic drug sensitivity prediction markers, which shows that our model may be used to predict the corresponding immunotherapy and chemotherapy efficacy to minimize adverse reactions.
As a statistical prediction model of cancer-specific survival probability, nomogram can combine the influence of several independent clinical variables to provide a personalized prognostic model for each patient (Duan et al., 2016). Due to the heterogeneity of cancer cells, patients will have different disease states at different times and lead to different clinical outcomes (Simson et al., 2007). Therefore, it is very necessary to construct  a nomogram with high prediction accuracy to accurately predict the prognosis of HNSCC patients, which is conducive to precise treatment and improves the treatment benefits of patients. After multiple methods of verification, the nomogram formed by our risk score satisfies this requirement well.
There are many prognostic models focusing on HNSCC. She et al. (2020) screened 27 prognostic immune-related genes to construct a signature that can effectively predict the prognosis of patients (She et al., 2020). In this article, the researchers have initially explored the role that immunity may play in the progression of HNSCC through functional analysis. Lu et al. (2019) paid close attention to the relationship between tumor-associated macrophages and HNSCC. In this study, 10 prognostic-related differentially expressed marker genes of macrophage were screened to construct a model to predict the prognosis of HNSCC patients with higher sensitivity (Lu et al., 2019). The author not only discovered the key transcription factors that regulate marker genes of macrophages and their regulatory relationship but also found four areas related to tumor-associated macrophages function in HNSCC through GSEA enrichment analysis: intercellular matrix remodeling, tumor killing, metabolic reprogramming, and tumor immunerelated pathways (Lu et al., 2019). Chen et al. (2021) adopted the novel weighted gene coexpression network analysis to identify 22 immune-related core genes and further screened 3 prognostic-related genes to construct a prognostic model. Not only that, this model shows good prognostic ability in both the training cohort and the external validation cohort. Our study focuses on the role of iron death in HNCC. This work not only used three sets to verify the ability of the model to predict the prognosis of patients but also explored the possible role of ferroptosis in the progression of HNSCC from the perspective of the immune microenvironment and TP53 mutation.
Although we have constructed a prognostic model with good clinical application value, and explored the role of ferroptosis in the development and prognosis of HNSCC from multiple perspectives such as immunity and mutation, and many valuable conclusions are obtained, there are still many shortcomings. First of all, because tobacco, HPV, and alcohol have a synergistic effect on HNSCC tumor progression, it seems essential to include tobacco and alcohol status in their research. However, it is difficult for us to obtain these data from the database and correct the influence of these factors. What is more, we could not find a dataset with both the transcriptome data of these 17 DE-FRG and survival data in other databases for external verification. To make up for this shortcoming, we split the HNSCC dataset of TCGA into three sets for multiple verification of the conclusions obtained. Second, although we used IHC images from clinical specimens from the HPA database to verify the differential expression of 11 PR-DE-FRGs at the protein level, the differential expression of 6 PR-DE-FRGs is still unknown. Due to various limitations of realistic conditions, it is difficult for us to collect clinical specimens for experiments to verify these results and conclusions. The good news is that our model performs well in all three data sets, but the reliability of the conclusions will still be questioned by future research results. Finally, some of the conclusions we got are only reasonable inferences based on the previous research conclusions and the results of our analysis. It only provides a new perspective for future research, and it is difficult to be treated as a conclusive conclusion.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ supplementary material.

AUTHOR CONTRIBUTIONS
XF and LL designed the research. XF and YO analyzed the data. XF, HL, XZ, and MC prepared the figures. XF, LL, DY, LZ, and QL wrote and revised the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This research was funded by the National Natural Science Foundation of China (81660444).