Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma.

Oral squamous cell carcinoma (OSCC) is a life-threatening disease, associated with poor prognosis and the absence of specific biomarkers. Studies have shown that the ferroptosis-related genes (FRGs) can be used as tumor prognostic markers. However, FRGs' prognostic value in OSCC needs further exploration. In our study, gene expression profile and clinical data of OSCC patients were collected from a public domain. We performed univariate and multivariate Cox regression analyses to construct a multigene signature. The Kaplan-Meier and receiver operating characteristic (ROC) methods were used to test the effectiveness of the signature, followed by the expression analysis of human leukocyte antigen (HLA) and immune checkpoints. The Cox regression analysis identified 4 hubs from 103 FRGs expressed in OSCC that were associated with overall survival (OS). A risk model based on the 4 FRGs was established to classify patients into high-risk and low-risk groups. Compared with the low-risk group, the survival time of the high-risk group was significantly reduced. According to the multivariate Cox regression analysis, the risk score acted as an independent predictor for OS. The accuracy of the 4 FRGs risk predictive model was confirmed by ROC curve analysis. Moreover, the low-risk group had the characteristics of higher expression of HLA and immune checkpoints, a lower tumor purity, and a higher immune infiltration, indicating a more sensitive response to immunotherapy. The novel FRGs-OSCC risk score system can be used to predict the prognosis of OSCC patients and their response to immunotherapy.

Background

Oral squamous cell carcinoma (OSCC) is the most common type of oral cancer, with more than 300 000 newly diagnosed cancer cases worldwide in 2018, according to the GLOBOCAN database. The OSCC predilection sites include the tongue, the alveolar, the mouth floor, the lips, and the buccal mucosa. The OSCC is characterized by high recurrence, metastatic, and mortality rates, especially in patients with a late diagnosis. In recent years, multidisciplinary collaborative diagnosis and treatment has been proposed for OSCC treatment, which included chemotherapy, biological therapy, and radiotherapy. Although multimodality therapies can improve the prognosis, the 5-year overall survival (OS) rate of OSCC patients remains stable, at approximately 56%, and the posttreatment local recurrence and distant metastasis rates are 25% to 50%.[3,4] Therefore, there is a need to find effective prognostic biomarkers that could guide these management decisions.

Ferroptosis is a new form of programmed cell death, characterized by iron overload and lipid peroxidation that cause lipid reactive oxygen species (ROS) accumulation.[5,6] Numerous studies demonstrated that ferroptosis is involved in cancer initiation, progression, and suppression. For instance, the tumor suppressor gene p53 may modulate the susceptibility of cancer cells to ferroptosis in a cell type–specific manner.[8,9] Artesunate, a clinically approved drug, can selectively kill OSCC cells by inducing ferroptosis. Recently, Kotaro Sato et al found that non-thermal plasma exposure kills OSCC cells through a specific mechanism that depends on ample catalytic Fe (II). There are 2 mixed forms of programmed cell death that are caused by this treatment method, including apoptosis and ferroptosis, which suggest that ferroptosis might be closely related to OSCC occurrence. Moreover, ferroptosis is also associated with the efficacy of immunotherapy.

It is well known that CD8+ T cells generally induce tumor cell death through the pore-forming protein-granzyme and the Fas/FasL pathways.[13,14] However, a new study showed that CD8+ T cells, which are activated by immunotherapy, augment ferroptosis in tumor cells, which contributes, therefore, to the antitumor efficacy of immunotherapy. This discovery provided important evidence for the correlation between ferroptosis and antitumor immunity. Besides, Lang et al also found that interferon-γ, which is produced by CD8+ T cells, has a synergistic effect with radiotherapy-activated ataxia-telangiectasia mutated, on promoting lipid oxidation and ferroptosis in tumor cells. These evidences indicate that the induction of ferroptosis is expected to enhance the antitumor efficacy of immunotherapy; however, it is not currently clear whether the immunotherapeutic targeting of ferroptosis is effective for OSCC patients.

The purpose of this study is to explore potential diagnostic and prognostic markers of OSCC and investigate their biological functions through bioinformatics analysis. We successfully constructed a novel prognostic model for OSCC, focusing on 4 FRGs that are mainly involved in the biological processes (BP) of immunity and glycolysis. Finally, we also discussed the prediction of a prognostic model on the sensitivity to immunotherapy.

Materials and Methods

Data collection

The mRNA expression profiles and the corresponding clinical characteristics of 273 OSCC patients were obtained from The Cancer Genome Atlas (TCGA) database (https://gdc-portal.nci.nih.gov/). The HTSeq-FPKM files of 273 oral samples (254 tumors and 19 controls) were retained, including the tongue, alveolus, buccal mucosa, soft and hard palate, oral cavity, and lips. In total, 108 driver genes and 111 marker genes were downloaded from FerrDb database (http://www.zhounan.org/ferrdb/). The genes tested only in mice (12 genes) and multi-annotated genes in both groups (15 genes) were filtered out. As a result, a total of 192 FRGs were obtained. The accession numbers of all samples were included in Supplementary Table 1. The names of public domains and the direct Web links were listed in Supplementary Table 2.

Cox risk regression establishment

A univariate Cox regression analysis was used to filter the prognostic-associated factors that were closely related to the OS of OSCC patients. Then, we performed a multivariate Cox regression analysis with a stepwise regression analysis to construct a risk model. Finally, 4 FRGs-OSCC were enrolled in a risk Cox regression, forming a risk formula that was determined by a linear combination of the 4 genes’ expression levels and weighted with the corresponding regression coefficients from the stepwise Cox regression model. The risk score was defined as

According to the median of the risk score, we divided the OSCC patients into 2 groups: a high-risk and a low-risk group. The Kaplan-Meier (K-M) and ROC analyses were performed, based on the risk score, using the survival R package. Moreover, univariate and multivariate Cox regression analyses were used to analyze whether the risk score was an independent prognostic factor. A nomogram was also established, based on the risk score, the pathological stage, and the M stage, to obtain survival rates of patients at 1, 3, and 5 years. We also used the K-M analysis to test the prognostic value of 4 single genes. Furthermore, randomized sampling method was used to obtain 4 random genes. A random gene signature was finally constructed to verify the validity of the 4-gene signature.

Identification of differentially expressed mRNAs

The limma R package was used to screen the differentially expressed genes (DEGs) between high-risk and low-risk samples, according to the thresholds of|log2 (fold change)| > 2.0 and P < .05.

Functional enrichment analysis

We performed a Gene Ontology (GO) enrichment analysis and a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the clusterProfiler R package. The GO results were composed of 3 parts: the BP, the cellular component (CC), and the molecular function (MF). Functional categories were considered when P value is less than .05.

Estimation of immune and stromal scores

ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) is a novel algorithm, based on ssGSEA (single sample Gene Set Enrichment Analysis) for predicting the level of tumor tissues’ infiltrating immune and stromal cells, based on gene expression profiles.[21,22] Herein, this method was applied to estimate the immune and stromal scores, for each OSCC patient from the 2 risk groups using the estimate R package. The tumor purity was inferred according to the formula derived from Prof. Yoshihara’s research.

Estimation of immune cell type fractions

CIBERSORT is a deconvolution algorithm that is used to characterize the cellular constitution of complex tissues.[23,24] The LM22 gene signature contains 547 gene expression signatures that can distinguish 22 human hematopoietic cell phenotypes, including natural killer (NK) cells, T cells, myeloid subsets, B cells, and plasma cells. The CIBERSORT R package and the txt files of LM22 are available on the CIBERSORT Web site (http://cibersort.stanford.edu/). We used the CIBERSORT method and LM22 to compare the proportions of 22 infiltrating immune cell types between the 2 risk groups. For each sample, the sum of all estimates of immune cell type fractions was equal to 1.

Differential analysis of immunotherapy between the high-risk and low-risk groups

We checked the expression of HLA genes and immune checkpoint molecules between the high-risk and low-risk groups. The HLA signatures were downloaded from Nomenclature (http://hla.alleles.org/genes/index.html). The immune checkpoint molecules included CD27, CD274, CTLA4, HAVCR2, ICOS, IDO1, LAG3, PDCD1, PDCD1LG2, and TIGIT. Tumor Immune Dysfunction and Exclusion (TIDE) is a computational framework for predicting the response of the immune checkpoint blockade response. The TIDE scores of OSCC patients were downloaded from the TIDE Web site (http://tide.dfci.harvard.edu), based on the uploaded transcriptome profiles.

Statistical analysis

The statistical analyses were undertaken using R v3.6.1 and Bioconductor (https://www.bioconductor.org/). All statistical tests were bilateral and P < .05 was considered significant. Besides, the differences were considered significant if *P < .05, **P < 0.01, ***P < .001, or ****P < .0001.

Results

A prognosis prediction model based on FRGs in OSCC

Using the comprehensive genome annotation files obtained from NCBI (Supplementary Table 2), 10 000 genes that were expressed in the tumor samples were identified by analyzing the transcriptome data of oral cancer of the tongue, alveolus, buccal mucosa, soft and hard palate, oral cavity, and lips, while the genes not expressed were eliminated. Subsequently, 192 FRGs were sorted out from the FerrDb database (Supplementary Table 3). The results of the 2 programs were integrated and only the intersection was selected to increase the specificity. As a result, 103 FRGs were expressed in the TCGA-OSCC (Figure 1A, Supplementary Table 4).

Figure 1.

Identification of the prognostic ferroptosis-related genes in the TCGA database. (A) Venn diagram to identify the ferroptosis-related genes (FRGs) which are expressed in oral squamous cell carcinoma (OSCC) samples. A total of 103 FRGs were expressed in TCGA-OSCC. (B) Forest plots showing the results of the multivariate Cox regression analysis between gene expression and overall survival. FTH1, FLT3, CDKN2A, and DDIT3 were identified as hub genes.

TCGA indicates The Cancer Genome Atlas.

To evaluate the prognostic value of the OSCC 103 FRGs, a univariate Cox regression analysis was performed on these genes and 5 variables (FTH1, FLT3, CDKN2A, DDIT3, and BNIP3; all P < .01) were identified to be significantly associated with OS (Table 1, Supplementary Table 5). Furthermore, we applied to the 5 FRGs a multivariate Cox regression with a stepwise regression. As a result, 4 FRGs (FTH1, FLT3, CDKN2A, and DDIT3) were harvested in the Cox regression (P < .05; Table 2, Figure 1B, Supplementary Table 6) and a risk formula was constructed with the expression levels of 4 genes and the corresponding regression coefficients: risk score = −0.0171 × CDKN2A + 0.0111 × DDIT3 − 1.2269 × FLT3 + 0.0014 × FTH1.

Table 1.

Univariate Cox regression of prognostic-related genes for OS.

ID	HR	HR.95L	HR.95H	P value
FTH1	1.001832	1.000681	1.002985	.00181
FLT3	0.189899	0.061876	0.582804	.003689
CDKN2A	0.981546	0.96921	0.994039	.003896
DDIT3	1.013436	1.0042	1.022756	.004273
BNIP3	1.021741	1.005669	1.038069	.007842

Abbreviations: HR, hazard ratio; OS, overall survival.

P value < .01 was considered statistically significant.

Table 2.

Multivariate Cox regression analysis results.

Gene ID	Coefficient	HR	HR.95L	HR.95H	P value
CDKN2A	–0.01707	0.983071	0.970235	0.996058	.01078
DDIT3	0.011122	1.011184	1.00103	1.021441	.030778
FLT3	–1.2269	0.2932	0.92297	0.931413	.037485
FTH1	0.001357	1.001358	1.000198	1.002518	.021707

Abbreviations: HR, hazard ratio; OS, overall survival.

P value < .05 was considered statistically significant.

The prognostic value of the ferroptosis-related signature in OSCC

We calculated the risk score of each patient by the previously described formula and classified the patients into a high-risk (n = 126, score ⩾ 1.0917) or a low-risk group (n = 126, score < 1.0917), based on the median of the risk score. The risk score ranking, the survival status of OSCC patients, and the heatmaps of 4 FRGs’ expressions are shown in Figure 2A to C. The abscissas of each graph represented risk scores of OSCC patients. According to the results of the heatmaps, patients in the high-risk group had much higher expression of DDIT3/FTH1 genes and lower expression of CDKN2A/FLT3 genes. Based on the K-M curves of the TCGA data set, there was an obvious difference between the 2 risk groups (P = 1.005e-04; Figure 2D). Furthermore, the ROC curve confirmed the predictive capacity of the 4-FRG risk signature. The area under the curve (AUC) for the risk signature was 0.749, indicating a considerable predicting power (Figure 2E). Regarding the prognostic value of 4 single genes, the K-M curve results showed that patients with high expression of FLT3 and low expression of FTH1 have better OS (P < .05; Supplementary Figure 1A and B). However, DDIT3 and CDKN2A are relatively poor in predicting the prognosis of patients (P > .05; Supplementary Figure 1C and D). In general, the 4-gene predictive model could reflect the prognosis of patients more effectively than a single gene. In addition, we used the sample function in R software and obtained 4 random genes from 103 FRGs, namely, CYBB, EGFR, ATM, and EGLN2. Subsequently, we performed a stepwise multivariate Cox analysis and 2 FRGs (EGFR and ATM) were finally used to construct a random gene signature (Supplementary Figure 2A). The risk score of each patient was calculated according to the previous formula and all patients were divided into high- and low-risk groups based on the median value of risk scores (Supplementary Figure 2B, C, and D). The K-M curve showed that there was no significant difference in OS between the 2 groups (Supplementary Figure 2E) and the AUC was 0.596 (Supplementary Figure 2F), indicating that the random gene signature is not reliable in predicting the prognosis of OSCC patients. These results showed that the 4-FRG signature could have better prediction performance from another perspective.

Figure 2.

Establishment of risk scores model. (A) The distribution of patients’ risk score. (B) The distribution of patients’ survival status. (C) Heatmaps of 4 ferroptosis-related genes’ expression. (D) Kaplan-Meier survival curves of overall survival between high-risk and low-risk groups. The abscissa represents the survival time and the ordinate represents the survival rate; red represents the high-risk group and purple represents the low-risk group. (E) ROC curves showed the predictive efficiency of the risk score model.

AUC indicates area under the curve; ROC, receiver operating characteristic.

The 4-FRG risk model was independent of conventional clinical characteristics

The correlation between the ferroptosis-related risk signature and the clinical traits was analyzed. Although the ferroptosis-related risk score was significantly related to histological grade and T stages, there was no correlation between sex, pathological stages, M stages, N stages, and risk score (Figure 3A and B, Supplementary Figure 3, P < .05). Meanwhile, a stratified analysis of the clinicopathological characteristics was further carried out and the results showed that the risk score level was closely related to prognosis (Supplementary Figure 4 and 5).

Figure 3.

Relationships between the risk score and clinicopathological parameters in oral squamous cell carcinoma (OSCC). (A) Clinical correlation analysis between risk scores and histological grades. (B) Clinical correlation analysis between risk scores and anatomical sizes of the primary tumors (P < .05). (C) Univariate analysis with Cox proportional hazard model of the association between clinicopathological variables and overall survival (OS) of patients with OSCC. (D) Multivariate analysis with Cox proportional hazard model of the association between clinicopathological variables and OS of patients with OSCC.

To analyze the relationship between the OS, the clinicopathological factors, and the ferroptosis-related risk signature, a univariate analysis was performed. Subsequently, a multivariate analysis was used to explore the independent prediction of the FRGs’ signature. The results of univariate and multivariate analyses with the Cox proportional hazard model showed that pathological stage, risk score, and M stage were independent prognostic indicators for OS (Figure 3C and D). The ferroptosis-related risk model also had a high accuracy in predicting the patients’ 5-year survival rate (Figure 4A to C).

Figure 4.

Construction of the nomogram model. (A) The nomogram model for predicting OS of oral squamous cell carcinoma (OSCC), each variable axis represented an individual risk factor, and the line drawn upward was used to determine the points of each variable. Then, the total points would be calculated to obtain the probability of 1-, 3-, and 5-year OS rate plotted on the 2 axes below. (B) The calibration plots for predicting patients’ 1-, 3-, or 5-year OS. The closer the slope is to 1, the more accurate the prediction is. The results show that the model has high accuracy in predicting the 5-year survival rate of patients with OSCC. (C) ROC analysis to evaluate the prognostic value of risk score.

AUC indicates area under the curve; OS, overall survival; ROC, receiver operating characteristic.

Functional enrichment analysis of DEGs between the 2 risk groups

To explore the DEGs between the 2 groups,|log2FC| > 2 and P < .05 were set as the screening criteria using the limma R package. As a result, 631 DEGs were identified, including 244 upregulated and 387 downregulated genes (Supplementary Table 7).

The biological functions of the 631 DEGs were speculated by GO enrichment and KEGG pathway analyses with the help of the clusterProfiler R package (Supplementary Table 8 and 9). The top 10 enriched GO terms and 11 pathways were listed in Figure 5. It was revealed that the 631 DEGs were significantly enriched in metabolic processes, including glucose catabolic process to pyruvate (GO:0061718), canonical glycolysis (GO:0061621), NADH regeneration (GO:0006735), and glycolytic process through fructose-6-phosphate (GO:0061615) in the biological process category. As for the CC category, the top 4 markedly enriched GO terms were chaperone complex (GO:0101031), chaperonin-containing T-complex (GO:0005832), collagen-containing extracellular matrix (GO:0062023), and eukaryotic 48S preinitiation complex (GO:0033290). The most significantly enriched MF terms included alcohol dehydrogenase (NADP+) activity (GO:0008106), and alditol: NADP+ 1-oxidoreductase activity (GO:0004032) (Figure 5A). Also, 11 enriched pathway terms, including Glycolysis/Gluconeo-genesis (hsa00010), Glutathione metabolism (hsa00480) and Hypoxia-inducible factor 1 signaling pathway (hsa04066), were explored by KEGG pathway analysis (Figure 5B).

Figure 5.

Representative results of GO and KEGG analyses. (A) The top 10 most significant Gene ontology terms, including biological processes, the cellular component, and molecular function. (B) The result of KEGG pathway analysis.

GO indicates Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Notably, in GO-BP, we could observe that many DEGs were significantly enriched in the entry “neutrophil activation involved in immune response” (P = .0025), indicating these DEGs may participate in the immune response by regulating the activation of neutrophils.[28,29] In addition, GO-MF has also enriched “NAD or NADP as acceptor” and “alcohol dehydrogenase (NADP+) activity” items. It was important to know that the decomposing substance (NAD+) in NADH can increase the activity of macrophages to activate other cells of the immune system and, finally, stimulate the entire immune system.[30,31] This evidence indicated that the DEGs between the 2 groups may be closely related to immunity. Therefore, we turned our attention to the difference in the tumor microenvironment (TME) between the 2 risk groups.

Differential analysis of the TME between the 2 risk groups

We calculated the immune and stromal scores, for each OSCC patient from the 2 risk groups through the estimate R package, finding that patients from the low-risk group tended to have higher immune scores. However, the stromal score between the 2 groups showed no significant difference (Figure 6A to C). It was also shown that the patients from the low-risk group had a higher immune cells content in the TME.

Figure 6.

Comparison of the ESTIMATE scores between high-risk group and low-risk groups. (A) Immune score of high-/low-risk groups. Immune score represents the infiltration of immune cells in tumor tissue. (B) Stromal score of high-/low-risk groups. Stromal score captures the presence of stroma in tumor tissue. (C) Estimate score of high-/low-risk groups. Estimate score infers tumor purity. (D) TME cell composition group by high/low risk, including 22 human hematopoietic cell phenotypes. Adjusted P values were shown as ns, *P < .05, **P < .01, ***P < .001, ****P < .0001.

ESTIMATE indicates Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data; ns, not significant; TME, tumor microenvironment.

After that, we estimated the relative infiltration of the 22 immune cells in the 2 risk groups by the CIBERSORT R package and LM22. The relative proportion of 22 immune cells in OSCC samples from 2 groups is significantly different. The presence of resting memory CD4+ T cells, M0 Macrophages, and M2 Macrophages was positively correlated with the level of risk score. By contrast, the presence of naive B cell, activated memory CD4+ T cells, regulatory T cells, T follicular helper cells, CD8+ T cells, M1 Macrophages, gamma delta T cells, and resting mast cells was negatively correlated with the level of risk score (Figure 6D). The unique differences in immune infiltration between high- and low-risk groups suggested that the 4-FRG model can be used to predict the prognosis of immunotherapy.

Differential analysis of the immunotherapy response between the 2 risk groups

Immunotherapy targeting immune checkpoints cytotoxic T lymphocyte antigen 4 (CTLA4) or programmed death 1/programmed cell death ligand 1 (PD/PDL1) has been successfully applied as a first-line treatment of OSCC.[32-34] However, it is disappointing that the response rates of this therapy are low, which suggests that only a subset of OSCC patients respond to immunotherapy, and that the efficacy could be enhanced by determining the type of immune checkpoint inhibitors.[35,36] Besides, the major histocompatibility class 1 (MHC1) complex is essential for presenting endogenous cellular antigens to circulating T cells and in initiating specific anticancer immune responses. In that case, we further investigated the HLA expression and immune checkpoints between the 2 different groups. The patients from the low-risk group possessed significantly higher expression of most HLA compared with the patients from the high-risk group (Figure 7A), proving that a higher immune status closely correlates with the OSCC prognosis. The expression of the immune checkpoints IDO1, LAG3, PDCD1, and TIGIT significantly increased in the low-risk patients’ group (Figure 7B). Based on these results, we speculated that OSCC patients with lower risk score might be promising candidates for immune checkpoint inhibitors. Furthermore, the TIDE score in the high-risk group was much higher (Figure 7C), suggesting that compared with the high-risk patients, the low-risk patients may be more sensitive to immunotherapy. This also confirmed the previous conclusion that the patients with a lower 4-FRG signature risk score were more suitable for immunotherapy.

Figure 7.

Distribution of immunotherapy response markers in high-/low-risk groups. (A) The distribution of HLA genes of high-/low-risk groups was shown in the box plot. (B) The distribution of immune checkpoint molecules of high-/low-risk groups was shown in the box plot. (C) The distribution of TIDE score in high-/low-risk groups.

HLA indicates human leucocyte antigen; ns, not significant; TIDE, Tumor Immune Dysfunction and Exclusion.

Discussion

Ferroptosis is a new type of programmed cell death involving lipid metabolism, iron metabolism, and amino acid metabolism. The metabolic level of cancer cells is much higher than that of normal cells, so a large amount of ROS accumulates in the cells, which makes regulation of ferroptosis an effective way to kill tumor cells. Evidences have shown that multiple cancer-related genes are involved in ferroptosis. P53 achieves bidirectional regulation of ferroptosis by affecting key proteins such as SLC7A11 and GLS2.[38-40] Myc, essential for metabolic reprogramming of tumor cells, inhibits ferroptosis in tumor cells by activating lymphoid-specific helicase. In addition, ferroptosis is also closely related to the TME. In the fatty acid–rich TME, CD8+ T cells will take up more polyunsaturated fatty acids through CD36, resulting in the accumulation of lipid peroxides and the upregulation of ferroptosis, which is detrimental to the antitumor ability. Hsieh et al found that zero-valent-iron nanoparticles could promote ferroptosis in lung cancer cells and also affect macrophage M1/M2 polarization and the antitumor function of CD8+ T cells. Therefore, finding ferroptosis-related genes (FRGs) related to OSCC prognosis is of great significance for precise targeted therapy of the disease.

To identify potential diagnostic and prognostic markers of OSCC, we constructed a prognostic model based on FRGs. By applying univariate Cox regression and multivariate Cox regression analyses to 103 FRGs that were downloaded from databases, we found 4 FRGs (FTH1, FLT3, CDKN2A, and DDIT3) that were related to the OS of OSCC patients. The FTH1 protein expression was significantly upregulated in breast cancer cells and the epigenetic silencing of FTH1 and TFRC that is induced by estrogen, reduced liver cancer cell growth, and survival. FLT3 is a receptor tyrosine kinase that plays a crucial role in the development of hematopoietic progenitor cells. Furthermore, FLT3 genetic alterations occurred in up to 30% of cases with acute myelogenous leukemia and patients with FLT3 mutations have poor outcomes. CDKN2A is a tumor suppressor gene that was reported to be frequently altered in OSCC progression. CDKN2A low gene expression is associated with the recurrence of disease in oral cancer patients and could be used as a prognostic marker for OSCC. DDIT3 is an endoplasmic reticulum stress-responsive transcription factor that plays an important role in apoptotic execution pathways that are induced by the endoplasmic reticulum stress. The speckle-type POZ protein contributes to prostate cancer by targeting DDIT3. Besides, DDIT3 acts as a transcription factor that enhances the expression of TNFRSF10A and TNFRSF10B, resulting in the initiation of ER stress-mediated apoptosis in human lung cancer cells. In summary, all 4 FRGs have been reported to be closely correlated with various cancers.

The prognostic model was constructed based on these 4 FRGs. The patients were separated into high-risk and low-risk groups, according to the threshold of the median risk score. As mentioned above, the patients with higher expression of CDKN2A/FLT3 genes and lower expression of DDIT3/FTH1 genes were more likely to be in the low-risk group, suggesting that FRGs CDKN2A/FLT3 genes may work as tumor suppressor genes. This result is consistent with previous research results.[44,45,49] The prognostic value of the 4-FRG risk signature was evaluated using the K-M and log-rank methods. There were significant differences in the survival curves of patients in the 2 groups. The prediction capability of the specificity and sensitivity of the FRG risk model was assessed by calculating the AUC of the risk score. Moreover, we found that the ferroptosis-related risk score was an independent prognostic indicator for OS when considering conventional clinical characteristics. The results indicated that this risk score model is a firm prognostic tool that can be used to classify patients and guide future targeted therapies.

For further understanding of the biological functions of DEGs, between different risk groups, we performed functional enrichment analyses. The results showed significant enrichment in processes, including NADH regeneration, glucose catabolic process to pyruvate, canonical glycolysis, and glycolytic process through fructose-6-phosphate, in the biological process category. The NADH regeneration is a metabolic process that consumes NAD+ to generate a pool of NADH, which is important to the immune system. It is indicated that NAD+ promotes the differentiation of CD4+ T cell without antigens. Furthermore, without relying on antigen-presenting cells, NAD+ can also regulate the fate of CD4+ T cell. Pro-inflammatory stimuli induce the NAD activation of macrophages and dendritic cells, resulting in a metabolic switch toward glycolysis,[28,29] while inflammatory macrophages depend on NAD+ salvage, resulting from ROS-mediated DNA damages. For KEGG, 11 pathways, including Glycolysis/Gluconeogenesis, Glutathione metabolism, and the HIF-1 signaling pathway, were identified. Several types of cancer, including OSCC, highly depend on glycolysis for ATP generation. Zheng et al found that zeste homolog2 can regulate STAT3 and FoxO1 signaling in human OSCC cells and promote invasion and tumor glycolysis. Another study demonstrated that circMDM2 could promote the proliferation and glycolysis of human OSCC cells by acting as ceRNAs to sponge miR-532-3p. The HIF-1 signaling pathway, as a cancer therapy glycolytic target, is involved in the regulation of glycolysis at preclinical and clinical stages.[53,54] Our results indicated that the DEGs, between the 2 different groups, may affect OSCC progression by altering these immunity-related BP or metabolic pathways.

We also focused on investigating the difference in the TME between the 2 different risk groups. To this end, we explored the correlation between OSCC and TME. The immune score and stromal score calculated based on the ESTIMATE algorithm can help quantify the immune and stromal components in tumors. The immune scores and ESTIMATE scores of the low-risk group were significantly higher than those in the high-risk group; however, no significant differences were found in the stromal scores. A high fraction of gamma delta T cells, macrophages M1, B cell naive, CD4 memory activated T cells, CD8 T cells, regulatory T cells, follicular helper T cells, and mast cells resting mainly infiltrated the tumors of the low-risk OSCC patients. A recent study suggested that the cell density of the high parenchymal CD8+, at the invading tumor edge, was associated with improved OS, and therefore could be used as an independent favorable prognostic marker for OSCC. Moreover, OSCC patients with high levels of CD4+CD25highCD127low regulatory T cells (Tregs) were found to have a better survival probability compared to patients with lower Tregs. This result indicated that immune cells might have an important effect on the OSCC TME. What’s more, we analyzed the expression of HLA-related genes, important to the immune system, in the 2 different risk groups, and found that the expression of most HLAs was significantly higher in the low-risk group, demonstrating that higher immune status was related to the prognosis of OSCC. The HLA molecules on the surface of tumor cells can help T cells recognize new antigens to create opportunities for anticancer immune responses. The expression of the immune checkpoints IDO1, LAG3, PDCD1, and TIGIT significantly increased in the low-risk group. Foy et al found that OSCC tissues are characterized by a higher level of intratumor T cells, overexpression of PD-L1 and IDO1, and a higher score of response signature to pembrolizumab, suggesting the inhibition of IDO1 and PD-L1 may have good clinical significance for OSCC. Another research indicated that several immune checkpoint receptors (TIM3, LAG3, IDO, PDL1, and CTLA4) could be considered as biomarkers that reflect the immune status in OSCC patients’ TME during nimotuzumab therapy. T cells from peripheral blood mononuclear cells, which were collected from OSCC donors, possessed a high expression level of TIGIT. Moreover, TIGIT blockade can promote the in vitro proliferative ability and effective cytokine secretion capacity of CD4+ T cells and CD8+ T cells isolated from OSCC patients. These results provided support for the hypothesis that OSCC patients with lower risk score (patients with higher expression of CDKN2A/FLT3 genes and lower expression of DDIT3/FTH1 genes) might respond better to the IDO1, LAG3, PDCD1, and TIGHT inhibitors.

For the first time, we assessed the effects of FRGs on the prognosis of OSCC and constructed an effective prognostic model to reveal the involved BP. We also proved that this model can be used as a criterion for determining whether a patient is suitable for immunotherapy.

Conclusions

In summary, we constructed an effective prognostic model based on 4 FRGs and proved that it has an independent correlation with the survival time of OSCC patients, which can provide useful information for the prediction of OSCC prognosis. However, more data from in vivo/in vitro experiments and clinical trials are needed to elucidate the mechanisms between FRGs and tumor immunity in OSCC.

Click here for additional data file.

Supplemental material, sj-tif-3-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-tif-4-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-tif-5-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-tiff-1-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-tiff-2-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-10-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-11-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-12-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-13-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-14-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-6-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-7-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-8-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights

Click here for additional data file.

Supplemental material, sj-xlsx-9-bbi-10.1177_11779322221115548 for Comprehensive Analysis of the Significance of Ferroptosis-Related Genes in the Prognosis and Immunotherapy of Oral Squamous Cell Carcinoma by Junhao Yin, Jiayao Fu, Yijie Zhao, Jiabao Xu, Changyu Chen, Lingyan Zheng and Baoli Wang in Bioinformatics and Biology Insights