An oxidative stress-related prognostic signature for indicating the immune status of oral squamous cell carcinoma and guiding clinical treatment

Oral squamous cell carcinoma (OSCC) is the eighth most common cancer worldwide and presents high mortality. Oxidative stress, caused by reactive oxygen species accumulation, plays a crucial role in tumorigenesis, cancer progression, and drug resistance. Nevertheless, the specific prognostic and clinical values of oxidative stress-related genes (OSGs) in OSCC remain unclear. Here, we developed an oxidative stress-related prognostic signature according to mRNA expression data from The Cancer Genome Atlas (TCGA) database and evaluated its connections with the prognosis, clinical features, immune status, immunotherapy, and drug sensitivity of OSCC through a series of bioinformatics analyses. Finally, we filtered out six prognostic OSGs to construct a prognostic signature. On the basis of both TCGA-OSCC and GSE41613 cohorts, the signature was proven to be an independent prognostic factor with high accuracy and was confirmed to be an impactful indicator for predicting the prognosis and immune status of patients with OSCC. Additionally, we found that patients with high-risk scores may obtain greater benefit from immune checkpoint therapy compared to those with low-risk scores, and the risk score presented a close interaction with the tumor microenvironment and chemotherapy sensitivity. The prognostic signature may provide a valid and robust predictive tool that could predict the prognosis and immune status and guide clinicians to develop personalized therapeutic strategies for patients with OSCC.

Introduction

Oral squamous cell carcinoma (OSCC) accounts for most head and neck squamous cell carcinomas, and is the eighth most common cancer worldwide (Bray et al., 2018; Huang L. et al., 2021). Tobacco consumption is thought to be a major etiological factor of OSCC (Das et al., 2007). Although the progress in treatment techniques of OSCC has been notable in recent decades, the overall 5-years survival rate and recurrence rate (approximately 50%) remain disappointing (Nör and Gutkind, 2018). Therefore, it is crucial to develop an efficient and personalized therapeutic strategy for patients with OSCC. Recently, extensive efforts have been dedicated to the identification of prognostic biomarkers or signatures for OSCC based on gene expression or DNA methylation (Huang Z. et al., 2021; Huang et al., 2021b; Zou et al., 2021), whereas their specific roles in guiding personalized treatment still need to be explored in depth.

Oxidative stress is characterized by the imbalance between oxidant and antioxidant production, which contributes to an excess of reactive oxygen species (ROS) and can activate proto-oncogenes and inactivate cancer suppressor genes (Gorrini et al., 2013; Biselli-Chicote et al., 2019). ROS have been identified as a potentially critical factor in driving tumorigenesis and cancer progression (Qiu et al., 2020). Patients with OSCC have been shown to present an elevated level of oxidative stress and a compromised capacity of antioxidant defense (Choudhari et al., 2014). Oxidative stress can induce oxidative damage of DNA and protein, enhancing lipid peroxidation and antioxidant defense disorders, which, if unrepaired, can promote the formation and progression of oral cancer (Toyokuni et al., 1995; Choudhari et al., 2014). Additionally, ROS production is involved in the development of oral cancer in chewers of tobacco and areca nuts (Stich and Anders, 1989), and strictly correlates with the clinical stage in patients with advanced cancer (Mantovani et al., 2002). Furthermore, oxidative stress, as an additional metabolic feature, plays a pivotal immunoregulatory role in the tumor microenvironment (TME) (Cubillos-Ruiz et al., 2015; Maj et al., 2017). Previous studies have demonstrated that oxidative stress could not only alter the phenotype and function of myeloid dendritic cells (DCs) in the TME (Cubillos-Ruiz et al., 2015) but also control tumor Treg cell functional behavior and temper the therapeutic efficacy of immune checkpoint therapy (Maj et al., 2017). Importantly, aberrant levels of ROS can profoundly affect the tumor heterogeneity by modifying the DNA structure of cancer cells, frequently leading to chemotherapeutic resistance (de Sá Junior et al., 2017). Nevertheless, the specific roles of oxidative stress genes (OSGs) in the prognosis, immune status, and chemotherapy response of OSCC remain largely unclear.

In this study, we filtered out six prognostic OSGs to construct a predictive signature according to mRNA expression data from The Cancer Genome Atlas (TCGA) database. Then, the prognostic value of the signature and its connection with clinical features were thoroughly explored in TCGA-OSCC cohort and validated in an independent OSCC cohort GSE41613. Additionally, this signature was shown to have close connections with immune status, immunotherapy response, and chemotherapy sensitivity. Overall, our results demonstrate the potential roles of OSGs in the prognosis, immune status, and drug response of OSCC, and provide a reliable tool for predicting the prognosis of patients with OSCC and guiding clinical treatment.

Materials and methods

Raw data collection

The RNA-sequencing and somatic mutation data of 330 OSCC samples and 32 normal oral tissues with corresponding clinical information were acquired from the TCGA GDC portal (https://portal.gdc.cancer.gov/repository). Additionally, gene expression profiles and clinical information of 97 patients in OSCC were obtained from the GSE41613 dataset in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). A total of 1,399 protein domains related to oxidative stress were downloaded from the GeneCards database (https://www.genecards.org/), with a relevance score of ≥7 (Qiu et al., 2020). The immunohistochemistry (IHC) validation data was obtained from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

Identification of differentially expressed OSGs and functional enrichment analysis

Firstly, we removed the batch effect between TCGA and GEO cohorts using “ComBat” function of “sva” R package. Then, based on the training dataset from TCGA database, we used the “limma” R package to identify the differentially expressed OSGs in OSCC samples and para-cancerous oral tissues via the Wilcoxon test, with a |log₂fold change (FC)| > 1 and a false discovery rate (FDR) < 0.05. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were used to analyze the functions and pathways associated with the differentially expressed OSGs using the R packages “clusterProfiler” and “enrichplot,” with an FDR <0.05.

Construction of an oxidative stress-related prognostic signature

Based on the expression of differentially expressed OSGs filtered above and the survival information of patients from TCGA cohort, univariate Cox analysis of overall survival (OS) was applied to screen the prognostic OSGs via the coxph function of “survival” R package, with p < 0.01. Next, to minimize the risk of overfitting, we constructed a prognostic model based on the prognostic OSGs using LASSO Cox regression analysis using the R packages “survival” and “glmnet.” According to this oxidative stress-related gene signature, we calculated the risk score of each patient with OSCC on the basis of the regression coefficient and expression level of each model gene. The median risk score was regarded as a boundary to classify patients with OSCC into low-risk (LRisk) and high-risk (HRisk) groups.

Efficacy evaluation of the gene signature

We used the R packages “survival” and “survminer” to explore the OS difference between the LRisk and HRisk groups and to draw Kaplan–Meier (KM) survival curves. Next, to evaluate the predictive accuracy and sensitivity of the gene signature, the “survival”, “survminer” and “timeROC” R packages were used to plot time-dependent receiver operating characteristic (ROC) curves and the area under the curve (AUC) values were calculated using the additive Aalen model in “timeROC” function of “timeROC” R package. The independent prognostic value of the signature was analyzed via univariate and multivariate Cox regression analyses. Additionally, we explored the correlations between the prognostic signature and clinical traits via the Wilcoxon test.

Gene Set Enrichment Analysis in different risk groups

We performed GSEA to define the activated pathways in different risk subgroups according to the expression of the model genes. The annotated gene set “c2. cp.kegg.v7.4. symbol.gmt” was used for reference. The number of permutations was set as 1,000 and the top five pathways in each risk group were obtained to plot the enrichment results.

TME and immunotherapy analysis

We compared the immune scores, stromal scores, and tumor purity between the LRisk and HRisk groups using the “ESTIMATE” R package. Next, the single-sample GSEA (ssGSEA) was used to estimate the infiltration levels of 16 immune cells and the enrichment scores of 13 immune-related functions in different risk groups. The different expressions of human leukocyte antigen (HLA) genes were explored in TCGA-OSCC and GEO-OSCC cohorts. Although immune checkpoint inhibitors (ICIs) can provide long-lasting clinical benefits to patients with cancer, only a fraction of patients fully respond to immunotherapy (Jiang et al., 2018). There is evidence indicating that the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm and the tumor mutation burden (TMB) score can serve as predictive markers for the efficacy of ICIs (Jiang et al., 2018; Choucair et al., 2020). Therefore, to define the correlation between the prognostic signature and immunotherapy response to ICIs, we calculated the TIDE score of each patient with OSCC from the TCGA cohort online (http://tide.dfci.harvard.edu/) and analyzed the TMB score based on somatic mutation data from TCGA database. On the basis of TCGA dataset, we simultaneously calculated the RNA stemness score (RNAss) and the DNA stemness score (DNAss) based on the transcriptome data and DNA methylation data, respectively, to assess the tumor stemness of each patient with OSCC (Malta et al., 2018).

Drug sensitivity analysis

To analyze the correlation between chemotherapeutic drug response and the risk signature, we estimated the half-maximal inhibitory concentration (IC50) of chemotherapeutic agents in TCGA cohort based on drug response data from the Cancer Genome Project (CGP) (Garnett et al., 2012) and the two largest publicly available screening efforts, the Genomics of Drug Sensitivity in Cancer (GDSC) (Yang et al., 2013) and the Cancer Therapeutics Response Portal (CTRP) (Basu et al., 2013). The CGP-derived drug response data was downloaded from the website https://osf.io/5xvsg/and processed via the “pRRophetic” package (Geeleher et al., 2014). Moreover, we retrieved GDSC2 and CTRP-derived drug response data from the website https://osf.io/c6tfx/and analyzed them using the “oncoPredict” R package (Maeser et al., 2021).

To predict potential small molecular drugs that could reverse the gene expression of high-risk patients with OSCC, the differentially expressed genes (DEGs; FDR <0.05 & |log₂FC| > 1) between the high- and low-risk patients were uploaded to the Connectivity Map (CMap, http://www.broad.mit.edu/cmap/) and some small molecule drugs related to the risk signature were obtained online using a modified Kolmogorov–Smirnov test. A negative enrichment score indicated an inhibiting effect on the expression of high-risk genes. Finally, we set p < 0.05, enrichment <0, mean < −0.4 and percentage non-null ≥ 75 as the cut-off criteria (Cheng et al., 2021). The 2D chemical structures of the selected small molecule drugs were obtained from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/).

Statistical analysis

The drug sensitivity analysis via the oncoPredict R package was implemented with R software (version 4.1.1; https://www.R-project.org) and all other statistical tests were analyzed using Perl software (version 5.32; https://strawberryperl.com/) and R software (version 4.0.3). The differential analysis between two continuous variables was performed via the Wilcoxon test, and correlation analysis between two continuous variables was estimated using the Spearman’s correlation coefficient. All results were taken as statistically significant when p < 0.05.

Results

Identification of OSCC samples and OSGs

We included 459 OSCC samples in our study, including 362 from TCGA cohort (330 tumor samples and 32 normal samples) and 97 from the GSE41613 cohort. After excluding patients with no survival information and follow-up time <30 days, 322 patients in TCGA-OSCC cohort and 96 patients in GEO-OSCC cohort remained for further analyses. Detailed clinical information of the included patients is shown in Supplementary Table S1. Among 1,399 OSGs, 1,108 genes were present in both cohorts (Supplementary Table S2).

Establishment of an oxidative stress-related prognostic signature

Of the 1,108 OSGs, 239 DEGs were obtained in OSCC tissues vs adjacent non-cancerous tissues, including 136 upregulated and 103 downregulated genes (Supplementary Table S3). As expected, GO enrichment analysis demonstrated that these differentially expressed OSGs were generally enriched in oxidative stress-related biological processes (BPs) and cytokine-related molecular functions (MFs). Furthermore, KEGG pathway enrichment analysis indicated that they were mainly associated with cancer and inflammation-related pathways (Figures 1A,B). Next, we filtered out eight differentially expressed OSGs associated with OS using univariate Cox regression analysis, among which HPRT1, ADA, CCNA2, PLAU, IL1A, VEGFA, and CXCL8 were risk-associated OSGs (p < 0.01, hazard ratio [HR] > 1) and CTLA4 was a protection-associated OSGs (p < 0.01, HR < 1) (Supplementary Figure S1). Finally, an oxidative stress-related prognostic model was constructed using LASSO Cox regression analysis based on TCGA-OSCC cohort and the expression of each model gene was contributing to the risk score independently with different weight according to their corresponding regression coefficient (Coef). The risk score of each patient in both TCGA-OSCC and GEO-OSCC cohorts was calculated according to the expression values of the model genes and their corresponding regression Coefs (Table 1) with the following formula: Risk score = HPRT1*Coef_HPRT1 + ADA*Coef_ADA + PLAU*Coef_PLAU + CTLA4*Coef_CTLA4 + VEGFA*Coef_VEGFA + CXCL8*Coef_CXCL8. As a result, patients in TCGA-OSCC cohort were divided into LRisk (n = 161) and HRisk (n = 161) groups according to the median risk score, which was set as the cut-off value of risk score (0.499883). The validation cohort comprised 47 low-risk patients and 49 high-risk patients. Most signature genes (including HPRT1, ADA, PLAU, and VEGFA) were validated with IHC data from the HPA database (Figure 1C).

FIGURE 1

FIGURE 1. The functional enrichment analysis of oxidative stress-related differentially expressed OSGs between OSCC samples and matched adjacent normal tissues and the immunohistochemistry of signature genes from the HPA database. (A) GO term enrichment analysis of differentially expressed OSGs. (B) KEGG pathway enrichment analysis of differentially expressed OSGs. (C) Representative images showing the expression of HPRT1, ADA, PLAU, and VEGFA in OSCC tissues vs. normal oral cavity mucosal tissues. OSGs, Oxidative stress-related genes, OSCC, Oral squamous cell carcinoma, GO, Gene ontology, KEGG, Kyoto Encyclopedia of Genes and Genomes, HPA: Human Protein Atlas.

TABLE 1

TABLE 1. The prognostic model genes and their risk coefficients.

Prognostic evaluation of the gene signature in the training and validation cohorts

To test the prognostic values of the gene signature, we performed KM survival analysis, time-dependent ROC exploration, and univariate and multivariate Cox regression analyses in both the training and validation cohorts. As expected, patients in the HRisk group had a significantly worse OS than those in the LRisk group according to the KM survival curves of the two cohorts (TCGA: Figure 2A, p < 0.001; GEO: Figure 2B, p < 0.05). The accuracy of the gene signature in survival prediction was explored with time-dependent ROC curves and their corresponding AUC values (Figures 2C,D, Supplementary Figures S1B,C). The AUC values of the risk score reached 0.688 in TCGA-OSCC cohort and 0.709 in GEO-OSCC cohort at 3 years, which were both higher than those of the clinical variables (Figures 2C,D). To test the independent prognostic value of the model, we generated univariate and multivariate Cox regression analyses of OS. On the basis of univariate Cox analysis, the risk score in both the training and validation cohorts was significantly associated with OS (TCGA: Figure 2E, HR = 3.224, 95% confidence interval [CI] = 2.119–4.904, p < 0.001; GEO: Figure 2F, HR = 4.724, 95% CI = 1.793–12.541, p < 0.01). After adjusting for other confounding clinical factors using multivariate Cox analysis, the risk score was found to be an independent prognostic predictor of patients with OSCC (TCGA: Figure 2G, HR = 2.81, 95% CI = 1.79–4.413, p < 0.001; GEO: Figure 2H, HR = 3.779, 95% CI = 1.435–9.954, p < 0.01).

FIGURE 2

FIGURE 2. Prognostic analyses of the conducted signature in OSCC based on TCGA and GEO cohorts. The Kaplan–Meier survival curves of the prognostic signature in TCGA (A) and GEO (B) cohorts. The ROC curves and AUC values of the signature and clinical features in TCGA (C) and GEO (D) cohorts. Univariate Cox regression analyses of the signature and clinical features in TCGA (E) and GEO (F) cohorts. Multivariate Cox regression analyses of the signature and clinical features in TCGA (G) and GEO (H) cohorts. OSCC, Oral squamous cell carcinoma, TCGA, The Cancer Genome Atlas, GEO, The Gene Expression Omnibus, ROC, Receiver operating characteristic, AUC, Area under curve.

Expression levels and clinical features underlying the prognostic signature

Among patients in the training and validation cohorts, we compared the expression levels of all model genes between the HRisk and LRisk groups. As expected, the expression levels of risk-associated model genes were all significantly upregulated in the HRisk group, while the protection-associated CTLA4 was notably overexpressed in the LRisk group on the basis of both TCGA-OSCC and GEO-OSCC cohorts (Supplementary Figure S2, all p < 0.001). Based on the optimal cut-off expression value, we performed survival analysis of each model gene using the “survival” and “survminer” R packages, which demonstrated that the expression levels of all model genes were significantly correlated with OS in TCGA-OSCC cohort Supplementary Figures S3A–F, p < 0.01), and similarly, except for VEGFA (Supplementary Figure S3L, p = 0.072), the expression of the other model genes was significantly correlated with OS in the GEO-OSCC cohort (Supplementary Figures S3G,F, p < 0.05). However, patients with lower expression of VEGFA had longer OS in the GEO-OSCC cohort, although the p-value was not significant, possibly due to the small number of patients with lower expression of VEGFA (n = 13). The above results suggest that the prognostic model genes could serve as potential therapeutic targets for patients with OSCC. According to KM survival analyses, we found that clinical stage, T stage, and N stage had significant prognostic values in OSCC. As shown in Supplementary Figure S4, patients with clinical stage I–II, T_1–2, or N₀ had significantly better OS than those with stage III–IV (p < 0.001), T_3–4 (p < 0.001), or N_1–3 (p < 0.01) in TCGA-OSCC cohort. Additionally, survival analyses in the GEO-OSCC cohort confirmed that the stage I–II group was significantly correlated with better OS compared to the stage III–IV group (p < 0.001) (there were no data on the clinical T and N stages of OSCC in the GEO cohort). Next, by analyzing the correlation between the clinical features and risk score via the Wilcoxon test, we confirmed that patients with stage III–IV or T_3–4 were significantly related to a higher risk score, while patients with stage I–II or T_1–2 understandably had a lower risk score in TCGA-OSCC cohort (Figures 3D,E, p < 0.001), which might determine the worse OS of patients in HRisk group. Moreover, patients with stage III–IV still presented a higher mean risk score than those with stage I–II in the GEO-OSCC cohort, although there was no significant difference in the risk score between patients with different clinical stages (Figure 3I, p = 0.058). Therefore, the identified signature could be involved in the occurrence and development of OSCC.

FIGURE 3

FIGURE 3. The correlations between clinical features and the prognostic signature. Different risk scores in OSCC patients with different age (A), gender (B), grade (C), clinical stage (D), T stage (E) and N stage (F) according to TCGA cohort. Different risk scores in OSCC patients with different age (G), gender (H), and clinical stage (I) on basis of the validation cohort. OSCC, Oral squamous cell carcinoma, TCGA, The Cancer Genome Atlas.

Functional enrichment analysis in the HRisk and LRisk groups

To explore the main functions enriched in different risk subgroups, we performed GSEA in both TCGA-OSCC and GEO-OSCC cohorts. The results indicated that pathways that were active in the HRisk group were mainly cell cycle-related, while those enriched in the LRisk group were mostly related to autoimmunity and cell adhesion in both the training and validation cohorts. The detailed results of GSEA are shown in Figure 4.

FIGURE 4

FIGURE 4. Gene Set Enrichment Analyses based on the conducted signature. (A) Active pathways in the HRisk and LRisk groups in TCGA cohort. (B) Active pathways in the HRisk and LRisk groups in GEO cohort. GEO: The Gene Expression Omnibus, TCGA, The Cancer Genome Atlas, HRisk, high-risk, LRisk, low-risk.

TME and immunotherapy response analysis

Next, we evaluated the correlation between the prognostic signature and the TME in TCGA-OSCC and GEO-OSCC cohorts. The ESTIMATE results demonstrated that patients in the LRisk group presented notably higher immune and stromal scores, but a lower level of tumor purity compared to those in the HRisk group (all p < 0.001, Figures 5E–J). Consistently, the risk score negatively correlated with both the immune score and stromal score, whereas it showed a positive relationship with the level of tumor purity (all p < 0.001, Supplementary Figure S5). Meanwhile, the tumor stemness of the patients with OSCC was estimated using DNAss and RNAss based on the TCGA dataset. Both DNAss (p < 0.001, Figure 5A) and RNAss (p < 0.01, Figure 5C) were obviously higher in the HRisk group compared to those in the LRisk group. Similarly, both DNAss (p < 0.001, Figure 5B) and RNAss (p < 0.01, Figure 5D) were positively related to the risk score. We then estimated the enrichment levels of immune cells and immune-related functions in different risk groups using ssGSEA to accurately assess their immune status. The LRisk group displayed significantly higher infiltration levels of some innate immune cells, including mast cells, neutrophils, natural killer (NK) cells, DCs, immature DCs (iDCs), activated DCs (aDCs), and plasmacytoid DCs (pDCs), as well as adaptive immune cells, including B cells, CD8⁺ T cells, helper T cells, regulatory T (Treg) cells, T helper 1 (Th1) cells, Th2 cells, T follicular helper cells (Tfh), and tumor-infiltrating lymphocytes (TILs). The LRisk group also presented higher enrichment scores of some immune-related functions, including cytolytic activity, checkpoint, promoting inflammation, type II IFN response, T cell co-stimulation, and HLA (Figures 6A,B). KM survival curves based on the optimal cut-off value demonstrated that immune cells (including DCs, aDCs, iDCs, mast cells, neutrophils, NK cells, B cells, and helper T and Treg cells) and immune functions (including cytolytic activity, checkpoint, type II IFN response, T cell co-stimulation, and HLA) showed a beneficial effect on the prognosis of patients with OSCC according to both the training and validation cohorts (all p < 0.05, Supplementary Figure S6), which may result in the better prognosis of the LRisk group. Considering HLA markers play a crucial role in anti-tumor immunity by driving antigen presentation (Anderson et al., 2021), we analyzed the differential expression of 24 HLA genes between two risk groups. The results demonstrated that the LRisk group had higher expression levels of most HLA genes according to both TCGA-OSCC and GEO-OSCC cohorts (Figures 6C,D). Overall, patients in the LRisk group showed more active immune activity, which may explain their better prognosis.

FIGURE 5

FIGURE 5. The tumor stemness and immune microenvironment between different risk subgroups. (A,B) Comparison of DNAss and RNAss in different risk subgroups. (C,D) The relationship between risk score and DNAss or RNAss based on TCGA cohort. (E–G) The immune score, stromal score and tumor purity between different risk groups in TCGA cohort. (H–J) The immune score, stromal score and tumor purity between different risk groups in GEO cohort. TCGA, The Cancer Genome Atlas, GEO, The Gene Expression Omnibus.

FIGURE 6

FIGURE 6. Analyses of immune cells and immune functions of the prognostic signature. Evaluation of infiltrating scores of 16 immune cells and activity of 13 immune-related pathways in different risk groups via ssGSEA according to TCGA cohort (A) and GEO cohort (B). Differentially expressed analysis of 24 HLA genes based on TCGA cohort (C) and GEO cohort (D). The p values were showed as: *p < 0.05; **p < 0.01; ***p < 0.001. ssGSEA: single-sample Gene Set Enrichment Analysis, TCGA, The Cancer Genome Atlas, GEO, The Gene Expression Omnibus, HLA, Human leukocyte antigen.

Next, to evaluate the predictive effect of the risk signature on the immunotherapy response, we explored the relationship between the risk score and the TIDE/TMB score. On the basis of the TIDE algorithm, the results revealed that patients in the HRisk group presented prominently lower TIDE scores (p < 0.001, Figure 7A), higher immune exclusion scores (p < 0.001, Figure 7B), and greater immunotherapy response (p < 0.001, Figure 7G). Meanwhile, the risk score was positively associated with the TIDE score (p < 0.001, Figure 7D) and negatively related to the immune exclusion score (p < 0.001, Figure 7E). To estimate the sensitivity of the risk score for predicting an immunotherapy response, a ROC curve was plotted using “pROC” R package and the AUC value was 0.766 (95% CI: 0.695–0.837, Figure 7H). Moreover, we found that the TMB score in the HRisk group was notably higher than that in the LRisk group (p < 0.01, Figure 7C), and the risk score had a positive relationship with the TMB score (p < 0.001, Figure 7F). In view of the above results, the risk signature is an appropriate predicting indicator of ICI treatment with high sensitivity.

FIGURE 7

FIGURE 7. Evaluation of the predictive effect of the prognostic signature on ICI treatment response. (A–C) The scores of TIDE, immune exclusion and TMB in different risk subgroups. (D–F) The relationship between the risk score and the score of TIDE, immune exclusion and TMB. (G) Predicted ICI treatment responses in different risk subgroups based on TCGA cohort using TIDE algorithm. (H) The ROC curve and AUC value to estimate the accuracy of the signature for predicting ICI treatment response. ICI, Immune checkpoint inhibitor, TIDE, Tumor Immune Dysfunction and Exclusion, TMB, Tumor Mutation Burden, TCGA, The Cancer Genome Atlas, ROC, Receiver operating characteristic, AUC, Area under curve.

Selecting appropriate chemotherapy drugs and uncovering potential small molecular drugs

To identify the relationship between drug response and the risk signature, we evaluated the differential drug response between HRisk and LRisk groups (p < 0.05 was considered significant) and performed the Spearman correlation test between the IC50 of selected chemotherapy drugs and the risk score (|R| > 0.2 and p < 0.05 were set as the threshold values). Eight common OSCC chemotherapy drugs (Cisplatin, Paclitaxel, Cytarabine, Docetaxel, Doxorubicin, Gemcitabine, Methotrexate, and 5-Fluorouracil) were selected for predicting the drug response on the basis of National Comprehensive Cancer Network (NCCN) guidelines Version 2.2021. Additionally, Gefitinib, an orally active selective EGFR (epidermal growth factor receptor) inhibitor, was explored based on previous studies (Brown et al., 2010; Tang et al., 2019). The detailed results of all selected drugs are shown in Supplementary Figure S7. Consequently, patients with low-risk scores showed a more sensitive response to CTRP-derived drugs (5-Fluorouracil and Gemcitabine, Figures 8I–L) and CGP-derived Gefitinib (Figures 8C,D), while patients with high-risk scores were more sensitive to GDSC2-derived Paclitaxel (Figures 8E,F) and CGP/GDSC2-derived Docetaxel (Figures 8A,B,G,H).

FIGURE 8

FIGURE 8. Drug sensitivity analyses on basis of TCGA cohort. The IC50 of some selected chemotherapeutic drugs, including CGP/GDSC2-derived Docetaxel (A–G), CGP-derived Gefitinib (C), GDSC2-derived Paclitaxel (E) and CTRP-derived drugs, 5-Fluorouracil and Gemcitabine (I,K) in different risk subgroups. The correlations between the risk score and some selected chemotherapeutic drugs, including CGP/GDSC2-derived Docetaxel (B–H), CGP-derived Gefitinib (D), GDSC2-derived Paclitaxel (F) and CTRP-derived drugs, 5-Fluorouracil and Gemcitabine (J,L). TCGA, The Cancer Genome Atlas, IC50, Half-maximal inhibitory concentration, CGP, the Cancer Genome Project, GDSC, Genomics of Drug Sensitivity in Cancer, CTRP, Cancer Therapeutics Response Portal.

To uncover novel drugs for treating patients with OSCC in the HRisk group, we analyzed the DEGs between the LRisk and HRisk groups using the CMap algorithm. CMap analysis was used to screen out 10 small molecule drugs (including etifenin, cortisone, sulfaguanidine, cyclopenthiazide, alcuronium chloride, xylometazoline, dextromethorphan, AH-6809, methylprednisolone, and oxybuprocaine) (Table 2), which were considered likely to reverse the expression of high risk-related genes and may be novel drugs showing anti-tumor effects in HRisk patients with OSCC. Furthermore, we obtained 2D chemical structures of the selected small molecule candidates online to facilitate future studies, which are shown in Supplementary Figure S8.

TABLE 2

TABLE 2. Potential small molecular drugs targeting high-risk OSCC patients.

Discussion

Oxidative stress, caused by ROS accumulation, plays a crucial role in cancer cells throughout the initiation, progression, metastasis, recurrence, and therapy of multiple tumors (Hayes et al., 2020). Accumulating evidence shows that the ROS levels not only correlate with tumor growth but also affect both the TME and the sensitivity of cancer cells to chemotherapeutic agents (Sosa et al., 2013; de Sá Junior et al., 2017). Additionally, oxidative stress-related gene signatures have been identified to be a reliable and efficient tool to predict the prognosis and progression of cancers (Qiu et al., 2020; Wu et al., 2021). Considering the above views, OSGs may represent valuable biomarkers for predicting the prognosis, immune status, and drug sensitivity of cancers and thus help clinicians define individual treatment plans for patients. However, the predictive effect of OSGs in OSCC remains unclear and is yet to be thoroughly investigated. Hence, in this study, we first filtered out 239 aberrantly expressed OSGs in OSCC and explored their potential functional pathways in OSCC. Next, we selected eight prognostic OSGs, according to which an oxidative stress-related signature with six OSGs was conducted. Moreover, the main functions active in different subgroups classified by the risk signature were identified using GSEA analysis. Finally, to further define the specific roles of OSGs in OSCC, we deeply evaluated the correlations between the oxidative stress-related signature and the prognosis, clinical features, immune status, immunotherapy, and drug sensitivity of OSCC.

The oxidative stress-related prediction signature comprised five risk-associated OSGs (HPRT1, ADA, PLAU, VEGFA, and CXCL8) and one protection-associated OSG (CTLA4), which were verified to be prognostic indicators in OSCC via univariate Cox regression and KM survival analyses. Moreover, these five risk OSGs were significantly upregulated in the HRisk group, while CTLA4, as a protection gene, was notably overexpressed in the LRisk group. HPRT1, a salvage enzyme involved in nucleotide production and recycling in cell cycle modulation, has been shown to promote proliferation and metastasis of head and neck squamous cell carcinoma (HNSCC) through direct interaction with STAT3 and has been implicated as a promising prognostic indicator and potential therapeutic target for HNSCC (Wang et al., 2021b). ADA, a housekeeping enzyme crucial in purine metabolism, makes a certain contribution to the regulation of inflammatory reactions and immune status (Bagheri et al., 2019), and its inhibitor has been shown to be obviously associated with reduced tumor size and decreased aggressiveness of cancer cells (Kutryb-Zajac et al., 2018; Bagheri et al., 2019). Additionally, Wang et al. reported that pre-operative serum ADA levels could be a reliable independent prognostic predictor of OSCC (Zhu et al., 2019). Previous studies demonstrated that the upregulated expression of PLAU (uPA), VEGFA, and CXCL8 (IL8) could promote the occurrence and progression of OSCC and might be important prognostic factors for patients with OSCC (Hundsdorfer et al., 2005; Sales et al., 2016; Magnussen et al., 2017; de Matos et al., 2019; Reyimu et al., 2021). The immune checkpoint inhibitory receptor CTLA4 has also been identified as a potential therapeutic target for OSCC (Bundela et al., 2014) and can enhance the therapeutic efficacy of anti-PD-1 immunotherapy in patients with HPV⁺ OSCC (Dorta-Estremera et al., 2019). Combining the results of the current and previous studies, we believe that these six OSGs could serve as reliable prognostic biomarkers and provide potential treatment targets for OSCC.

The prediction signature, composed of six prognostic OSGs with different weighting coefficients, was verified to be an effective prognostic indicator of OSCC according to the training and validation cohorts using univariate Cox regression and KM survival analyses. Additionally, patients with OSCC in the HRisk group had significantly shorter OS than those in the LRisk group. Multivariate Cox regression analysis indicated that the six-gene prediction signature was an independent prognostic predictor of OS of patients with OSCC. Furthermore, the results of ROC curves and AUC values validated the high prediction accuracy of this signature. In terms of the relationship between the risk score and clinical features, patients with tumor stage III–IV or T_3–4 were significantly associated with a higher risk score, which suggested that the risk score increases with the progression of OSCC, showing the poorer prognosis of patients with the higher risk score.

We next performed GSEA analysis in both the training and validation cohorts to define the specific signal pathways involved in the oxidative stress-related signature. The results demonstrated that cell cycle-related pathways such as the cell cycle, spliceosome, and base excision repair were apparently activated in the HRisk group. Moreover, oxidative stress-mediated ROS production exerts a key influence on cell cycle dysregulation by incorporating phosphorylation, ubiquitination, and receptor activation, which can contribute to aberrant cell proliferation and promote tumor progression (Verbon et al., 2012; Ahmad et al., 2020). Therefore, the activated cell cycle-related pathways may relate to the poor prognosis of patients with OSCC in the HRisk group. Meanwhile, autoimmunity-related pathways (i.e., autoimmune thyroid disease and systemic lupus erythematosus) were enriched in the LRisk group, illustrating that the patients with OSCC in the LRisk group might present an active immune state.

Growing evidences suggest that the TME plays a critical role in carcinogenesis, tumor progression, and survival, among which the immune microenvironment serves as a determinative factor (Hinshaw and Shevde, 2019). Thus, we next explored the relationship between the constructed signature and the immune status of OSCC. According to the ESTIMATE algorithm, a higher risk score was notably correlated with a lower immune score based on TCGA-OSCC and GEO-OSCC cohorts, suggesting that the high-risk score indicated an immune-suppressive state of patients with OSCC. Additionally, patients in the HRisk group displayed significantly higher tumor immune evasion scores compared to low-risk patients, and the immune evasion score showed a positive correlation with the risk score on the basis of the TIDE results in TCGA-OSCC cohort, further demonstrating the poor prognosis and immune-suppressive state of patients in the HRisk group. When we estimated the correlation between the prognostic signature and immune functions using ssGSEA, we found that compared to patients in the HRisk group, those in the LRisk group presented higher enrichment scores of cytolytic activity, promoting inflammation, the type II IFN response, and HLA. Of note, the abovementioned immune functions were all positively associated with a favorable prognosis in patients with OSCC according to both the training and validation cohorts via KM survival analysis, which contributes to the beneficial prognosis of low-risk patients. Rooney et al. reported that the cytolytic activity was associated with a modest but significant pan-cancer survival benefit and was connected to counter-regulatory immune responses (Rooney et al., 2015). Consistently, high cytolytic activity in tumor-free tongue tissue conferred improved prognosis in patients with tongue squamous cell carcinoma (Gu et al., 2019). Type II IFN (IFNγ) is thought to play a crucial role in cancer immunosurveillance, with the ability to promote anti-tumor immunity by increasing tumor immunogenicity (Dunn et al., 2006). Naturally, HLA markers are essential for antigen presentation and display a pivotal role in anti-tumor immunity by enhancing immunosurveillance and preventing immune escape. Additionally, high HLA class I expression in OSCC shows a significantly positive connection with a better prognosis (Koike et al., 2020; Anderson et al., 2021). Besides, it is worth emphasizing that various immune cells in the TME participate in anti-tumor immune responses as main components. Hence, we simultaneously calculated the enrichment level of immune cells in different risk subgroups using ssGSEA. The results showed that patients in the LRisk group displayed significantly higher infiltration of DCs, NK cells, mast cells, neutrophils, B cells, helper T cells, and Treg cells, each of which had a positive relationship with favorable prognosis in patients with OSCC. As the most potent antigen-presenting cells, DCs underpin the successful generation of anti-tumor immune responses by initiating and regulating innate and adaptive immune responses in the TME, and thus, targeting DCs is a promising strategy to improve the efficacy of current immunotherapies (Verneau et al., 2020). NK cells, with the potent ability to kill tumor cells, induce remodeling of the oral TME via IFN-γ and TNF-α, as well as prevent tumor growth and metastasis (Jewett et al., 2018). Furthermore, both CD103⁺ DCs and activated NK cells have been shown to have a favorable prognosis in OSCC (Hadler-Olsen and Wirsing, 2019; Xiao et al., 2019). Interestingly, Tregs were significantly connected to an improved prognosis of OSCC in our study, despite being recognized as immunosuppressive cells in numerous cancers. Among patients with OSCC, a high level of infiltrated Tregs has been proven to be notably associated with a lower frequency of lymph node metastasis and prolonged OS (Bron et al., 2013). Moreover, Hanakawa et al. suggested that Tregs in the TME may prevent tumor cell invasion and metastasis by inhibiting inflammatory processes (Hanakawa et al., 2014). In brief, the immune-related results in our study expectedly cohere with those in previous studies and the improved prognosis of patients in the LRisk group may be a result of the active immune state and increasing anti-tumor immune responses in these patients.

In addition to the TME, as the central component of cancers, tumor cells naturally play a vital role in tumorigenesis. Tumor cells present distinct phenotypic states with different functional attributes, among which, cancer stem cells (CSCs) possess the principal properties of self-renewal, clonal tumor initiation capacity, and clonal long-term repopulation potential (Plaks et al., 2015). From this perspective, CSCs can facilitate the initiation and progression of tumors, which may induce drug resistance (Plaks et al., 2015; Malta et al., 2018). Thus, we evaluated the correlation between the prognostic signature and the tumor stemness of OSCC. Based on DNAss and RNAss, patients with OSCC in the HRisk group presented higher tumor stemness than those in the LRisk group, and the risk score showed a positive connection with tumor stemness. These results suggested that this oxidative stress-related signature could predict tumor progression and invasion for patients with OSCC and simultaneously explain the poor prognosis of high-risk patients.

Presently, personalized medicine, concentrating on designing specific therapeutics that are best suited for an individual patient based on genome information, is considered to be the future direction in oncotherapy (Jackson and Chester, 2015). To verify whether the prognostic signature could guide clinicians to develop an effective personalized treatment decision for patients with OSCC, we thoroughly investigated the correlation between the constructed signature and response to immunotherapy or chemotherapy. Immunotherapy, particularly ICI treatment, offers a reliable alternative or adjunctive therapy to conventional therapies for refractory patients with OSCC (Dorta-Estremera et al., 2019). The TIDE score has been identified as a potent biomarker to predict the tumor response to anti-PD1 or anti-CTLA4 (Jiang et al., 2018), while the TMB score has also been validated as an available tool to help oncologists select patients who may benefit from ICIs (Choucair et al., 2020). Encouragingly, the signature in our study was not only closely correlated with both TIDE and TMB scores but was also efficient in predicting the immunotherapy response with high accuracy (AUC = 0.766). Next, we comprehensively evaluated the efficacy of chemotherapy drugs in different risk subgroups and their relationship to the risk score on the basis of three public drug sensitivity databases (CGP, GDSC, and CTRP) (Garnett et al., 2012; Basu et al., 2013; Yang et al., 2013). The IC50 values of eight common chemotherapy drugs for OSCC (Cisplatin, Paclitaxel, Cytarabine, Docetaxel, Doxorubicin, Gemcitabine, Methotrexate, and 5-Fluorouracil), on the basis of the NCCN guidelines Version 2.2021, and one novel anti-cancer drug, Gefitinib, based on previous studies, were calculated and estimated. Consequently, we found that patients with low-risk scores were more sensitive to 5-Fluorouracil, Gemcitabine, and Gefitinib, while those with high-risk scores showed more sensitive responses to paclitaxel and docetaxel, which indicated that the prognostic signature could facilitate personalized chemotherapy decisions. Furthermore, 5-Fluorouracil, Gemcitabine, and Gefitinib generate ROS, which can cause cancer cell death by inducing oxidative damage, whereas CSCs can increase cellular adaptive responses to ROS to result in chemoresistance (Okon et al., 2015; Blondy et al., 2020; Xue et al., 2020). As mentioned above, OSCC patients with high-risk scores presented notably elevated tumor stemness, which may explain their higher resistance to 5-Fluorouracil, Gemcitabine, and Gefitinib. Moreover, paclitaxel and docetaxel are both anti-cancer agents belonging to the taxane family and can inhibit cancer cell proliferation by inducing cell cycle arrest (Ashrafizadeh et al., 2021). On the basis of the GSEA above, cell cycle-related pathways, such as cell cycle, spliceosome, and base excision repair, were apparently activated with increased risk scores, which shed light on the higher sensitivity of high-risk patients to Paclitaxel and Docetaxel. To summarize, considering the strong correlation of the prognostic signature with both immunotherapy and chemotherapy, we present a valid and robust tool that can guide clinicians to make effective personalized treatment decisions for patients with OSCC with different risk levels.

Despite considerable strides in treatment regimens and new drugs for OSCC, the survival rate has been poor and unsatisfactory in recent decades (Bray et al., 2018; Nör and Gutkind, 2018). Thus, the development of novel drugs for OSCC is still necessary. The new use of old drugs is more cost-effective compared to the exploration and development of novel drugs (Cheng et al., 2021). CMap, a public database for predicting small molecule drugs according to transcriptional expression data, has been used to predict and screen potential novel drug candidates for cancers in previous studies (Cheng et al., 2021; Wang M. et al., 2021). Therefore, to uncover novel drug candidates for patients with OSCC with high-risk scores, we evaluated DEGs between the two risk subgroups and finally screened out 10 potential therapeutic agents, including Dextromethorphan and AH-6809, in the CMap database, all of which have been barely explored in OSCC. Dextromethorphan is a safe FDA-approved drug with few undesirable side effects and usually serves as an effective antitussive agent (Silva and Dinis-Oliveira, 2020). Wang et al. found that Dextromethorphan and Metformin, at their pharmacological doses, could synergistically repress nicotine-enhanced cancer-initiating cell properties and halt tumor progression by directly targeting CHRNA7 to inhibit JAK2/STAT3/SOX2 signaling in esophageal squamous cell carcinoma and perhaps other nicotine-sensitive cancer types (Wang et al., 2021a). Additionally, AH-6809 was reported to inhibit the proliferation of non-small cell lung cancer cells by antagonizing prostaglandin receptors (Casibang and Moody, 2002). Hence, the results of CMap analysis may provide some new promising drug candidates for high-risk patients with OSCC and identify a viable direction for the future research on chemotherapy drugs.

Conclusion

In summary, we established a novel prognostic signature with six OSGs. On the basis of both TCGA-OSCC and GSE41613 cohorts, the signature was proven to be an independent prognostic factor with high accuracy, as well as an impactful indicator for predicting the prognosis and immune status of patients with OSCC. Meanwhile, the constructed signature demonstrated that patients with high-risk scores might benefit more from ICI treatment compared to those with low-risk scores, and the risk score presented a close interaction with the TME and chemotherapy sensitivity. Our findings may also provide valuable new insight into the roles of oxidative stress in the prognosis, TME, immune status, immunotherapy response, and chemotherapy sensitivity of patients with OSCC. The oxidative stress-related signature may provide a valid and robust tool that can not only efficiently predict the prognosis and immune status but also guide clinicians to develop effective personalized therapeutic strategies for patients with OSCC.

Data availability statement

All datasets presented in this study are included in the article/Supplementary Material.

Author contributions

WL and SH selected the subject and designed the research approach; WL and YW constructed the risk model via bioinformatic analysis; YW performed the statistical analysis; TZ drew the figures and made the tables; WL and YW wrote sections of the manuscript. All authors approved the final version of the manuscript.

Funding

This work was supported by funding from the Shandong Provincial Natural Science Foundation (ZR2021MH270) and Jinan Clinical Medical Science and Technology Innovation Plan (202134035).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.977902/full#supplementary-material

Supplementary Figure S1 | A forest plot of the univariate Cox regression analysis with the prognostic ORGs (A) and the 1-, 3- and 5-year ROC curves and AUC values of the signature in TCGA cohort (B) and GEO cohort (C). OSGs, Oxidative stress-related genes, ROC, Receiver operating characteristic, AUC, Area under curve, TCGA, The Cancer Genome Atlas, GEO, The Gene Expression Omnibus.

Supplementary Figure S2 | Expression levels of model genes between low- and high-risk subgroups in TCGA (A) and GEO (B) cohorts. *** p < 0.001. GEO, The Gene Expression Omnibus, TCGA, The Cancer Genome Atlas.

Supplementary Figure S3 | Survival analyses of model genes in OSCC using Kaplan-Meier algorithm based on TCGA (A-F) and GEO (G-L) cohorts. OSCC, Oral squamous cell carcinoma, GEO, The Gene Expression Omnibus, TCGA, The Cancer Genome Atlas.

Supplementary Figure S4 | Survival analyses of clinical features in OSCC using Kaplan-Meier algorithm based on TCGA and GEO cohorts. TCGA cohort: (A) age, (B) gender, (C) grade, (D) clinical stage, (E) T stage and (F) N stage. GEO cohort: (G) age, (H) gender, and (I) clinical stage. OSCC, Oral squamous cell carcinoma, GEO, The Gene Expression Omnibus, TCGA, The Cancer Genome Atlas.

Supplementary Figure S5 | Correlations between the risk score and the levels of immune score, stromal score and tumor purity according to TCGA (A-C) and GEO (D-E) cohorts. GEO, The Gene Expression Omnibus, TCGA, The Cancer Genome Atlas.

Supplementary Figure S6 | Immune cells and immune-related functions with prognostic values in OSCC on basis of TCGA and GEO cohorts. OSCC, Oral squamous cell carcinoma, GEO, The Gene Expression Omnibus, TCGA, The Cancer Genome Atlas.

Supplementary Figure S7 | The IC50 of selected chemotherapeutic agents between different risk subgroups of OSCC and correlations between the risk score and the IC50 of selected chemotherapeutic agents in TCGA cohort based on CGP, GDSC2 and CTRP databases. OSCC, Oral squamous cell carcinoma, IC50, Half-maximal inhibitory concentration, TCGA, The Cancer Genome Atlas, CGP, the Cancer Genome Project, GDSC, Genomics of Drug Sensitivity in Cancer, CTRP, Cancer Therapeutics Response Portal.

Supplementary Figure S8 | The 2D chemical structures of the selected small molecule candidates for treating OSCC patients using CMap algorithm online. OSCC, Oral squamous cell carcinoma, CMap, the Connectivity Map.

References

Ahmad, I., Fakhri, S., Khan, H., Jeandet, P., Aschner, M., and Yu, Z. L. (2020). Targeting cell cycle by β-carboline alkaloids in vitro: Novel therapeutic prospects for the treatment of cancer. Chem. Biol. Interact. 330, 109229. doi:10.1016/j.cbi.2020.109229

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, P., Aptsiauri, N., Ruiz-Cabello, F., and Garrido, F. (2021). HLA class I loss in colorectal cancer: Implications for immune escape and immunotherapy. Cell. Mol. Immunol. 18, 556–565. doi:10.1038/s41423-021-00634-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashrafizadeh, M., Mirzaei, S., Hashemi, F., Zarrabi, A., Zabolian, A., Saleki, H., et al. (2021). New insight towards development of paclitaxel and docetaxel resistance in cancer cells: EMT as a novel molecular mechanism and therapeutic possibilities. Biomed. Pharmacother. 141, 111824. doi:10.1016/j.biopha.2021.111824

PubMed Abstract | CrossRef Full Text | Google Scholar

Bagheri, S., Saboury, A. A., and Haertlé, T. (2019). Adenosine deaminase inhibition. Int. J. Biol. Macromol. 141, 1246–1257. doi:10.1016/j.ijbiomac.2019.09.078

PubMed Abstract | CrossRef Full Text | Google Scholar

Basu, A., Bodycombe, N. E., Cheah, J. H., Price, E. V., Liu, K., Schaefer, G. I., et al. (2013). An interactive resource to identify cancer genetic and lineage dependencies targeted by small molecules. Cell 154, 1151–1161. doi:10.1016/j.cell.2013.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Biselli-Chicote, P. M., Lotierzo, A. T., Biselli, J. M., Paravino, É. C., and Goloni-Bertollo, E. M. (2019). Atorvastatin increases oxidative stress and inhibits cell migration of oral squamous cell carcinoma in vitro. Oral Oncol. 90, 109–114. doi:10.1016/j.oraloncology.2019.01.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Blondy, S., David, V., Verdier, M., Mathonnet, M., Perraud, A., and Christou, N. (2020). 5-Fluorouracil resistance mechanisms in colorectal cancer: From classical pathways to promising processes. Cancer Sci. 111, 3142–3154. doi:10.1111/cas.14532

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Ca. Cancer J. Clin. 68, 394–424. doi:10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

Bron, L., Jandus, C., Andrejevic-Blant, S., Speiser, D. E., Monnier, P., Romero, P., et al. (2013). Prognostic value of arginase-II expression and regulatory T-cell infiltration in head and neck squamous cell carcinoma. Int. J. Cancer 132, E85–E93. doi:10.1002/ijc.27728

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, T., Boland, A., Bagust, A., Oyee, J., Hockenhull, J., Dundar, Y., et al. (2010). Gefitinib for the first-line treatment of locally advanced or metastatic non-small cell lung cancer. Health Technol. Assess. 14, 71–79. doi:10.3310/hta14suppl2/10

PubMed Abstract | CrossRef Full Text | Google Scholar

Bundela, S., Sharma, A., and Bisen, P. S. (2014). Potential therapeutic targets for oral cancer: ADM, TP53, EGFR, lyn, CTLA4, skil, ctgf, CD70. PloS One 9, e102610. doi:10.1371/journal.pone.0102610

PubMed Abstract | CrossRef Full Text | Google Scholar

Casibang, M., and Moody, T. W. (2002). AH6809 antagonizes non-small cell lung cancer prostaglandin receptors. Lung cancer 36, 33–42. doi:10.1016/s0169-5002(01)00476-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Y., Hou, K., Wang, Y., Chen, Y., Zheng, X., Qi, J., et al. (2021). Identification of prognostic signature and gliclazide as candidate drugs in lung adenocarcinoma. Front. Oncol. 11, 665276. doi:10.3389/fonc.2021.665276

PubMed Abstract | CrossRef Full Text | Google Scholar

Choucair, K., Morand, S., Stanbery, L., Edelman, G., Dworkin, L., and Nemunaitis, J. (2020). Tmb: A promising immune-response biomarker, and potential spearhead in advancing targeted therapy trials. Cancer Gene Ther. 27, 841–853. doi:10.1038/s41417-020-0174-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Choudhari, S. K., Chaudhary, M., Gadbail, A. R., Sharma, A., and Tekade, S. (2014). Oxidative and antioxidative mechanisms in oral cancer and precancer: A review. Oral Oncol. 50, 10–18. doi:10.1016/j.oraloncology.2013.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cubillos-Ruiz, J. R., Silberman, P. C., Rutkowski, M. R., Chopra, S., Perales-Puchalt, A., Song, M., et al. (2015). ER stress sensor XBP1 controls anti-tumor immunity by disrupting dendritic cell homeostasis. Cell 161, 1527–1538. doi:10.1016/j.cell.2015.05.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, S., Kar Mahapatra, S., Gautam, N., Das, A., and Roy, S. (2007). Oxidative stress in lymphocytes, neutrophils, and serum of oral cavity cancer patients: Modulatory array of L-glutamine. Support. Care Cancer 15, 1399–1405. doi:10.1007/s00520-007-0266-3

PubMed Abstract | CrossRef Full Text | Google Scholar

de Matos, F. R., Santos, E. M., Santos, H. B. P., Machado, R. A., Lemos, T., Coletta, R. D., et al. (2019). Association of polymorphisms in IL-8, MMP-1 and MMP-13 with the risk and prognosis of oral and oropharyngeal squamous cell carcinoma. Arch. Oral Biol. 108, 104547. doi:10.1016/j.archoralbio.2019.104547

PubMed Abstract | CrossRef Full Text | Google Scholar

de Sá Junior, P. L., Câmara, D. A. D., Porcacchia, A. S., Fonseca, P. M. M., Jorge, S. D., Araldi, R. P., et al. (2017). The roles of ROS in cancer heterogeneity and therapy. Oxid. Med. Cell. Longev. 2017, 2467940. doi:10.1155/2017/2467940

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorta-Estremera, S., Hegde, V. L., Slay, R. B., Sun, R., Yanamandra, A. V., Nicholas, C., et al. (2019). Targeting interferon signaling and CTLA-4 enhance the therapeutic efficacy of anti-PD-1 immunotherapy in preclinical model of HPV(+) oral cancer. J. Immunother. Cancer 7, 252. doi:10.1186/s40425-019-0728-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunn, G. P., Koebel, C. M., and Schreiber, R. D. (2006). Interferons, immunity and cancer immunoediting. Nat. Rev. Immunol. 6, 836–848. doi:10.1038/nri1961

PubMed Abstract | CrossRef Full Text | Google Scholar

Garnett, M. J., Edelman, E. J., Heidorn, S. J., Greenman, C. D., Dastur, A., Lau, K. W., et al. (2012). Systematic identification of genomic markers of drug sensitivity in cancer cells. Nature 483, 570–575. doi:10.1038/nature11005

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One 9, e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorrini, C., Harris, I. S., and Mak, T. W. (2013). Modulation of oxidative stress as an anticancer strategy. Nat. Rev. Drug Discov. 12, 931–947. doi:10.1038/nrd4002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, X., Boldrup, L., Coates, P. J., Fahraeus, R., Wang, L., Wilms, T., et al. (2019). High immune cytolytic activity in tumor-free tongue tissue confers better prognosis in patients with squamous cell carcinoma of the oral tongue. J. Pathol. Clin. Res. 5, 240–247. doi:10.1002/cjp2.138

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadler-Olsen, E., and Wirsing, A. M. (2019). Tissue-infiltrating immune cells as prognostic markers in oral squamous cell carcinoma: A systematic review and meta-analysis. Br. J. Cancer 120, 714–727. doi:10.1038/s41416-019-0409-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanakawa, H., Orita, Y., Sato, Y., Takeuchi, M., Ohno, K., Gion, Y., et al. (2014). Regulatory T-cell infiltration in tongue squamous cell carcinoma. Acta Otolaryngol. 134, 859–864. doi:10.3109/00016489.2014.918279

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayes, J. D., Dinkova-Kostova, A. T., and Tew, K. D. (2020). Oxidative stress in cancer. Cancer Cell 38, 167–197. doi:10.1016/j.ccell.2020.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinshaw, D. C., and Shevde, L. A. (2019). The tumor microenvironment innately modulates cancer progression. Cancer Res. 79, 4557–4566. doi:10.1158/0008-5472.Can-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, L., Yu, X., Jiang, Z., and Zeng, P. (2021a). Novel autophagy-related gene signature investigation for patients with oral squamous cell carcinoma. Front. Genet. 12, 673319. doi:10.3389/fgene.2021.673319

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Z., Lan, T., Wang, J., Chen, Z., and Zhang, X. (2021b). Identification and validation of seven RNA binding protein genes as a prognostic signature in oral cavity squamous cell carcinoma. Bioengineered 12, 7248–7262. doi:10.1080/21655979.2021.1974328

PubMed Abstract | CrossRef Full Text | Google Scholar

Hundsdorfer, B., Zeilhofer, H. F., Bock, K. P., Dettmar, P., Schmitt, M., Kolk, A., et al. (2005). Tumour-associated urokinase-type plasminogen activator (uPA) and its inhibitor PAI-1 in normal and neoplastic tissues of patients with squamous cell cancer of the oral cavity - clinical relevance and prognostic value. J. Craniomaxillofac. Surg. 33, 191–196. doi:10.1016/j.jcms.2004.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Jackson, S. E., and Chester, J. D. (2015). Personalised cancer medicine. Int. J. Cancer 137, 262–266. doi:10.1002/ijc.28940

PubMed Abstract | CrossRef Full Text | Google Scholar

Jewett, A., Kos, J., Fong, Y., Ko, M. W., Safaei, T., Perišić Nanut, M., et al. (2018). NK cells shape pancreatic and oral tumor microenvironments; role in inhibition of tumor growth and metastasis. Semin. Cancer Biol. 53, 178–188. doi:10.1016/j.semcancer.2018.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Koike, K., Dehari, H., Shimizu, S., Nishiyama, K., Sonoda, T., Ogi, K., et al. (2020). Prognostic value of HLA class I expression in patients with oral squamous cell carcinoma. Cancer Sci. 111, 1491–1499. doi:10.1111/cas.14388

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutryb-Zajac, B., Koszalka, P., Mierzejewska, P., Bulinska, A., Zabielska, M. A., Brodzik, K., et al. (2018). Adenosine deaminase inhibition suppresses progression of 4T1 murine breast cancer by adenosine receptor-dependent mechanisms. J. Cell. Mol. Med. 22, 5939–5954. doi:10.1111/jcmm.13864

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeser, D., Gruener, R. F., and Huang, R. S. (2021). oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief. Bioinform. 22, bbab260. doi:10.1093/bib/bbab260

PubMed Abstract | CrossRef Full Text | Google Scholar

Magnussen, S. N., Hadler-Olsen, E., Costea, D. E., Berg, E., Jacobsen, C. C., Mortensen, B., et al. (2017). Cleavage of the urokinase receptor (uPAR) on oral cancer cells: Regulation by transforming growth factor - β1 (TGF-β1) and potential effects on migration and invasion. BMC cancer 17, 350. doi:10.1186/s12885-017-3349-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Maj, T., Wang, W., Crespo, J., Zhang, H., Wang, W., Wei, S., et al. (2017). Oxidative stress controls regulatory T cell apoptosis and suppressor activity and PD-L1-blockade resistance in tumor. Nat. Immunol. 18, 1332–1341. doi:10.1038/ni.3868

PubMed Abstract | CrossRef Full Text | Google Scholar

Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., et al. (2018). Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell 173, 338–354. e315. doi:10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantovani, G., Macciò, A., Madeddu, C., Mura, L., Gramignano, G., Lusso, M. R., et al. (2002). Quantitative evaluation of oxidative stress, chronic inflammatory indices and leptin in cancer patients: Correlation with stage and performance status. Int. J. Cancer 98, 84–91. doi:10.1002/ijc.10143

PubMed Abstract | CrossRef Full Text | Google Scholar

Nör, J. E., and Gutkind, J. S. (2018). Head and neck cancer in the new era of precision medicine. J. Dent. Res. 97, 601–602. doi:10.1177/0022034518772278

PubMed Abstract | CrossRef Full Text | Google Scholar

Okon, I. S., Coughlan, K. A., Zhang, M., Wang, Q., and Zou, M. H. (2015). Gefitinib-mediated reactive oxygen specie (ROS) instigates mitochondrial dysfunction and drug resistance in lung cancer cells. J. Biol. Chem. 290, 9101–9110. doi:10.1074/jbc.M114.631580

PubMed Abstract | CrossRef Full Text | Google Scholar

Plaks, V., Kong, N., and Werb, Z. (2015). The cancer stem cell niche: How essential is the niche in regulating stemness of tumor cells? Cell stem Cell 16, 225–238. doi:10.1016/j.stem.2015.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, X., Hou, Q. H., Shi, Q. Y., Jiang, H. X., and Qin, S. Y. (2020). Identification of hub prognosis-associated oxidative stress genes in pancreatic cancer using integrated bioinformatics analysis. Front. Genet. 11, 595361. doi:10.3389/fgene.2020.595361

PubMed Abstract | CrossRef Full Text | Google Scholar

Reyimu, A., Chen, Y., Song, X., Zhou, W., Dai, J., and Jiang, F. (2021). Identification of latent biomarkers in connection with progression and prognosis in oral cancer by comprehensive bioinformatics analysis. World J. Surg. Oncol. 19, 240. doi:10.1186/s12957-021-02360-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G., and Hacohen, N. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 160, 48–61. doi:10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Sales, C. B., Buim, M. E., de Souza, R. O., de Faro Valverde, L., Mathias Machado, M. C., Reis, M. G., et al. (2016). Elevated VEGFA mRNA levels in oral squamous cell carcinomas and tumor margins: A preliminary study. J. Oral Pathol. Med. 45, 481–485. doi:10.1111/jop.12398

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, A. R., and Dinis-Oliveira, R. J. (2020). Pharmacokinetics and pharmacodynamics of dextromethorphan: Clinical and forensic aspects. Drug Metab. Rev. 52, 258–282. doi:10.1080/03602532.2020.1758712

PubMed Abstract | CrossRef Full Text | Google Scholar

Sosa, V., Moliné, T., Somoza, R., Paciucci, R., Kondoh, H., and Me, L. L. (2013). Oxidative stress and cancer: An overview. Ageing Res. Rev. 12, 376–390. doi:10.1016/j.arr.2012.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Stich, H. F., and Anders, F. (1989). The involvement of reactive oxygen species in oral cancers of betel quid/tobacco chewers. Mutat. Res. 214, 47–61. doi:10.1016/0027-5107(89)90197-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, X., He, J., Li, B., Zheng, Y., Li, K., Zou, S., et al. (2019). Efficacy and safety of Gefitinib in patients with advanced head and neck squamous cell carcinoma: A meta-analysis of randomized controlled trials. J. Oncol. 2019, 6273438. doi:10.1155/2019/6273438

PubMed Abstract | CrossRef Full Text | Google Scholar

Toyokuni, S., Okamoto, K., Yodoi, J., and Hiai, H. (1995). Persistent oxidative stress in cancer. FEBS Lett. 358, 1–3. doi:10.1016/0014-5793(94)01368-b

PubMed Abstract | CrossRef Full Text | Google Scholar

Verbon, E. H., Post, J. A., and Boonstra, J. (2012). The influence of reactive oxygen species on cell cycle progression in mammalian cells. Gene 511, 1–6. doi:10.1016/j.gene.2012.08.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Verneau, J., Sautés-Fridman, C., and Sun, C. M. (2020). Dendritic cells in the tumor microenvironment: Prognostic and theranostic impact. Semin. Immunol. 48, 101410. doi:10.1016/j.smim.2020.101410

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Du, L., Xiong, X., Lin, Y., Zhu, J., Yao, Z., et al. (2021a). Repurposing dextromethorphan and metformin for treating nicotine-induced cancer by directly targeting CHRNA7 to inhibit JAK2/STAT3/SOX2 signaling. Oncogene 40, 1974–1987. doi:10.1038/s41388-021-01682-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Wang, Y., Han, N., Wang, X., and Ruan, M. (2021b). HPRT promotes proliferation and metastasis in head and neck squamous cell carcinoma through direct interaction with STAT3. Exp. Cell Res. 399, 112424. doi:10.1016/j.yexcr.2020.112424

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Zhong, B., Li, M., Wang, Y., Yang, H., and Du, K. (2021c). Identification of potential core genes and pathways predicting pathogenesis in head and neck squamous cell carcinoma. Biosci. Rep. 41, BSR20204148. doi:10.1042/bsr20204148

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Wang, L., Wen, Z., and Yao, J. (2021). Integrated analysis identifies oxidative stress genes associated with progression and prognosis in gastric cancer. Sci. Rep. 11, 3292. doi:10.1038/s41598-021-82976-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, Y., Li, H., Mao, L., Yang, Q. C., Fu, L. Q., Wu, C. C., et al. (2019). CD103(+) T and dendritic cells indicate a favorable prognosis in oral cancer. J. Dent. Res. 98, 1480–1487. doi:10.1177/0022034519882618

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, D., Zhou, X., and Qiu, J. (2020). Emerging role of NRF2 in ROS-mediated tumor chemoresistance. Biomed. Pharmacother. 131, 110676. doi:10.1016/j.biopha.2020.110676

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of drug sensitivity in cancer (GDSC): A resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Z., Zhou, J., Xiong, X., Liu, D., Zheng, Y., Ding, X., et al. (2019). Prognostic value of serum liver enzymes in oral and oropharynx squamous cell carcinomas. J. Oral Pathol. Med. 48, 36–42. doi:10.1111/jop.12803

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, C., Huang, D., Wei, H., Wu, S., Song, J., Tang, Z., et al. (2021). Identification of immune-related risk signatures for the prognostic prediction in oral squamous cell carcinoma. J. Immunol. Res. 2021, 6203759. doi:10.1155/2021/6203759

PubMed Abstract | CrossRef Full Text | Google Scholar