Identification of Ferroptosis-Related lncRNA Pairs for Predicting the Prognosis of Head and Neck Squamous Cell Carcinoma

Background Ferrogenesis was strongly associated with tumorigenesis and development, and activating the ferrogenic process was a novel regimen in treating cancer, especially conventional treatment-resistant cancers. The purpose of the article was to construct a ferroptosis-related long noncoding RNAs (FRlncRNAs) signature, regardless of expression levels to effectively predict prognosis and immunotherapeutic response for head and neck squamous cell carcinoma (HNSCC). Methods The RNA-seq data for HNSCC and corresponding clinical information were obtained in the TCGA database, and ferroptosis-related genes (FRGs) were extracted in the ferroptosis database. On this basis, differentially expressed FRlncRNAs (DEFRlncRNAs) pairs were identified through coexpression analysis, differential expression analysis, and a fresh pairing algorithm. Then, a risk assessment model was established with univariate Cox, LASSO, and multivariate Cox regression analysis. Finally, we evaluated the model from various aspects, including survival status, clinicopathological characteristics, infiltration status of immune cells, immune functions, chemotherapeutic sensitivity, immune checkpoint inhibitors (ICIs)-related molecules, and N6-methyladenosine (m6A) mRNA status. Result We established a signature of 11-DEFRlncRNA pairs related to the prognosis of HNSCC that had AUC values above 0.75 in the one-, three-, and five-year ROC curves, underscoring the high susceptibility and specifiability of predicting HNSCC prognosis. Survival rates were remarkably higher for the low-risk patients than for the high-risk patients, and the signature was significantly correlated with survival, clinical, T, and N stages. Finally, immune cell infiltration status, immune functions, chemotherapeutic sensitivity, and expression levels of ICIs-related and m6A-related molecules were statistically different among different groups. Conclusion Our study established a novel lncRNA signature, which is independent of specific expression levels, could predict patient prognosis, and might have promising clinical applications in HNCSS.


Introduction
Worldwide, the incidence and mortality of head and neck cancer are estimated at 930,000 and 470,000, respectively [1]. Head and neck squamous cell carcinoma (HNSCC) constitutes the majority of pathological types of head and neck cancers, and its major risk contributors are alcohol consumption, cigarette smoking, and human papilloma virus infection. Despite significant advances in treatment with surgery, radiation, chemotherapy, targeted therapies, and immunotherapy, the mortality rate of HNSCC remains high [2]. High insensitivity or resistance to chemotherapy is a major cause of death in patients with advanced HNSCC [3]. Significantly, ferroptosis inducers may be an effective weapon in the treatment of various chemotherapy-resistant tumors, including HNSCC [4].
Cancer cells are usually characterized by a defect in the cell death executioner mechanism, which is a main cause of resistance to treatment. Compared to normal cells, cancer cells have an increased need for iron in order to promote growth, and this dependence on iron leads to cancer cells being more susceptible to iron-catalyzed necrosis called ferroptosis [5]. Abundant evidence has found that ferrogenesis is strongly associated with tumorigenesis and development, and activating the ferrogenic process is a new treatment regimen for cancer, especially conventional treatment-resistant cancers [6][7][8]. Long noncoding RNAs (lncRNAs) are the most common RNAs, which are more than 200 nucleotides in length and have no protein-coding ability [9]. A study illustrated that LINC00336 functioned through interaction with ELAVL1 as a significant inhibitor of ferroptosis in oncogenesis [10]. Mao et al. indicated that P53RRA promoted ferroptosis in cancer by nuclear sequestration of p53 [11]. One recent study revealed that GABPB1-AS1 was upregulated by erastin, inhibiting peroxidase gene expression and accumulating reactive oxygen species and cancer cell death, indicating that GABPB1-AS1 might have an essential molecular function in ferroptosis with hepatocellular carcinoma cells [12]. Some other research has suggested that lncRNAs exert their antitumor effects by modulating ferroptosis [13,14].
Consequently, the identification of ferroptosis-related lncRNAs (FRlncRNAs) has significant implications for elucidating the specific mechanism of oncogenesis and predicting the prognosis of HNSCC. A variety of studies identified different FRlncRNAs and established signatures to predict patient prognosis with malignancies [15][16][17][18]. However, these promising signatures have a few intrinsic shortcomings that might not be clinically applicable for translation. ese signatures briefly incorporate transcriptomic data and clinical data based on sequencing or microarrays. However, because of varied platforms, detection, and batch technologies, the expression level of individual genes has great variability. In this scenario, the utilization of these models would be severely limited, and they would be prone to biased diagnostic results, so their diagnostic accuracy would not be sufficient to translate into practical applications [19,20]. Here, we built a 0-or-1 matrix based on differentially expressed FRlncRNAs (DEFRlncRNA) pairs and replaced particular transcriptome expression values with dichotomous aggregated values to eliminate bias in the transcriptome expression values obtained under diverse situations [21]. e study established a novel prognostic model on the basis of DEFRlncRNA pairs that is independent of expression level. We then assessed the model's predictive power, tumor immune infiltration, N6-methyladenosine (m 6 A) mRNA status, chemotherapeutic efficacy, and immune checkpoint inhibitors (ICIs)-related molecules. In conclusion, this signature can accurately predict patient prognosis and characterize diverse immune landscapes, which is a promising prognostic biomarker.

Collection of Data and Identification of DEFRlncRNAs.
e RNA-seq data of HNSCC patients and corresponding clinical characteristics were extracted from e Cancer Genome Atlas (TCGA (https://tcga-data.nci.nih.gov)). Tissue sources of HNSCC in the TCGA database include mainly the oral cavity, tonsils, pharynx, and larynx. By removing duplicate or severely missing data (unknown or 0-day follow-up time and unknown survival status), valid clinical data were obtained. On the basis of GTF files in the Ensembl database, we added annotations to these RNA-seq data and then obtained the expression profiles of mRNA and lncRNA [22]. Ferroptosis-related genes (FRGs, S1 Table) were obtained from the ferroptosis database (FerrDb; https://www. zhounan.org/ferrdb) and were applied to define FRlncRNAs through a coexpression strategy [23]. Some lncRNAs were recognized as FRlncRNAs by criteria of a p value smaller than 0.001 and cor greater than 0.4. Finally, we utilized the R package Limma to distinguish DEFRlncRNAs from FRlncRNAs, and statistical significance was assumed as FDR < 0.05 and |log2FC| ≥ 1.5.

Pairing DEFRlncRNAs.
e DEFRlncRNAs were separately paired, and a 0-or-1 matrix was created. If the first DEFRlncRNA expression is lower than the second one, it is scored as 0; otherwise, it is scored as 1. Since some pairs with no certain class do not precisely predict patient prognosis, these pairs were filtered out, when the amounts with an expression quantity of 0 or 1 made up less than 20% and more than 80% of all pairs, respectively.

Establishment of the Risk Model.
Univariate Cox, LASSO, and multivariate Cox regression analyses were utilized to construct a prognostic model. e concrete risk score on each patient was computed, and the risk score formula was as follows: the time-dependent receiver-operating characteristic (ROC) curves and the area under the curve (AUC) were obtained to assess the predictive ability of the model for survival. On the basis of the maximum inflection point of the Akaike information criterion (AIC) values in the five-year ROC curve, we used it to be a cut-off point for classifying patients into high-or low-risk groups. e ROC curve and decision curve analysis (DCA) were also used to assess the precision of the model compared to the traditional clinical features. e comparison between the established models was utilized to assess the forecasting ability of the novel model [16,[24][25][26][27][28].

Validation of the Prognostic Model.
e Kaplan-Meier analysis and log-rank test were applied to assess the difference in survival between different groups on the basis of the R package Limma. Next, the chi-squared test was applied to show the association between the signature and clinical characteristics. en, the Wilcoxon signedrank test was utilized to compute the risk score differences among different groups of these clinical characteristics. Also, univariate and multivariate Cox analyses were applied to confirm that the risk score was an independent predictor of clinical prognosis. Finally, we developed a nomogram integrating the signature and other clinical 2 Journal of Oncology characteristics to predict the one-, three-, and five-year survival rates of patients.

Evaluation of Immune Cell Infiltration.
To obtain accurate immune infiltration status, the widely accepted approaches to estimate immune cell infiltration were used, including TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC.

Exploration of the Risk Assessment
Model. e singlesample gene set enrichment analysis (ssGSEA) was utilized to evaluate the difference in immune function among diverse groups. en, we investigate the expression levels of ICIsrelated molecules and m 6 A-related genes and the half-inhibitory concentration (IC50) of chemotherapeutic agents in different groups. Figure 1 depicts a flowchart of the complete procedure. e RNA-seq data and clinical information contained 44 normal samples and 501 tumor samples in TCGA-HNSCC. Table 1 shows the clinical features of the patients. Next, according to GTF files from Ensembl, the data were annotated and coexpression analysis between FRGs and lncRNAs was done. Finally, 1344 FRlncRNAs were distinguished by the criteria of p value less than 0.001 and cor more than 0.4 (Table S2), and 196 FRlncRNAs were identified as DEFRlncRNAs by the standard of FDR < 0.05 and | log2FC| ≥ 1 (Table S3), among which 174 genes were upregulated and 22 genes were downregulated ( Figure 2(a)).

Construction of the Risk Assessment Model. 13,444
DEFRlncRNA pairs were acquired through the pairing of these DEFRlncRNAs. Next, 2,753 DEFRlncRNA pairs were distinguished as statistically significant through univariate Cox regression analysis (Table S4). Next, 19 candidate DEFRlncRNA pairs were determined via the LASSO regression analysis, whose coefficient profiles and a partial likelihood deviation plot are presented in Figures 2(b) and 2(c). Finally, we utilized multivariate Cox regression analysis to distinguish 11 DEFRlncRNA pairs and establish the prognostic model ( Figure 2(d)).
To evaluate the predictability of the model for survival, ROC curves, AUC, and DCA were assessed. e ROC curves indicated the high sensitiveness and specification of the 11-DEFRlncRNA-pair model for one-, three-, and five-year survival prediction, with all AUC values exceeding 0.75 (Figure 3(a)). Also, this signature has a larger AUC value compared to other established signatures ( Figure S1). Subsequently, in light of the maximum inflection point (0.302) of AIC values on the five-year ROC curve, we categorized patients into high-or low-risk groups (Figure 3(b)). e AUC and DCA for five-year survival prediction demonstrated that the signature was more precise than other traditional clinicopathological features (Figures 3(c) and 3(d)).

Assessment of the Prognostic Risk Model.
e profiles of risk score and survival status were displayed in the different groups of HNSCC patients, suggesting that the low-risk patients had better survival status than the high-risk patients (Figure 3(e)). With patients separated by the 11-DEFRlncRNA-pair signature, the low-risk patients showed substantially longer survival than the high-risk patients (Figure 3(f )). e association of the risk score of the model with the clinicopathological features of the patients was performed by chi-squared tests. e strip chart and scatter diagrams acquired indicated that survival status, clinical, T, and N stages were significantly related to the risk score (Figures 4(a)-4(e)). Univariate Cox analysis demonstrated that clinical elements, such as risk score (p value <0.001), age (p value<0.001), and clinical stage (p value<0.001), were significantly related to prognosis (Figure 4(f )), and multivariate Cox analysis demonstrated that risk score (p value<0.001), age (p value <0.001), and clinical stage (p value <0.001) were also individual prognostic risk elements (Figure 4(g)). Detailed information is provided in Supplementary Table S5. e nomogram combining clinicopathological characteristics and the 11-DEFRlncRNA-pair signature was reliable and accurate and can be utilized in the survival prediction of HNSCC patients ( Figure 5).  XCELL, and EPIC algorithms is presented in Figure 6. e findings revealed that the high-risk group was positively related to tumor-infiltrating immune cells, including B cells, CD4+ T cells, CD8+ T cells, NK cells, macrophages, monocytes, and myeloid dendritic cells. e detailed Spearman correlation analysis is presented in Supplementary Table S6.

Association between the Risk Model and Other
Biomarkers. Association between immune cell subsets and immune functions according to ssGSEA suggested that costimulation of APC, immune checkpoint response, cytolytic response, HLA, promotion of inflammation, costimulation of T cells, coinhibition of T cells, and type II INF response were observed as statistical differences among the different groups (Figure 7(a)). We discovered that the IC50 of methotrexate, gemcitabine, and docetaxel was statistical differences between the different groups; however, the difference in IC50 for cisplatin and paclitaxel was small (Figure 7(b)). We also explored if the model was associated with ICIs and revealed remarkable differences in the expression of CTLA-4 (p < 0.001), PDCD1 (p < 0.001), LAG3 (p < 0.01), TIGIT (p < 0.001), and BTLA (p < 0.001) among others, between the different groups (Figure 7(c)). e comparison of m 6 A-related mRNAs in different subgroups showed that the expression levels of METTL14 (p < 0.001), YTHDC2 (p < 0.001), RBM15 (p < 0.01), YTHDF2 (p < 0.01), YTHDC1 (p < 0.001), and HNRNPC (p < 0.05) were significant (Figure 7(d)).

Discussion
Ferroptosis is a unique pattern of cell death that has received extensive interest, especially in the area of tumorigenesis and therapies [29]. Several research studies have mainly concentrated on developing DEFRlncRNA signatures to assess the prognosis of tumor patients [16,18,[30][31][32]. Nevertheless, the vast majority of signatures are established on the basis of the expression levels of quantitative transcripts. e study was motivated by the idea of establishing an immunerelated lncRNA-pair signature and tried to establish a rational signature with two FRlncRNA combinations, which are not affected by their precise expression levels in the signature [33].
In the study, we developed for the first time a prognostic signature based on FRlncRNA pairs in HNSCC by structuring a 0 or 1 matrix. Moreover, the signature showed superior diagnostic accuracy and was predictive of patient prognosis, which has clinical application in HNCSS. e model can guide clinicians to choose appropriate chemotherapeutic agents for the treatment of HNSCC by comparing the sensitivity of the different groups to commonly administered drugs. e expression levels of ICIs-related molecules and m 6 A-related genes were statistically different in the different groups, which could provide a therapeutic theory for different immune and targeted therapies for HNSCC patients.
Ferroptosis is a novel mode of programmed cell death that may offer a novel way of antitumor treatment. Nonetheless, the current study has many shortcomings and limitations. Firstly, our sample size is relatively small and the normal to tumor sample counts are not proportional. Secondly, the results may be biased as the majority of samples from TCGA are nonmetastatic.
irdly, our signature needs further validation using external validation to be more convincing.
In general, our research identified a new lncRNA signature that was unaffected by expression level, which might be used to predict patient prognosis and has potential clinical uses in HNCSS. e putative regulatory mechanisms by which lncRNA regulates ferroptosis and how it impacts the therapeutic efficacy of ferroptosis inducers are likely to be explored in future investigations. We expect that the practicability of the model can be verified in further clinical studies.

Data Availability
e datasets generated and analyzed in the current study can be obtained in the TCGA repository at https://tcga-data.nci. nih.gov/tcga.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
Congxiao Wu designed the study and drafted the manuscript. Fei Liu and Haixu Chen performed the data analysis. Qin Liu, Chao Song, and Kang Cheng collected the data and conducted literature search. Zhiyong Gao and Cheng Fan supervised and revised the manuscript. All authors read and approved the final manuscript. Acknowledgments e data provided by the TCGA public database are greatly appreciated.

Supplementary Materials
Supplementary Figure 1: a comparison of one-, three-, and five-year ROC curves with other established models indicates the superiority of this risk model. Supplementary