The datasets used in the present study are all available on public databases. The RNA-sequencing (RNA-seq) expression data and the corresponding clinical information of HNSCC samples were obtained from the TCGA GDC portal (https://portal.gdc.cancer.gov/repository) and the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm. nih. gov/geo/). A total of 259 FRGs were downloaded from the ferroptosis database (FerrDb; http://www.zhounan.org/ferrdb) (Zhou and Bao, 2020). The validation data of mRNA expression and DNA copy number was retrieved from the Oncomine database (https://www.Oncomine.org/) (Rhodes et al., 2004). Besides, we also obtained immunohistochemistry (IHC) validation data from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

Based on the TCGA-HNSCC dataset, the “limma” R package was used to identify the differentially expressed FRGs in HNSCC tissues vs. adjacent non-cancerous tissues via Wilcoxon test, with a false discovery rate (FDR) < 0.01. The Bioconductor packages “clusterProfiler” and “enrichplot” were then used for gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the differentially expressed FRGs. Meanwhile, univariate Cox analysis of overall survival (OS) was performed to screen prognostic FRGs; the cutoff p-value was defined as 0.001. Next, LASSO Cox regression analysis was applied to construct a prognostic signature based on the above genes using the R package “glmnet.” Based on the established risk model, we calculated the risk score of each patient and identified the median risk score of all HNSCC samples. The risk score was calculated as follows: Risk score = ∑ j = 1 n C o e f j ∗ X j , where Coefj denotes the regression coefficient and Xj denotes the normalized expression level of each FRG (Miao et al., 2020). The patients were assigned to low or high-risk group according to the median risk score.

We performed principal-component analysis (PCA) based on the expression profiles of all genes and the prognostic signature genes in both TCGA cohort and GEO cohort using the “stats” R package. Next, we used the “survival” package in R to compare the OS between the low and high-risk groups and to plot Kaplan–Meier survival curves. Univariate and multivariate Cox regression analyses were used to determine whether the risk signature could act as an independent prognostic indicator. Additionally, the time-dependent Receiver Operating Characteristic (ROC) curves of clinical characteristics and the risk signature were drawn with the R package “survival ROC” and the Area Under the Curve (AUC) values at 1 year were calculated. Meanwhile, ROC curves of the risk signature at 1 year, 3 and 5 years were also drawn. The risk signature was then validated in the GSE65858 cohort.

In order to predict the prognosis of patients with HNSCC more accurately, we draw a nomogram with age, gender, stage, T stage, N stage, and risk score using “rms” R package. Meanwhile, time-dependent calibration curves were used to evaluate the accuracy of the predictive nomogram at 1, 2 and 3 years.

Next, to explore the signaling pathways related to the risk signature, we performed gene set enrichment analysis (GSEA) based on the model gene expression between the low and high-risk subgroups. We set the number of permutations as 1,000 and chose the top five results in each group to build an enrichment plot.

The stromal, immune, and ESTIMATE scores were compared between low and high-risk groups with the “ESTIMATE” R package. We then used the single-sample GSEA (ssGSEA) with “GSVA” and “GSEABase” R packages to evaluate some immune-related characteristics (including the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways) between different risk groups (Liang et al., 2020). Meanwhile, the abundance of 22 immune cells was estimated via CIBERSORT algorithm (https://cibersort.stanford.edu/) to further compare the different immune infiltration levels between low and high-risk groups. Besides, we analyzed the different expression levels of immune checkpoints including PD-1, CTLA4, LAG3, TIGIT, and BTLA between low and high-risk groups.

In order to explore the clinical significance of the constructed prognostic model for HNSCC treatment, “pRRophetic” R package was used to calculate the half-maximal inhibitory concentration (IC50) of common chemotherapeutic drugs in TCGA cohort. According to National Comprehensive Cancer Network (NCCN) guidelines Version 2.2021, Cisplatin, Paclitaxel, Docetaxel, Doxorubicin, Etoposide, Gemcitabine, Methotrexate, and Cytarabine were main chemotherapeutic agents for head and neck cancers. Besides, on the basis of previous studies (Tang et al., 2019; Gulati et al., 2020), IC50 of Gefitinib and Metformin was also explored in different risk groups. Then, we analyzed the correlations between the expression of prognostic risk genes and the resistance/sensitivity of pan-cancer cells to chemotherapeutic drugs based on the CellMiner database (https://discover.nci.nih.gov/cellminer), which is an open-access Web interface containing molecular and pharmacological data for the NCI-60 cancerous cell lines (a panel of 60 diverse human cancer cell lines) (Reinhold et al., 2012). And we totally chose 218 drugs approved by FDA from this database.

All statistical analyses were performed using R software (version 4.0.3; https://www.R-project.org) and Perl software (version 5.32.0.1-64bit; https://strawberryperl.com/). In this study, p < 0.05 was known as “statistically significant”, p < 0.01 was regarded as “more statistically significant” and p < 0.001 was taken as “most statistically significant”.

A total of 545 HNSCC samples, including 501 tumor samples and 44 normal samples, from the TCGA dataset were enrolled. Meanwhile, RNA-seq and clinical data of 270 HNSCC samples from GEO-HNSCC cohort was used as the validation dataset. Patients with a follow-up time <60 days and those with no survival information were excluded. Finally, 486 patients from TCGA-HNSCC cohort and 266 patients from GEO-HNSCC cohort were included in the analyses, the detailed clinical characteristics of whom were listed in Table 1. Besides, among 259 FRGs, 45 FRGs that were just tested in non-human species were excluded and a total of 187 genes were identified in the above two cohorts.

We filtered out 124 differentially expressed genes (DEGs) in HNSCC tissues vs. adjacent nontumorous tissues, 94 of which were upregulated and 30 were downregulated. As expected, GO and KEGG pathway enrichment analyses showed that these DEGs were mainly enriched in iron-related and metabolism-related molecular functions and ferroptosis-related and cancer-related pathways (Figures 1A,B). Moreover, univariate Cox regression analysis identified five differentially expressed FRGs related to OS of HNSCC, among which FTH1, BNIP3, TRIB3, and SLC2A3 were high-risk FRGs (p < 0.001, hazard ratio [HR] > 1) and CDKN2A was low-risk FRGs (p < 0.001, hazard ratio [HR] < 1) (Figure 1C). We excluded CDKN2A from the following research considering that CDKN2A was overexpressed in HNSCC samples, but it seemed unreasonable that up-regulated CDKN2A was a favorable factor for the OS of HNSCC samples according to univariate Cox analysis.

Based on the four high-risk FRGs obtained above, we constructed a ferroptosis-related prognostic model via LASSO Cox regression analysis and calculated the regression coefficient of each model gene (Supplementary Table S1). The 4-gene signature was constructed on genes including FTH1, BNIP3, TRIB3, and SLC2A3. The risk score of each patient was calculated according to the expression and regression coefficient of model genes. Patients in the training cohort were then divided into a high-risk group (n = 243) and a low-risk group (n = 243) based on the median risk score. In the validation cohort, patients were also classified into high-risk (n = 165) and low-risk (n = 101) groups.

The signature genes were validated using mRNA expression and DNA copy number data from the Oncomine database. The levels of mRNA expression and DNA copy number of FTH1, BNIP3, TRIB3, and SLC2A3 were all significantly elevated in HNSCC samples compared with those in normal samples (Table 2, Supplementary Figure S1), which were consistent with our results. Besides, the signature genes were also validated with IHC data from the HPA database (Figure 2).

The transcriptional levels of FTH1, BNIP3, TRIB3, and SLC2A3 were all significantly upregulated in the high-risk group compared to the low-risk group according to both training and validation cohorts (Supplementary Figure S2). Based on TCGA-HNSCC (Figures 3A,B) and GEO-HNSCC (Figures 4A,B) datasets, the PCA before and after establishment of the prognostic risk model indicated that patients in different risk subgroups were distributed into two directions well. The Kaplan–Meier survival curves showed the patients with HNSCC in the low-risk group had a significantly better OS than those in the high-risk group (TCGA-HNSCC cohort: Figure 3E, p < 0.001; GEO-HNSCC cohort: Figure 4E, p < 0.05). In the TCGA-HNSCC cohort, the 5-years survival rate of the high-risk group was 0.375 (95% CI: 0.300–0.469), while that of the low-risk group was 0.563 (95% CI: 0.477–0.663) (Figure 3E). Meanwhile, the 5-years survival rate of the high-risk group was 0.402 (95% CI: 0.290–0.559), whereas that of the low-risk group was 0.682 (95% CI: 0.562–0.829) (Figure 4E) in the GEO-HNSCC cohort. Apparently, patients in high-risk group had a lower 5-years survival rate than those in low-risk group. Univariate Cox regression analysis of OS indicated that several clinical characteristics, including clinical stage (p < 0.001), T stage (p < 0.001), and N stage (p < 0.001), as well as the risk score (TCGA-HNSCC: p < 0.001; GEO-HNSCC: p < 0.01), were effective prognostic indicators for patients with HNSCC (3C, 4C). Moreover, multivariate Cox regression analysis of OS demonstrated that the risk score was an independent prognostic predictor for HNSCC patients (TCGA-HNSCC: p < 0.001; GEO-HNSCC: p < 0.01) (3D, 4D). In order to explore the sensitivity and specificity of clinical characteristics and the risk signature with regard to survival prediction, we drew ROC curves and then calculated the AUC values (3F, 4F; Supplementary Figures S3A, S3B). The AUC values at 1 year in training and validation cohorts were 0.664 and 0.679, respectively, and the prognostic accuracy of this signature was higher than that of all six clinical characteristics (3F, 4F). Considering that HNSCC contains multiple tumors at different anatomical sites, we analyzed the correlation between anatomical sites and the risk score using Kruskal test (Supplementary Figures S4E) according to HNSCC-TCGA cohort. The risk score was shown to be significantly associated with different anatomical sites (p < 0.05).

The nomogram is an effective method to predict the onset, progression or prognosis of diseases by integrating multiple risk factors. We successfully constructed a nomogram based on seven risk factors, including age, gender, grade, stage, T stage, N stage, and the ferroptosis-related signature, to predict 1, 2, and 3-years OS in TCGA-HNSC cohort (Figure 5A). Each risk factor had its own point and all contributed to the total point of each patient, according to which we got to know the 1, 2, and 3-years OS probabilities of patients. Calibration curves indicated that the predicted 1, 2, and 3-years OS probabilities were all well consistent with the actual ones (Figures 5B–D).

The GSEA was performed to explore the active pathways or functions enriched in the low and high-risk groups according to TCGA-HNSCC cohort (Figure 6A) and GEO-HNSCC cohort (Figure 6B). Results of GSEA in the training and validation cohorts were basically the same and the detailed information was shown in Figures 6A,B. Briefly speaking, pathways enriched in the high-risk group were mainly energy metabolism-related, while among the low-risk group, the most enriched pathways were closely related to immunity.

To explore relationship between the prognostic risk model and immune status, the stromal, immune and ESTIMATE scores, immune cells, immune-related functions/pathways, and immune checkpoints were estimated in different risk groups. According to TCGA-HNSCC and GEO-HNSCC cohorts, the immune scores in low-risk group were remarkably higher than that in high-risk group (p < 0.001, 7A, 8A). Among GEO-HNSCC cohort, patients in low-risk group had higher ESTIMATE scores (p < 0.05, Figure 8C), indicating a higher level of tumor purity. Although there was no significant difference in ESTIMATE scores between the two risk groups in TCGA-HNSCC cohort (p = 0.064, Figure 7C), the low-risk group still had a higher mean value of ESTIMATE scores than the high-risk group. Next, the ssGSEA and CIBERSORT algorithm were used to evaluate the enrichment scores of immune cells and immune-related functions in different risk groups. We found that the infiltration levels of B cells, CD8⁺ T cells, mast cells, NK cells, iDCs, pDCs, helper T cells, follicular helper T cells, Tfh, Th1 cells, Th2 cells, TILs, M0 macrophage cells, eosinophils and activated mast cells were significantly different between low and high-risk groups in both TCGA-HNSCC (Figures 9A,C) and GEO-HNSCC cohorts (Figures 9B,D). Interestingly, the low-risk group was found to have a higher infiltration level of regulatory T cell (Tregs) as shown in Figures 9A,D, which might appear contradictory to the immunosuppressive nature of Tregs. Furthermore, compared with the high-risk group, the scores of check-point, cytolytic activity, HLA, MHC class I, T cell co-stimulation, T cell co-inhibition, and type II IFN response were all elevated in the low-risk group in both the training and validation cohorts (Figures 9A,B). To further evaluate the availability of the constructed risk model in immunotherapy, we explored immune checkpoints of HNSCC (Kok, 2020) in different risk groups and discovered that the expression levels of PD-1(PDCD1), CTLA4, LAG3, TIGIT, and BTLA were all significantly upregulated in the low-risk group compared with the high-risk group according to TCGA-HNSCC and GSE65858 cohorts (all p < 0.01, Figures 7D–H, 8D–H).

In order to investigate the correlation of the prognostic signature with efficacy of chemotherapy in HNSCC, we used IC50 to predict the treatment response to common chemotherapeutic drugs in TCGA cohort. We discovered that the IC50 of Cisplatin (p < 0.01, Figure 10A), Gemcitabine (p < 0.05, Figure 10F) and Cytarabine (p < 0.05, Figure 10H) was significantly higher in the low-risk group, whereas, the high-risk group had a higher IC50 of Paclitaxel (p < 0.05, Figure 10B), Doxorubicin (p < 0.001, Figure 10D) and Etoposide (p < 0.001, Figure 10E), which indicated that the risk signature could be an indicator for predicting sensitivity of chemotherapeutic drugs. However, there was no significant difference in the IC50 of Docetaxel (Figure 10C) and methotrexate (Figure 10G). In addition, the high-risk group was more sensitive to some novel anti-cancer drugs, such as Gefitinib (p < 0.001, Figure 10I) and Metformin (p < 0.001, Figure 10J). Based on the CellMiner database, we next explored the correlation between expression of signature genes and the resistance/sensitivity of pan-cancer cells to chemotherapeutic drugs. As a result, all signature genes were significantly associated with sensitivity of some chemotherapeutic drugs (p < 0.05, Supplementary Table S2). For example, the expression of FTH1 had a positive correlation with drug resistance of cancer cells to Arsenic trioxide, Tamoxifen and Raltitrexed (Figures 11A–C). Increased expression of BNIP3 was correlated with increased drug sensitivity of cancer cells to Cisplatin, Carboplatin and Gemcitabine (Figures 11D–F), whereas it was negatively correlated with drug sensitivity of cancer cells to Sunitinib, Palbociclib, and Trametinib (Figures 11G–I). Upregulated TRIB3 was associated with increased drug sensitivity of cancer cells to Imiquimod, Vismodegib and umbralisib (Figures 11J–L). As for risk gene SLC2A3, the expression had a positive correlation with drug sensitivity of cancer cells to Trametinib, but showed a negative correlation with drug sensitivity of cancer cells to Palbociclib and Carfilzomib (Figures 11M–O).