Identification of a Ferroptosis Gene Set That Mediates the Prognosis of Squamous Cell Carcinoma of the Head and Neck

Squamous cell carcinoma of the head and neck (HNSCC) is one of the six most common malignancies. HNSCC has both a high incidence and poor prognosis, and its prognostic factors remain unclear. Ferroptosis is a newly discovered form of programmed cell death that is iron-dependent. Increasing evidence indicates that targeting ferroptosis may present a new form of anti-tumor treatment. However, the prognostic value of ferroptosis-related genes (FRGs) in HNSCC is unclear. This study was designed to identify molecular markers associated with ferroptosis that influence prognosis in patients with HNSCC. We used HNSCC tumor and normal data from The Cancer Genome Atlas (TCGA) to identify prognosis-related FRGs. An FRG-based prognostic risk score was constructed, and its prognostic value for patients with HNSCC was evaluated using receiver operating characteristic curve (ROC) and nomogram analyses. The model was validated using the Gene Expression Omnibus (GEO) database. Univariate Cox regression analysis in patients with HNSCC revealed 11 FRGs that were significantly associated with overall survival (OS). We constructed a ferroptosis risk score model based on five genes and divided the patients into different risk groups based on its median value. Kaplan-Meier curve analysis showed that patients with a higher ferroptosis risk score had shorter OS (TCGA training set: P < 0.001, TCGA validation set: P < 0.05,GEO validation set: P < 0.001), and Gene Expression Profiling Interactive Analysis (GEPIA) further verified the relationships between these five genes and prognosis in patients with HNSCC. Multivariate Cox regression analysis showed that the risk score remained an independent predictor of OS after the exclusion of clinical confounders (HR > 1, P < 0.01). Significant differences in gene function enrichment analysis and immune cell infiltration status were identified between the two groups. The prognostic model can be used to predict the prognosis of patients with HNSCC. Moreover, the five FRGs may affect ferroptosis in HNSCC and thereby represent potential treatment targets. These results provide new directions for HNSCC treatment.


INTRODUCTION
Squamous cell carcinoma is a common malignant tumor of the head and neck (HNSCC). The incidence of HNSCC is trending upwards, with approximately 600,000 new cases diagnosed yearly (Torre et al., 2015;Shield et al., 2017;Solomon et al., 2018). In the past 40 years, despite continuous advances in radiotherapy equipment and surgical techniques, the 5year survival rate of patients with HNSCC has remained at approximately 50%, largely due to its high recurrence rate and the disability rate associated with radiotherapy and surgical treatment (Koyfman et al., 2019). Clinical studies have shown that adjuvant therapy for advanced HNSCC using platinumbased chemotherapy is effective and can significantly improve patient survival rates. However, because some patients are insensitive to chemotherapy or develop drug resistance after treatment, the overall 5-year survival rate associated with chemotherapy remains low (Vidal et al., 2017;Wang et al., 2020). The only molecular targeted therapy that is approved for the treatment of HNSCC is cetuximab, but its clinical response rate is modest (only 10-15%) (Vermorken et al., 2007;Hammerman et al., 2015;Mandal et al., 2016). The complex mutation landscape of HNSCC may be the main reason for the low response rate of targeted therapy, as most tumors are caused by multiple genetic drivers that have evolved further after multiple rounds of chemoradiotherapy. This highlights the importance of targeted therapy selection for tumors. Irondependent cell death, called ferroptosis, had been identified as a form of programmed cell death. A growing number of studies suggest that ferroptosis may be a new cancer therapy strategy, potentially avoiding the need to target complex and redundant molecular pathways (Latunde-Dada, 2017;Yu et al., 2017;Hirschhorn and Stockwell, 2019). In addition to the traditional ferroptosis inducer erastin , sorafenib also induces ferroptosis. Ferroptosis also participates in the CD8 + T cell-mediated antitumor effects induced by immunotherapy . In addition, some cancerrelated genes, including TP53 and ALOX12, may be important ferroptosis drivers (Jiang et al., 2015;Chu et al., 2019). In contrast, GPX4 and FSP1 negatively regulate ferroptosis (Imai et al., 2017;Bersuker et al., 2019;Doll et al., 2019;Seibt et al., 2019). Numerous studies note the significant enrichment of ferroptosis-related genes (FRGs) in HNSCC, suggesting that ferroptosis targeting may be a new approach for treating advanced HNSCC (Ye et al., 2019;Liu et al., 2020;Yu et al., 2020). However, further advances in ferroptosis-targeting strategies for HNSCC require a detailed understanding of FRG expression patterns in different types of HNSCC. We analyzed The Cancer Genome Atlas (TCGA) database transcriptome data of 506 tumors to comprehensively describe differences in FRG transcription levels in HNSCC, and to provide a theoretical basis for ferroptosis-targeting HNSCC treatments. The objective of this study was to analyze FRG transcriptional variations in HNSCC to elucidate FRG expression profiles and assess the association between FRGs and clinical prognosis. Increasing evidence indicates that ferroptosis is also involved in immunotherapy. We therefore compared the immune infiltration status of different ferroptosis risk groups. Understanding the level of ferroptosis in HNSCC may help to develop a theoretical basis for ongoing basic research in ferroptosis oncology, accelerate the development of therapeutic ferroptosis targets in HNSCC, and guide clinical research.

Collection and Processing of TCGA Cohort and GEO Cohort Data
Gene expression data of RNAseq counts of 546 patients with HNSCC were downloaded from UCSC Xena. Related clinical data and survival data were also retrieved from the UCSC Xena (Table 1). Gene expression data were normalized using the variance stabilization transformation method of the "DEseq2" R software package. Gene expression annotation information was downloaded from the Ensembl website. 1 As the downloaded data was logarithmic in value, it needed to be restored to its original count value prior to comparative analysis. We used the "Deseq2" R software package to analyze the genetic difference between HNSCC tumor tissues and normal tissues to identify differentially expressed genes (DEGs). RNAseq data of 97 additional tumor samples and clinical information was download from the Gene Expression Omnibus (GEO) database 2 and analyzed with the "limma" package for internal standardization. Gene sequencing data annotation information was downloaded from the Bioconductor with the R package "hgu133plus2 GPL570 platform."

Identification of FRGs and Immune Cell Infiltration Status
A table of FRGs, containing six datasets, was obtained from the FerrDb web portal. 3 We analyzed drivers, suppressors, and markers of ferroptosis in this study, and 259 genes were considered for further experiments (Supplementary Table 1).
Single-sample gene set enrichment analysis (ssGSEA) was used to evaluate the immune infiltration status of patients in different risk groups (Supplementary Table 2).

Construction of a Prognostic Model Containing FRGs and Verification of Its Predictive Ability
Differentially expressed genes (DEGs) between HNSCC and normal tissues were identified with the "Deseq2" R package in the TCGA cohort (padj < 0.05 and|log2FC| ≥ 1). Online software 4 and the "heatmap" R package were used to generate the Venn diagram and heatmap, respectively. The STRING online database was used for protein interaction network analysis, which was visualized with Cytoscape (version 3.7.1). Most plots in this study were drawn using ggplot2. To investigate the association of clinical pathologic characteristics and FRGs with overall survival (OS) in HNSCC, we randomly divided 466 patients with complete clinical information from TCGA into a training set (N = 326) and a validation set (N = 140) using a ratio of 7:3. In the training set, univariate Cox regression analysis was used to identify FRGs significantly related to prognosis, and P-values were adjusted with Benjamini and Hochberg analysis. To further identify those FRGs that were independent prognostic factors for HNSCC patients, a multivariate Cox regression analysis was used to exclude the effect of clinical confounders. We evaluated the impact of each factor on prognosis by calculating the hazard ratio (HR) and 95% confidence interval (CI). To further investigate the effect of differential FRGs on HNSCC prognosis, we constructed a risk score model including five FRGs. The risk score was calculated as risk score = (coefi × expi). The samples in the training set, validation set, and GEO validation set were divided into two different risk groups based on the median value of the risk score. Receiver operating characteristic (ROC) curve and nomogram analyses were used to assess the ability of the risk score to predict prognosis. Correlation analysis was performed using Pearson tests. A Kaplan-Meier curve for OS was generated, and the log-rank test was used to identify differences in OS between different risk groups. The "prcomp" function of the "stats" R package was used for principal component analysis (PCA). The Gene Expression Profiling Interactive Analysis (GEPIA) online database 5 was used to further validate differences in gene expression between cancer and normal tissues and the gene expression correlations with prognosis. 5 http://gepia.cancer-pku.cn/

Enrichment Analysis of DEGs
Based on the DEGs (identified with criteria padj < 0.05 and |log2FC| ≥ 1) identified from different risk groups, the R package of "clusterProfiler" was used to perform gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses. Differences in 16 immune cells and 13 immune-related pathways were analyzed in different risk groups using ssGSEA and the "gsva" R package.

Nomogram Construction and Evaluation
Based on the coefficients identified in the prognostic model, we developed a nomogram using the "rms" R package. The accuracy of the model prediction was determined by observing the difference between the observed nomogram and the ideal nomogram. The closer the observed nomogram was to the ideal nomogram, the higher the model prediction accuracy.

qRT-PCR
To further verify the expression levels of five signature genes in HNSCC tumor and normal tissues, we collected tumor tissues from patients with HNSCC who underwent oral squamous cell carcinoma (OSCC) resection at our hospital from January 2020 to March 2021. This study has been reviewed and approved by the Internal Review Committee of our hospital. Total RNA was extracted from HNSCC specimens and normal tissues using the RNA-Quick Purification kit (EScience, Shanghai, China). HiScript R II Q RT SuperMix (Vazyme) was used to synthesize complementary DNA (cDNA). The ChamQ TM SYBR R qPCR Master Mix kit (Vazyme) was used to detect mRNA levels of targeted genes. The human gene primers used were (5 -3 ): GAPDH F-primer: CTCCTGCACCACCAACTGCT, GAPDH R-primer: GGGCCATCCACAGTCTTCTG; SLC7A11 F-primer: TCTCCAAAGGAGGTTACCTGC, SLC7A11 R-primer: AGACTCCCCTCAGTAAAGTGAC; TRIB3 F-primer: TACCTGCAAGGTGTACCCC, TRIB3 R-primer: GGTCCGAGTGAAAAAGGCGTA; AURKA F-primer: GAGGTCCAAAACGTGTTCTCG, AURKA R-primer: ACAGGATGAGGTACACTGGTTG; CAV1 F-primer: GCGACCCTAAACACCTCAAC, CAV1 R-primer: ATGCCGTCAAAACTGTGTGTC; AKR1C3 F-primer: GTCATCCGTATTTCAACCGGAG, AKR1C3 R-primer: CCACCCATCGTTTGTCTCGTT.
GAPDH was used to calibrate the relative expression levels of target genes. The primers were synthesized by Synbio Technologies (Suzhou, China). 2 − Ct was used to evaluate relative target gene expression.

Statistical Analysis
Statistical analysis of qRT-PCR results was performed using Prism GraphPad software V8.0 (LaJolla, CA, United States). Gene expression differences between the two groups were calculated with Student's t-test. Correlations between the expression of AKR1C3, CAV1, AURKA, and TRIB3 in HNSCC tissues and clinicopathologic parameters were measured with the chi-square test of SPSS version 23.0. The criterion for determining significant differences was P < 0.05. Survival differences between different risk groups were assessed with Log-rank test. Univariate and Multivariate Cox regression analyses were used to evaluate the relationship between clinical factors and patients' OS. The ssGSEA score for immune cell infiltration status was evaluated in the two different risk groups by Mann-Whitney test. All statistical data analyses were performed using R software (version 4.0.3).

RESULTS
The overall study design is shown in Figure 1. RNAseq expression data from 466 patients with HNSCC was obtained from the TCGA database and 97 OSCC patients from GEO (GSE41613) FIGURE 1 | Flow chart of overall study design and analysis.
Frontiers in Genetics | www.frontiersin.org were ultimately analyzed. The specific clinical characteristics of the patients are shown in Table 1.

Identification of DEGs and Prognosis-Related FRGs in the TCGA Cohort
RNA expression data was downloaded from TCGA and included 502 HNSCC tumor samples and 44 normal samples. Gene expression values in HNSCC tumor and normal tissues were compared using the "DEseq2" R package. We obtained 4,795 DEGs with the criteria of |log2FC ≥ 1 and padj < 0.05. A series of FRGs (56/259, 21.6%) were significantly differentially expressed between tumor and normal tissues (Figures 2A,B). To further identify the roles of these FRGs in HNSCC, univariable Cox regression analysis was used to identify prognosis-related genes and string was used to perform protein interaction network analysis (Figures 2C-E). We then used the "igraph" R package to construct a core network to identify critical modules and to further screen hub genes. The results indicated that EGFR, SLC2A1, SLC7A5, SLC7A11, AKR1C3, and CAV1 are hub genes ( Figure 2F).

Construction of a Prognostic Model Including FRGs With the TCGA Training Set
To further explore prognosis-associated FRGs and ferroptosisrelated therapeutic targets in HNSCC, we downloaded 466 HNSCC samples with complete clinical information and mRNA expression data from the TCGA database. We then randomly divided the 466 samples into training (n = 326) and validation (n = 140) sets. Survival analysis using Kaplan-Meier and Univariate Cox regression analyses revealed that the expression levels of the 11 DEGs were significantly related to prognosis in the TCGA training set ( Figure 2C). Subsequently, stepwise multivariate Cox regression analysis was used to identify the optimal genes for predicting prognosis. Finally, based on the median value, we constructed a prognostic risk score model containing five genes (SLC7A11, TRIB3, AURKA, CAV1, and AKR1C3) ( Figure 3A). The model indicated that patients with lower risk scores had better prognosis than those with a higher risk scores ( Figure 3E). Risk scores were calculated as follows: Based on the median value of the risk score (Figure 3C), the TCGA training set was divided into high-(n = 163) and low-risk (n = 163) groups. Patients in the high-risk group had a shorter OS (high-risk group: 3-year survival rate: 0.4762, 95% CI: 0.3942-0.575; 5-year survival rate: 0.3953, 95% CI: 0.3075-0.508; low-risk group: 3-year survival rate: 0.768, 95% CI: 0.698-0.844, 5-year survival rate: 0.661, 95% CI: 0.561-0.779, P < 0.001) than patients in the low-risk group (Figure 3B). PCA showed that patients in the different risk groups were distributed into two clear clusters ( Figure 3D). Survival analysis indicated that a higher risk score was associated with higher mortality than a lower risk score. To evaluate the predictive efficacy of our risk prediction model, we performed time-dependent ROC curve analysis and found that the AUC values were 0.659, 0.727, and 0.68 at 1, 3, and 5 years, respectively ( Figure 3F). Heatmap analysis showed that high-risk group patients had higher TRIB3, CAV1, SLC7A11, AURKA, and AKR1C3 expression levels than did low-risk group patients ( Figure 3G). Detailed information about coefficient profile is summarized in Supplementary Table 3.

Risk Model Validation Using the TCGA Validation Set and the GEO Validation Set
To assess the robustness of the FRG prognostic model, model performance was assessed in the TCGA validation set and in the GEO validation set. A total of 237 patients were used for model validation. For the TCGA validation set 140 patients were divided equally into two different risk groups (Figure 4E). A similar analysis showed that, in the TCGA validation set, patients in the high-risk group had a significantly shorter OS time than those in the low-risk group ( Figure 4A). Additionally, higher risk scores were associated with a worse prognosis ( Figure 4C). Heatmap results showed that patients in the high-risk group had higher TRIB3, CAV1, SLC7A11, AURKA, and AKR1C3 expression levels than patients in the low-risk group ( Figure 4G). The GEO external validation results were consistent with the TCGA results (Figures 4B,D,F,H).

Verification of the Expression of the Five FRGs and Prognosis in the GEPIA Online Database
We used GEPIA to evaluate the expression levels of the five signature genes in HNSCC and to further determine the accuracy of the FRG-based prognostic model. Consistent with R analysis results, except for AKR1C3 and SLC7A11, expression levels of the remaining three signature genes were significantly higher in HNSCC tissues than in normal tissues (Figures 5A-E). CAV1 and TRIB3 expression significantly correlated with OS ( Figures 5I,J), while AURKA expression was closely related to patient disease-free survival ( Figure 5H). There is no significant correlation between the expression of AKR1C3 and the OS of HNSCC patients (Figure 5G). Although there was no significant correlation between SLC7A11 and prognosis ( Figure 5F), a large number of studies have suggested that SLC7A11 is associated with a poor prognosis. Further research is needed to clarify the relationship between SLC7A11 and the prognosis of HNSCC patients.

Calculating the Independent Prognostic Value of the FRG Prognostic Model and Constructing a Nomogram to Evaluate Patient Prognosis
To explore whether the FRG prognostic model is an independent prognostic factor for HNSCC, we used univariate and multivariate Cox regression analyses. Univariate Cox regression analysis revealed that the risk scores in the TCGA training set, TCGA validation set, and GEO validation set were significantly     Figures 6B,D,F). These results confirmed the independent ability of the risk model to predict prognosis in patients with HNSCC (Figure 6). Based on the multivariate Cox analysis, a nomogram was constructed to predict the 1-, 3-, and 5-year OS of patients with HNSCC with the TCGA training set ( Figure 7A). The calibration charts for 1-and 3-year OS showed good predictive value in the TCGA training set and the GSE41613 validation set (Figures 7B,C).

GO and KEGG Enrichment Analyses of the TCGA Training and Validation Sets
We used the "ClusterProfile" R package to conduct GO and KEGG enrichment analyses of identified DEGs to verify the role of ferroptosis-related pathways in different risk groups. Our results showed that key oxidative stress pathways involved in regulating ferroptosis pathways were enriched in both the TCGA training set and the TCGA validation set (Figures 8A-D). For example, heme binding and iron ion binding were enriched in both the TCGA training set and the TCGA validation set (padj < 0.05, Figures 8A,C). KEGG analysis of the TCGA validation set showed some immunerelated pathways were enriched, including pathways related to cytokine-cytokine receptor interactions, IL-17 signaling, and primary immunodeficiency ( Figure 8D). GO enrichment analysis of the TCGA training set and TCGA validation set revealed increased cytokine activity (padj < 0.05, Figures 8A,C).
Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the enrichment scores for different immune cell subsets, related functions, and pathways in different risk groups. Surprisingly, in the TCGA training set, B cells, CD8 T cells, induced dendritic cells (iDCs), mast cells, plasmacytoid dendritic cells (pDCs), T follicular helper (Tfh) cells, T helper 1 (Th1) cells, T helper 2 (Th2) cells, and tumor-infiltrating lymphocytes (TILs) were significantly differently enriched between the lowrisk and high-risk groups (padj < 0.05, Figures 9A,B). KEGG pathway enrichment analysis revealed that the enrichment of T cell co-stimulation significantly differed between the two risk groups (padj < 0.05, Figure 9B). The terms T cell co-inhibition, inflammation promotion, cytolytic activity, checkpoints, and CCR and APC co-inhibition were obviously enriched in the low-risk group (padj < 0.05, Figure 9B). TCGA validation set and GEO validation set analyses confirmed the differences in B cells, CD8 T cells, iDCs, mast cells, Th2 cells, TILs, cytolytic activity, and T cell co-stimulation between the two risk groups (padj < 0.05, Figures 9C-F).
IHC results showed that the CAV1, AKR1C3, and TRIB3 protein expression levels were significantly higher in HNSCC tumor tissues (P < 0.05), but the difference in AURKA was not significant (P = 0.144). This may be due to the small sample size (Figures 11A-H and Table 2).

DISCUSSION
We systematically calculated the differential expression of FRGs collected from a ferroptosis database in the cancer tissues of patients with HNSCC. We further analyzed the relationship between these differentially expressed FRGs and OS in patients with HNSCC to provide insights into the study of ferroptosis in HNSCC. Ferroptosis is a type of programmed cell death caused by imbalanced iron metabolism, which leads to a reduction in, or the disappearance of, mitochondrial crista and rupture of the mitochondrial outer membrane (Xie et al., 2016;Conrad et al., 2018). Interestingly, ferroptosis has been implicated in several different types of cancer, including lung cancer, melanoma, and hepatocellular carcinoma (Alvarez et al., 2017;Basit et al., 2017;Nie et al., 2018;Ashrafizadeh et al., 2019;Bai et al., 2019;Gai et al., 2020). Recently, an increasing number of basic studies have indicated that ferroptosis caused by iron overload plays an important role in reversing chemoradiotherapy resistance and improving the efficacy of immunotherapy (Sleire et al., 2015;Li et al., 2020). Here, Cox regression analyses were used to construct a ferroptosis-related prognosis model consisting of five genes (SLC7A11, AURKA, TRIB3, AKR1C3, and CAV1), and internal and external validation was conducted. GO and KEGG enrichment analyses of genes that were differentially expressed in the different risk groups showed significant enrichment in oxidative stress-related pathways associated with ferroptosis, as well as in immune-related pathways. These results are consistent with those of prior work. SLC7A11 is a wellknown ferroptosis-regulating factor, and the erastin ferroptosisinducer targets SLC7A11 to inhibit the cystine glutamate transport system to ultimately inhibit ferroptosis. Recent studies have shown that SLC7A11 expression is significantly lower in HPV-positive tumors than in HNSCC HPV-negative tumors, and that HNSCC with lower SLC7A11 expression is more sensitive to ferroptosis than HNSCC with higher SLC7A11 expression (Hémon et al., 2020). AURKA participates in tumor development by interfering with mitosis and other signaling pathways. AURKA is an important member of the AURK kinase family, which is involved in the occurrence of various tumors and is closely related to tumor prognosis (Bavetsias and Linardopoulos, 2015;Xie et al., 2017). Gomaa et al. (2019) found that AURKA was a direct target of miR-4715-3p and that inhibiting AURKA or modulating miR-4715-3p could inhibit GPX4 expression, leading to disruption of the intracellular lipid peroxidation clearance barrier and inducing ferroptosis.
TRIB3 is an important stress-related gene that can induce the expression of various stimulating factors and participates in the regulation of transforming growth factor β, phosphatidyl inositol 3 kinase, and SRC (the characteristic activated protein kinase signaling pathways) (Tomcik et al., 2016). The integrated stress response induces ATF4-dependent transcription during periods of negative feedback regulation. ATF4 promotes TRIB3 expression, and TRIB3 interacts with ATF4 and suppresses ATF4 transcriptional activity. TRIB3 can bind to Akt and inhibit its activation to block insulin signaling. TRIB3 may directly bind to and mask the "ThR-308" phosphorylation site in AKT1. TRIB3 can also interact with the P65 RELA NF-κB transactivator to inhibit phosphorylation. This inhibits the transcriptional activity of P65 RELA, leading to decreased ATF4 expression and increased erastin-induced cell death. ATF4 cooperates with HSPA5 to prevent GPX4 protein degradation and the induction of ferroptosis (Wu et al., 2003;Bi et al., 2015;Aimé et al., 2020). AKRs are a family of enzymes responsible for inhibiting the cytotoxic potential of aldehydes and ketones by converting them into corresponding alcohols. Thus, ferroptosis resistance may be related to the detoxification activity of AKRs, as it leads to a reduction in lipid peroxides, which are key operators of the ferroptosis process (Gagliardi et al., 2019). CAV1 protects hepatocytes from ferroptosis in autoimmune hepatitis (Deng et al., 2020). With the exception of SLC7A11, the mechanisms of the genes in ferroptosis are not clearly known. Therefore, whether these genes regulate ferroptosis and control the prognosis of patients with HNSCC is unclear. In 2019, Professor Wang et al. (2019) proposed, for the first time, that CD8 T cells partially induce tumor cell death by regulating ferroptosis , but the specific mechanism involved remains unclear. In this study, we conducted GO and KEGG enrichment analyses of the DEGs in the different risk groups and found that many immune cells and immune functions were enriched. This result confirms that ferroptosis is associated with immunity. Previous studies have suggested that CD8 + T cells can secrete IFNγ, promote the production of intracellular lipid peroxides, induce ferroptosis, and increase the efficacy of tumor immunotherapy. In this study, we found that the number of CD8 + T cells were significantly higher in the low-risk group than in the high-risk group, which is consistent with the present study. B cells have a variety of immune response functions. Tumor-infiltrating B lymphocytes can promote the T cell response and directly kill tumor cells by secreting immunoglobulins, thereby inhibiting tumor progression. B cells were significantly enriched in the low-risk group. These results indicate that B cells may also mediate ferroptosis to inhibit tumor progression. The number of iDCs was also significantly higher in the low-risk group. iDCs demonstrate effective antigen presentation following stimulation by cytotoxic T lymphocytes, thereby activating the immune response, which may also be related to ferroptosis. However, it is undeniable that these results still need further in-depth research for confirmation. There are several limitations to this study. First, we built a prognostic model based on the TCGA public database and verified it with the GEO database. The number of samples was limited, and most of the samples were from European and American white populations. There was a paucity of samples from Asians or other races. More data are needed to verify the clinical applicability of the risk model. Second, the inherent weakness of only considering factors related to ferroptosis to establish a prognostic model is inevitable because many significant prognostic-related genes in HNSCC may have been excluded. In addition, it is worth noting that the correlation between the ferroptosis risk score and immunity has not been verified at the cellular level. In summary, our study constructed a prognostic model based on five FRGs. This model was proven to be independently related to OS in both the training set and the validation set, providing a new molecular target for the prediction of HNSCC prognosis. The mechanism underlying the relationship between the FRGs and the infiltration status of immune cells in HNSCC is still unclear, and further research is needed.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The study was approved by the Internal Review Committee of The Third Affiliated Hospital of Kunming Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
CL and XW: study design, drafting the manuscript, interpretation of data, and statistical processing. CL, XW, RQ, and ZZ: acquisition and analysis of data. CS: reviewing and approving the final version of the manuscript. All authors read and approved the final manuscript.