A Novel Ferroptosis-Related Gene Signature to Predict Prognosis in Patients with Head and Neck Squamous Cell Carcinoma

The clinical TNM staging system is currently used to evaluate the prognosis of head and neck squamous cell carcinoma (HNSCC). The 5-year survival rate for patients with HNSCC is less than 50%, which is attributed to the lack of reliable prognostic biomarkers. Ferroptosis-related genes (FRGs) regulate cancer initiation and progression. Therefore, we analyzed the correlation between FRGs and the clinical outcomes of patients with HNSCC. A typical prognostic model of FRGs for HNSCC was constructed using bioinformatics tools and data from public databases, including The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and GeneCards. The model was generated based on the following six FRGs: ATG5, PRDX6, OTUB1, FTH1, SOCS1, and MAP3K5. The accuracy of model prediction was analyzed systematically. The overall survival (OS) of the high-risk group was significantly lower than that of the low-risk group. The AUC for 1-year, 3-year, and 5-year survival were 0.645, 0.721, and 0.737, respectively, in the training set (TCGA cohort) and 0.726, 0.620, and 0.584, respectively, in the validation set (GSE65858). The multivariate Cox regression analysis revealed that the risk score was an independent prognostic factor for HNSCC. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses revealed that six FRGs were enriched in the ferroptosis pathway. A novel FRG prognostic signature model was established for HNSCC. The findings of this study reveal that FRGs are potential biomarkers for HNSCC.


Introduction
Head and neck carcinoma (HNC), which is one of the top 10 tumors, accounts for 3% of all cancer cases. Each year, 900,000 new HNC cases and 500,000 HNC-related deaths are reported [1]. Head and neck squamous cell carcinoma (HNSCC) is the most common histological subtype of head and neck tumors [2]. Risk factors for HNSCC include smoking, drinking, and human papillomavirus (HPV) infection [3]. The current method for confirming the diagnosis of HNSCC remains pathologic histological examination [4]. Due to lack of indicators for early diagnosis, HNSCC is not easily detected and diagnosed, and approximately 60% of patients with HNSCC are at an advanced stage by the time they receive treatment [5]. The main treatment methods for focal or locally limited HNSCC are resection, radiotherapy, and systemic therapy [4]. Currently, the clinical TNM staging system is used to evaluate HNSCC prognosis. However, the predictive effect is not satisfactory, and the survival rate of patients with HNSCC is less than 50% [6]. With the advent of precision medicine, patients need personalized treatment, and studies of drug treatments targeting biomarkers have shown good results [7].
The genetic, biochemical, and morphological characteristics of ferroptosis are classified under the category of iron-dependent cell death and accumulation of superoxide lipids [8]. Ferroptosis has both advantages and disadvantages because it can promote tumor progression by upregulating DNA replication or inhibit tumor progression by enhancing cell death [9][10][11]. The induction of iron death has emerged as a promising therapeutic approach for triggering cancer cell death, particularly in malignancies that are resistant to conventional treatment [12,13]. Few studies have been conducted on iron-induced death in head and neck cancers, and even fewer studies have suggested that ferroptosis is associated with HNSCC pathogenesis [14,15]. Therefore, there is an urgent need to identify reliable biomarkers of ferroptosis.
The tumor microenvironment (TME) plays an important role in modern cancer treatment [16]. Tumorinfiltrating immune cells (TICs) are indispensable components of the TME [17]. In particular, the composition and activity of TICs at the tumor site are considered prognostic factors for many cancers [18]. The TME of HNSCC contains an increased number of immune infiltrates. The clinical efficacy of cancer immunotherapy against HNSCC has been demonstrated recently [19]. Previous studies reported that iron metabolism markedly influenced the TME and that the iron requirements of cancer cells were higher than those of nontumorous cells [20]. Therefore, the prognostic value of FRG in HNSCC must be evaluated.
In this study, correlations between the expression profiles of FRGs and clinical outcomes in patients with HNSCC were analyzed using clinical datasets curated in public databases. The Cox regression analysis was performed to analyze the prognostic value of six FRGs in HNSCC. Additionally, the correlation between these six FRGs and the TME of HNSCC was examined. The findings of this study reveal that FRGs are potential biomarkers for HNSCC.

Materials and Methods
2.1. Data Acquisition. RNA sequence and clinical information of 494 HNSCC cases were collected from TCGA database using the R package "TCGAbiolinks." A total of 44 normal samples and 8 samples with missing survival data were excluded. A validation dataset (GSE65858) containing information on 270 patients with HNSCC was downloaded from the GEO website (https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE65858) [21]. Four cases with missing tumor site information were excluded from the analysis; hence, only 266 patients with complete clinical data were included in our study. The keyword "ferroptosis" was used to identify 103 genes and download their data (S1) from the GeneCards database (http://www.genecards.org/) [22]. Since all data in our study were collected from public online databases, approval from the Ethics Committee or written informed consent from contributors was not required.

2.2.
Establishing the Prognostic Signature of FRGs. The 96 FRGs in HNSCC (S2) were identified from the intersection of the 103 FRGs from the GeneCards database and all genes from TCGA dataset, and the expression data of the 96 FRGs were subsequently combined with clinical data from 494 HNSCC patients. The univariate Cox regression analysis identified eight genes associated with survival of patients with HNSCC. Meanwhile, the multivariate Cox regression model identified six prognostic ferroptosis-related genes. The genes associated with the risk of developing HNSCC were determined based on hazard ratios (HRs) as follows: HR > 1, risk-associated genes; HR < 1, protective genes. Six FRGs were evaluated using a linear combination of the Cox regression coefficients (β). The best model was chosen based on the lowest Akaike information criterion (AIC) value. The risk score for all patients was calculated using Equation (1): where Coef ðiÞ represents the estimated coefficient of each FRG and x ðiÞ represents the expression of each FRG.
2.3. Evaluation of the Prognostic Model. All patients in the training cohort (TCGA cohort) were divided into high-risk and low-risk groups based on the cut-off values of median risk scores. Overall survival (OS) was evaluated using the Kaplan-Meier survival curves and compared using the logrank test. The risk curve and scatter plot of the risk score were generated using the R software package "pheatmap." Principal component analysis (PCA) was performed to visualize sample distribution. The R package "timeROC" was used to generate time-dependent receiver operating characteristic (ROC) curves, and area under the curve (AUC) plots were generated for the 1-year, 3-year, and 5-year survival rates to assess the predictive performance of the risk scoring model. Univariate and multivariate Cox regression analyses were used to assess the applicability of the prognostic model independent of other clinicopathological factors of patients with HNSCC, including age, sex, grade, clinical stage, Tstage, and risk score. N-stage and M-stage were not analyzed because these data were missing for some of the patients. To assess the net benefit to patients, we used the R package "stdca.R" for decision curve analysis (DCA) and plotted DCA curves for 1 year, 3 years, and 5 years. Finally, we validated model accuracy using the validation cohort (GSE65858).

Stratified Analyses of the Expression Levels of Six FRGs.
Differential expression of six ferroptosis-related genes was analyzed between high-risk and low-risk groups from TCGA cohort using the "ggpubr" R package.  Figure 2(a)). A prognostic model was constructed for six FRGs based on the expression levels of eight FRGs and the multivariate Cox regression analysis. The best prognostic gene signature was selected based on the lowest AIC value ( Table 2). Among them, ATG5, PRDX6, OTUB1, and FTH1 were classified as risk-associated genes (HR > 1), whereas SOCS1 and MAP3K5 were classified as protective genes (HR < 1) (Figure 2(b)). Kaplan-Meier survival curves were plotted based on expression levels (high-expression and low-expression) of six FRGs. The OS of patients with HNSCC in the ATG5, PRDX6, OTUB1, and FTH1 highexpression groups was lower than that of patients in the low-expression groups (Figures 2(c)-2(f)). Conversely, the OS of patients with HNSCC in the SOCS1 and MAP3K5 low-expression groups was lower than that of patients in the high-expression groups (Figures 2(g) and 2(h)).  Patients with high-risk scores exhibited decreased survival and increased death rates (Figure 3(i)). The prognostic model was validated using the GSE65858 cohort (n = 266). The OS of the high-risk group (n = 133) was significantly lower than that of the low-risk group (n = 133) (Figure 3(b)). The 1-year, 3-year, and 5-year OS rates of the high-risk group were 78.95%, 54.33%, and 36.81%, respectively, while those of the low-risk groups were 92.5%, 72%, and 51.50%, respectively. PCA was performed to demonstrate the significant risk score distribution differences between the low-risk and high-risk groups in the validation cohort (Figure 3

Stratified Analysis of Six FRGs Based on Clinical
Characteristics. The six FRGs were differentially expressed between the high-risk and low-risk groups. The expression of ATG5, PRDX6, OTUB1, and FTH1 was upregulated in the high-risk group (Figures 6(a), 6(b), 6(d), and 6(e)), whereas that of SOCS1 and MAP3K5 was downregulated in the high-risk group (Figures 6(c) and 6(f)). The correlations between the six FRGs and clinical parameters such as age, P16 status, HPV-ISH result, tumor grade, TNM stage, and N-stage were analyzed. The expression levels of PRDX6 in patients aged ≥ 65 years were higher than those in patients aged < 65 years (Figure 7(a)). Compared with those in the P16-positive group, the expression levels of SOCS1 and PRDX6 were upregulated in the P16-negative groups (Figure 7(b)). Additionally, the expression levels of SOCS1 and ATG5 in the HPV-ISH-positive group were higher than those in the HPV-ISH-negative group (Figure 7(c)). The expression level of ATG5 varied between G1 and G3 grades, while that of FTH1 varied between G1 and G2 grades, as well as between G1 and G3 grades. OTUB1 was differentially expressed between G1 and G4 grades, whereas PRDX6 was differentially expressed between G1 and G3 grades. The expression level of SOCS1 varied between G1 and G2 grades, as well as between G1 and G3 grades (Figure 7(d)). The expression level of SOCS1 in stage III tumors was higher than in stage I tumors (Figure 7(e)). The expression level of ATG5 varied between N0 and N1 stages, as well as between N0 and N2 stages. FTH1 and MAP3K5 were differentially expressed between N0 and N2 stages (Figure 7(f)).    (Figure 8(a)). The KEGG analysis revealed that ferroptosis was the most significantly enriched pathway (Figure 8(c)). The PPI network revealed correlations between MAP3K5, SOCS1, and PRDX6. FTH1, OTUB1, and ATG5 did not form part of the network (Figure 8(b)).

Correlation between Ferroptosis and TIC Infiltration in TCGA Cohort.
To investigate the correlation between ferroptosis and TIC infiltration, the "CIBERSORT" algorithm was used to analyze the relative proportion of immune cells in the top 22 HNSCC samples (Figure 9(a)). The association between the 22 TIC proportions and ferroptosis was repre-sented using a heatmap (Figure 9(b)). A scatter plot showing the association between the expression of FRGs and the proportion of 22 TICs in HNSCC samples revealed a positive association between ATG5 and resting CD4 memory T cells (Figure 10(a)). Additionally, resting dendritic cells, M1 macrophages, activated natural killer cells, and activated CD4 memory T cells were negatively correlated with FTH1, whereas macrophage M0 was positively correlated with FTH1 (Figures 10(b)-10(f)). Naive B cells and resting mast cells were positively correlated with MAP3K5 expression (Figures 10(g) and 10(h)). Furthermore, naive B cells, activated CD4 memory T cells, and CD8 T cells were positively correlated with SOCS1, whereas resting dendritic cells were negatively correlated with SOCS1 (Figures 10(i)-10(l)).

Discussion
The incidence of HNSCC is increasing worldwide [1]. The 5year survival rate of patients with HNSCC is less than 50% [6], which is attributed to lack of reliable prognostic biomarkers [25]. Recent studies have established a correlation between molecular markers such as autophagy genes, immune genes, autophagy-associated long noncoding RNAs (lncRNAs), and immune-related lncRNAs, and HNSCC prognosis [26][27][28][29], which may aid in determining clinical outcomes.
As ferroptosis is reportedly involved in both cancer progression and cancer suppression [30], it can be a novel therapeutic target for tumors. As such, prognostic ferroptosis-related signature genes have been established for various tumors, including hepatocellular carcinoma [31], glioma [32], uveal melanoma [33], and clear cell renal cancer [34]. Low concentrations of PTX and RSL3 synergistically suppress hypopharyngeal squamous carcinoma by inducing ferroptosis [15]. SLC7A11 is a biomarker and therapeutic target for HPV-positive HNSCC [14]. Additionally, ferroptosis enhances the clinical growth-inhibitory efficacy of PDT against oral tongue squamous cell carcinoma [35]. Previous studies have mainly focused on HNSCC pathogenesis, as well as on the effects of drugs and surgical treatment of HNSCC. However, the correlation between ferroptosis     Several studies have also suggested a close correlation between ferroptosis and tumor immunotherapy. Ferroptosis plays a critical role in the efficacy of tumor immunotherapy [36]. In cancer immunotherapy, CD8(+) T cells exert antitumor effects by promoting tumor ferroptosis [37]. Additionally, a chemically programmed vaccine targeting ferroptosis and immunity has been developed [38]. Thus, the ferroptosis pathway is a potential therapeutic target for tumors. However, further studies are needed to examine the role of ferroptosis in the immune environment of cancer. In this study, TCGA dataset was analyzed using the CIBERSORT algorithm to demonstrate the correlation between 22 TICs and FRGs.

Disease Markers
Six ferroptosis-related signature genes were subjected to GO and KEGG analyses, and a PPI network was constructed. The correlation between these genes and 22 TICs in HNSCC was examined. SOCS1 has been identified as an immune-related prognostic protective gene for HNSCC [39]. One study in an Indian population suggested that SOCS1 may mediate immunosuppression in HPV-positive tumors [40]. Treatment with resveratrol suppressed HNSCC proliferation through SOCS1 [41], suggesting that SOCS1 is a potential therapeutic target for HNSCC. This is in agreement with the findings of the present investigation, which demonstrated that the expression of SOCS1 was upregulated in HPV-positive patients in both the P16-positive and HPV-ISH-positive groups in TCGA cohort. The expression of SOCS1 was positively correlated with the proportion of naive B cells, activated CD4 T cells, and CD8 T cells. In contrast, the expression of SOCS1 was negatively correlated with the proportion of resting dendritic cells. These findings demonstrate the role of SOCS1 in HNSCC.
ATG5 has previously been reported to be associated with autophagy. Recent studies have reported that ATG5 regulates the migration, invasion, and apoptosis of HNSCC through autophagy [42,43]. Additionally, ATG5 determines the sensitivity of HNSCC to radiotherapy and chemotherapy [44][45][46][47]. The expression of ATG5 was upregulated in the HPV-positive group, indicating that it regulates radiotherapy sensitivity of HPV-positive tumors. GO analysis revealed that ATG5 was enriched in cellular components containing autophagosomes. The expression of ATG5 was upregulated in advanced-grade and N-stage tumors.
Radiotherapy exerts growth-inhibitory effects against HNSCC by promoting redox sensitivity through MAP3K5 [48,49]. The expression of MAP3K5 was downregulated in advanced N-stage tumors, which is consistent with the downregulated expression of MAP3K5 in the high-risk score group and increased survival rate in the MAP3K5 highexpression group. MAP3K5 was positively correlated with the proportion of naive B cells and resting mast cells.
PRDX6 suppresses apoptosis in HNSCC by exerting antioxidant effects [50]. The expression of PRDX6 was upregulated in the high-risk score group. Additionally, the PRDX6 high-expression group exhibited decreased OS. Furthermore, PRDX6 expression was upregulated in the subgroup with poor histopathological differentiation. These results suggest that PRDX6 is an oncogenic factor in HNSCC.
The correlation between HNSCC, OTUB1, and FTH1 has not previously been reported. OTUB1 is involved in ubiquitination in breast cancer [51]. GO enrichment analysis of six FRGs revealed that the terms biological processes and molecular functions also include ubiquitin-related mechanisms. Survival analysis revealed that the OTUB1 highexpression group exhibited low OS and that the expression of OTUB1 was upregulated in the high-risk score group. The results of this study suggest that OTUB1 is a risk factor for HNSCC, although further experimental studies are needed.
The National Center for Biotechnology Information database classifies FTH1 as a pseudogene. However, the role of FTH1 in tumor development has recently been reported. FTH1 functions as a neoplastic suppressor in non-small cell lung cancer [52], breast cancer [53], and ovarian cancer [54]. However, FTH1 functions as a tumor suppressor in metastatic melanoma [55]. In this study, FTH1 was identified as a risk factor for HNSCC. The OS of the FTH1 highexpression group was low. Additionally, the expression of FTH1 was upregulated in the high-risk score group.

Disease Markers
This study had some limitations. This was a retrospective study and established gene signatures based on data from public databases. Thus, the gene signature requires further validation in prospective studies and multicenter clinical trials. Additionally, the mechanisms underlying the association between ferroptosis-related genes and tumor immunity in HNSCC remain poorly understood and warrant further investigation.

Conclusions
This study identified and validated a clinical prognostic model based on six FRGs, which served as independent prognostic factors for patients with HNSCC. These genes may also serve as potential prognostic biomarkers for HNSCC.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
The authors have declared no conflict of interest. 14 Disease Markers