- Open Access
Prognostic biomarkers and therapeutic targets in oral squamous cell carcinoma: a study based on cross-database analysis
Hereditas volume 158, Article number: 15 (2021)
Oral squamous cell carcinoma (OSCC) is a malignant cancer, the survival rate of patients is disappointing. Therefore, it is necessary to identify the driven-genes and prognostic biomarkers in OSCC.
Four Gene Expression Omnibus (GEO) datasets were integratedly analyzed using bioinformatics approaches, including identification of differentially expressed genes (DEGs), GO and KEGG analysis, construction of protein-protein interaction (PPI) network, selection of hub genes, analysis of prognostic information and genetic alterations of hub genes. ONCOMINE, The Cancer Genome Atlas (TCGA) and Human Protein Atlas databases were used to evaluate the expression and prognostic value of hub genes. Tumor immunity was assessed to investigate the functions of hub genes. Finally, Cox regression model was performed to construct a multiple-gene prognostic signature.
Totally 261 genes were found to be dysregulated. 10 genes were considered to be the hub genes. The Kaplan-Meier analysis showed that upregulated SPP1, FN1, CXCL8, BIRC5, PLAUR, and AURKA were related to poor outcomes in OSCC patients. FOXM1 and TPX2 were considered as the potential immunotherapeutic targets with future clinical significance. Moreover, we constructed a nine-gene signature (TEX101, DSG2, SCG5, ADA, BOC, SCARA5, FST, SOCS1, and STC2), which can be utilized to predict prognosis of OSCC patients effectively.
These findings may provide new clues for exploring the molecular mechanisms and targeted therapy in OSCC. The hub genes and risk gene signature are helpful to the personalized treatment and prognostic judgement.
Oral squamous cell carcinoma (OSCC) is the major type of head and neck squamous cell carcinoma (HNSCC) . Despite great works have been made on early screen and personalized treatment for cancers, OSCC still is a challengeable disease and has brought seriously economic and medical burden . Risk factors, including smoking, drinking, and HPV infections, are closely associated with the development of OSCC [3,4,5]. Nevertheless, the mechanisms of OSCC are still unclear. Moreover, most OSCC patients can’t be screened early due to the lack of available diagnostic markers. In addition, due to the drug resistance, some patients with OSCC might suffer from cancer recurrence. Thus, identifying the novel biomarkers and effective targets is of great importance to OSCC research and management.
Recently, the gene chip assays and second-generation gene sequencing have been extensively applied in scientific researches [6, 7]. These approaches could effectively screen the key genes that influence the cancer development or progression . To date, numerous studies have used the second-generation gene sequencing or gene chip assays to explore the key genes in OSCC . However, the findings might be inconsistent due to the tumor heterogeneity, localization, HPV-related link, and microbiota. So far, there are few reliable markers and therapeutic targets for OSCC. The integrated bioinformatics analysis may solve these problems and eventually find more convincing results since it uses several bioinformatics methods and integrates the data from different gene profiles .
In this study, we utilized microarray data of OSCC tissues and normal oral tissues in GEO and TCGA databases to screen the hub genes and biomarkers. We purposed to screen the hub genes, important GO terms, and significant pathways in the OSCC progression, thus helping to reveal the mechanisms of this disease. We also hope to pick out the risk genes and construct a multiple-gene prognostic signature, which can be implemented for prognostic judgment in OSCC.
Material and methods
Data collection and procession
To acquire the OSCC mRNA expression datasets, the keywords: “Expression profiling by array” “Homo sapiens”, and “Oral squamous cell carcinoma”, were searched in GEO database (http://www.ncbi.nlm.nih.gov/geo/) [11, 12]. After a systematic review, OSCC and non-tumor oral tissues gene expression profiles of GSE23558 , GSE30784 , GSE74530 , and GSE37991  were selected and downloaded through getGEO function in “GEOquery” package . The criteria for selecting the datasets as following: (1) dataset contains both OSCC tissues and non-tumor tissues; (2) samples size was 10 or more. The data were pre-processed as the previous studies [18, 19]. Briefly, the raw data of gene chips were normalized by “limma” package in R software (Version: 3.6.3). Moreover, the “sva” package was utilized to remove the batch effect. The detail information of the four GEO datasets was presented in Table 1. The RNA-seq data and clinicopathological data of 502 OSCC patients and 44 normal samples were also downloaded from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/) database. Data procession of TCGA database was performed as the previous studies [20, 21].
Identification of DEGs
The limma package was utilized to screen the DEGs between OSCC and normal tissues in GEO datasets and TCGA dataset . For selecting the DEGs in each GEO dataset, |logFC| >1 and adjust P-value <0.05 were set as the cut-off criteria . Then, the overlapping DEGs was screened using a Venn tool (http://bioinfogp.cnb.csic.es/tools/venny/). The data profile of GSE37991 was used as a reference to construct the heatmap and to present the distribution of DEGs.
Functional enrichment analysis for DEGs
GO and KEGG analysis of overlapping DEGs were conducted via the DAVID database (https://david.ncifcrf.gov/), and FunRich tool (FunRich 3.0) as descried previously [12, 23,24,25]. We submitted the overlapping DEGs into the above databases or software. The top five GO terms for biological process (BP), cellular component (CC), and molecular function (MF) were illustrated as bar charts . The KEGG results were visualized by “clusterProfiler” package, the top 10 KEGG pathways of upregulated DEGs and downregulated DEGs were showed as bubble charts, respectively .
Construction of protein-protein interaction (PPI) network
PPI networks are formed by proteins due to the existence of biochemical or electrostatic forces . Here, the Search Tool for the Retrieval of Interacting Genes (STRING 11.0) database (https://string-db.org/cgi/input.pl) was applied to establish PPI networks [12, 27]. Cytoscape software (v3.6.1) was used to illustrate the PPI networks and the cut-off criteria were set as confidence score ≥ 0.4, maximum number of interactors = 0 [12, 28]. Gene clusters visualization are conducted as our previous studies [12, 28]. Gene clusters were screened with the following criteria: MCODE scores >10 and number of nodes >10 . 10 genes were considered to be the hub genes according to connectivity degree [12, 28, 29].
Hub genes validation
As for hub genes validation, we assessed the levels both in mRNA level and protein level. Specifically, ONCOMINE (https://www.oncomine.org) database was utilized to evaluate the mRNA expression of selected hub genes . The gene rank means that the median rank of one searched gene across the selected analyses . We also used GEPIA database(http://gepia.cancer-pku.cn/index.html) to validate the mRNA levels of hub genes [12, 32]. The Human Protein Atlas database collected more than 11,200 unique proteins , we thus used it to evaluate the protein levels of hub genes in OSCC tissues and normal control tissues.
Genetic alterations and survival analysis
The cBio Cancer Genomics Portal (http://www.cbioportal.org/) is an online database which enables us to compare the genetic alterations of the hub genes in OSCC . Subsequently, survival analysis was performed in the Kaplan Meier-plotter website (http://kmplot.com/analysis/index.php), which could assess the effect of 54,675 genes on survival using 18,674 cancer samples from GEO, European Genome-phenome Archive (EGA) and TCGA database .
Analysis of potential target genes
TIMER website (https://cistrome.shinyapps.io/timer/) was applied to investigate the expression levels of the selected target genes in different human tumors . The UALCAN website (http://ualcan.path.uab.edu/analysis.html) was utilized to assess the factors which are associated with the expression levels of the target genes in OSCC .
Correlation analysis of genes expression and immune cell infiltration and immune checkpoints
TIMER website was applied to analyze the genes expression data for FOXM1 and TPX2 in TCGA OSCC samples and its correlation with tumor infiltration of six immune cell types (B cells, CD4+T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) and five immunological checkpoints (CD274, CTLA4, PDCD1, PDCD1LG2, and TOX) .
Construction of prognostic gene signature
Univariate and multivariate Cox regression analyses were conducted to investigate the relationships between the expression of 261 overlapping DEGs and OSCC patients’ survival in TCGA dataset. Nine DEGs were identified as the prognostic indicators in OSCC. These nine DEGs were utilized to construct a prognostic signature as previous study . Each OSCC patient received a risk score according to the formula: risk score = (Coefficient gene 1 × Expression gene 1) + (Coefficient gene 2 × Expression gene 2) +…+ (Expression gene n × Coefficient gene n). Finally, the Kaplan-Meier curve and receiver operating characteristic (ROC) curve were used to evaluate the efficiency of the gene signature.
Identification of the independent prognostic indicators in OSCC
Univariate and multivariate Cox regression analyses were conducted to identify the independent prognostic indicators (including age, gender, grade, stage, T, M, N, and risk score) for OSCC patients.
The R statistical package (R version 3.6.3) and Perl language and were utilized to conduct the statistical tests and graphics unless otherwise stated. P <0.05 was regarded as statistical significance.
DEGs involved in OSCC
The four datasets (GSE23558, GSE30784, GSE74530, and GSE37991), including 240 OSCC tissues and 95 non-tumor tissues, were included in this study. We extracted 3383, 2532, 2106, and 2527 DEGs from GSE23588, GSE30784, GSE37991, and GSE74530, respectively (Fig. 1a-d). Totally 261 overlapping DEGs were screened from the 4 datasets (Fig. 1e; Table 2), including 135 upregulated DEGs and 126 downregulated DEGs. The data profile of GSE37991 was used as a reference to construct the heatmap and to show the differential distribution of DEGs (Fig. 1f).
Functional enrichment analysis
Functional enrichment analysis was analyzed by DAVID database and FunRich tool (FunRich 3.0). As shown in Figure S1a, in BP group, upregulated DEGs were mainly involved in energy pathways, metabolism, and lipid storage, whereas the downregulated DEGs were related to signal transduction, cell communication, and immune response (Figure S1b). For CC group, upregulated DEGs were mainly associated with exosomes, extracellular space, and endoplasmic reticulum membrane, whereas downregulated DEGs were involved in extracellular, extracellular region, and extracellular space. As for MF, upregulated DEGs were correlated to fucosyltransferase activity, catalytic activity, and cytokine activity. The downregulated DEGs were significantly connected with to extracellular matrix structural constituent, metallopeptidase activity, and cytokine activity. According to the GO results, subsequent studies on OSCC should focus on how to balance the cellular microenvironment or reshape the cellular physiological functions of cancer cells.
KEGG enrichment showed that upregulated DEGs were mainly involved in pathway in cancer, PI3K-Akt pathway, ECM receptor interaction, and focal adhesion (Figure S1c). Previous studies have showed that molecular pathways involved in OSCC are complex . The PI3K-Akt and Wnt/β-catenin signaling were demonstrated to be the three major interlinked pathways involved in the molecular pathogenesis of OSCC . Here, we found that pathway in cancer and PI3K-Akt pathway were associated with the upregulated DEGs. Therefore, monitoring these signaling pathways may aid to decide appropriate therapeutic approaches in OSCC patients. The downregulated DEGs were enriched metabolic pathway, serotonergic synapse, tyrosine metabolism, and arachidonic acid metabolism (Figure S1d). Studies have showed that metabolic alterations might provide energy and nutrients for sustaining the cancer proliferation and growth . In this study, downregulated DEGs were closely related to metabolic pathway, tyrosine metabolism, and arachidonic acid metabolism. Therefore, these metabolism-targeted pathways may help us to improve the treatments efficiency.
PPI network analysis
Totally 261 overlapping DEGs were mapped into PPI network via STRING database and Cytoscape software, (Fig. 2a). The top two important clusters were picked out by MCODE plug-in in Cytoscape software (Cluster 1, MCODE score =15.2; Cluster 2, MCODE score = 9.789) (Fig. 2b-c). Cluster 1 consists of 16 nodes and 144 edges (Table S1), which are mainly related to cell cycle, DNA replication, and cellular senescence. Cluster 2 consists of 20 nodes and 93 edges (Table S2), which are mainly enriched in IL-17 signaling pathway, influenza A, and complement and coagulation cascades. Using cytoHubba software, 10 genes (FN1, CXCL8, CXCL10, SPP1, FOXM1, AURKA, ISG15, PLAUR, TPX2, and BIRC5) were selected as hub genes according to the connectivity degree (Table 3). According to Figure S2a, the hub genes could interact with each other and they might be the driven-genes in OSCC development and progression. KEGG pathways analysis showed that the significantly enriched terms for the hub genes were RIG-I-like receptor signaling pathway, Toll-like receptor signaling pathway, IL-17 signaling pathway, Influenza A, Cellular senescence, Chemokine signaling pathway, and Cytokine-cytokine receptor interaction (Figure S2b).
Validation the expression of hub genes
ONCOMINE, TCGA, and The Human Protein Atlas databases were used to validate the expression of hub genes. Firstly, ONCOMINE was used to perform a meta‑analysis to compare the mRNA levels of FN1, CXCL8, CXCL10, SPP1, FOXM1, AURKA, ISG15, PLAUR, TPX2, and BIRC5 between OSCC and non-tumor oral tissues. As showed in Figure S3, the mRNA expression levels of FN1 (Figure S3a), CXCL8 (Figure S3b), CXCL10 (Figure S3c), SPP1 (Figure S3d), FOXM1 (Figure S3e), AURKA (Figure S3f), ISG15 (Figure S3g), PLAUR (Figure S3h), TPX2 (Figure S3i), and BIRC5 (Figure S3j) were significantly upregulated in OSCC tissues compared to those in normal oral tissues (P < 0.05). In addition, the median rank of FOXM1 was the highest among the top 10 hub genes in OSCC tissues (Figure S3e). The findings from GEPIA database also demonstrated that the 10 hub genes were significantly higher in OSCC tissues than those of non-tumor tissues (Fig. 3a). These results were consistent with the observed in GEO datasets. Due to the lack of CXCL8, CXCL10, FOXM1 and PLAUR information in The Human Protein Atlas dataset, their protein expression level was not analyzed (Fig. 3b). The results indicated that the protein expressions of AURKA, BIRC5, FN1, ISG15, SPP1 and TPX2 were overexpressed in OSCC compared with the control samples (Fig. 3b).
Genetic alterations and survival analysis of hub genes
Kaplan Meier-plotter website was applied to evaluate the prognostic potential of the 10 hub genes. The results demonstrated that high levels of SPP1 (HR=1.45(1.09-1.92), P=0.01), FN1 (HR=1.49(1.14-1.96), P=0.0038), CXCL8 (HR=1.52(1.13-2.04), P=0.0048), BIRC5 (HR=1.33(1.01-1.76), P=0.043), PLAUR (HR=1.39(1.06-1.82), P=0.016), and AURKA (HR=1.44(1.05-18.2), P=0.022) were associated with poor OS in OSCC patients (Fig. 4). The remaining four genes present similar trends, but not statistically significant. Moreover, the genetic alterations were enquired by cBioPortal. Figure 5a and b showed the alteration state of 10 genes. These 10 genes were changed in 217 (44%) of 496 sequenced patients (496 total), and that FOXM1 and TPX2 were changed most often (15% and 14%), including amplification, missense mutation, and deep deletion. Figure 5c illustrated the network established by the 10 hub genes and their 47 neighbor genes. Besides, drugs targeting the hub genes were presented. According to Fig. 5c, only AURKA, BIRC5, and PLAUR were utilized as chemotherapy targets for cancer treatment presently. Therefore, we supposed that the other seven genes (FN1, CXCL8, CXCL10, SPP1, FOXM1, ISG15, and TPX2) might be the novel targets for OSCC treatment.
Biological functions of FOXM1 and TPX2 in tumors
As FOXM1 and TPX2 were changed most often in OSCC, we thus chose FOXM1 and TPX2 as the target genes to conduct the following studies. Firstly, TIMER website was used to explore the expression levels of FOXM1 in several tumors and corresponding normal tissues. FOXM1 is upregulated in a variety of tumors, including BLCA, BRCA, CHOL, COAD, ESCA, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, PPAD, READ, STAD, THCA, and UCEC (Fig. 6a). This suggests that FOXM1 may also play an oncogenic role in other tumors. Moreover, by using the UALCAN website, we found that FOXM1 is highly expressed in OSCC tissues (Fig. 6b). FOXM1 has differential expression in patients with different tumor stages (Fig. 6c), genders (Fig. 6d), races (Fig. 6e), and molecular subtypes (Fig. 6f). Similarly, TPX2 is upregulated in a variety of human tumors, including BLCA, BRCA, CHOL, COAD, ESCA, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, PPAD, READ, STAD, THCA, and UCEC (Fig. 7a). We also found that TPX2 is highly expressed in OSCC tissues via the UALCAN website (Fig. 7b). In addition, TPX2 has differential expression in patients with different tumor stages (Fig. 7c), genders (Fig. 7d), races (Fig. 7e), and molecular subtypes (Fig. 7f).
FOXM1 and TPX2 act as the immune-associated genes in OSCC
The TIMER website was applied to analyze the relationship between FOXM1, TPX2, and immune cell infiltration. The results showed that both FOXM1 and TPX2 are involved in the infiltration of B cells, CD4+ T cells, CD8+ T cells, Macrophage cells and Neutrophils cells in OSCC (Fig. 8a, c). Moreover, we further analyzed the co-expression relationship of FOXM1, TPX2, and immune checkpoint-related genes CD274, CTLA4, PDCD1, PDCD1LG2, and TOX. Interesting, we found that FOXM1 was markedly correlated with PDCD1, PDCD1LG2, CD274, TOX, and CTLA4 (Fig. 8b). And TPX2 has a significant co-expression relationship with PDCD1LG2, CD274, and TOX (Fig. 8d). These findings suggest that FOXM1 and TPX2 may act as the immune-related therapeutic targets in OSCC.
Construction of the risk gene signature for OSCC patients
To construct a promising risk gene signature using the 261 overlapping DEGs in OSCC, univariate analysis was performed to identify the DEGs related to the prognosis of OSCC patients using the TCGA dataset. Totally 42 overlapping DEGs was markedly associated with the prognosis of OSCC patients (Table S3). A stepwise multivariate Cox regression was then conducted to construct the risk gene signature. Nine candidate genes (TEX101, DSG2, SCG5, ADA, BOC, SCARA5, FST, SOCS1, and STC2) were identified as the significant prognostic indicators for OSCC patients (Table S4). Each OSCC patient was assigned a risk score calculated as follows: risk score = (-0.32654 × expression value of TEX101) + (0.10257 × expression value of DSG2) + (0.15022 × expression value of SCG5) + (0.18621 × expression value of ADA) + (-0.29365 × expression value of BOC) + (-0.26023 × expression value of SCARA5) + (0.10796 × expression value of FST) + (-0.26676 × expression value of SOCS1) + (0.15679 × expression value of STC2).
The OSCC patients were divided into high-risk and low-risk group according to the median value of risk score (Fig. 9a-c). Figure 9a-b shows that the survival time of OSCC patients decreases along with the rising of risk score. As shown in Fig. 9d, the Kaplan-Meier curves shows that OSCC patients with lower risk score present a longer survival time than those with high risk score (P=1.059e-09). Figure 9e shows that the risk score curve presents a good feasibility in predicting the patients’ survival with AUC of 0.685.
The risk score, stage and N were independent prognostic indicators in OSCC
As shown in Fig. 10a, the risk score was significantly associated with the poorer OS in OSCC (HR=1.871, 95% CI: 1.455-2.405, P < 0.001). Moreover, stage (HR=1.895, 95% CI: 1.287-2.790, P < 0.01), N (lymph nodes) (HR=1.398, 95% CI: 1.152-1.697, P < 0.001), and T (primary tumor) (HR=1.389, 95% CI: 1.082-1.783, P < 0.01) were also related to the OS. Then, all these factors were entered into multivariate Cox analysis. The risk score (HR=2.048, 95% CI: 1.546-2.711, P < 0.001), stage (HR=1.706, 95% CI: 1.040-2.799, P < 0.05), and N (HR=1.368, 95% CI: 1.086-1.726, P < 0.01) were still identified as the independent prognostic indicators for worse OS in OSCC patients (Fig. 10b).
Nowadays, the incidence rate of OSCC is still increasing quickly . It is estimated that over 354,864 new cases and 177,384 deaths occurred in 2018 . Compared with the other researches that only explored a single cohort or several genes, our study used several databases to integratedly investigate the key genes, pathways, biomarkers and risk gene signature in OSCC development. In this study, totally 261 overlapping DEGs (135 upregulated DEGs and 126 downregulated DEGs) were identified. KEGG pathway enrichment analysis shows that pathway in cancer, PI3K-Akt signaling pathway, ECM-receptor interaction, focal adhesion, metabolic pathways, serotonergic synapse, and tyrosine metabolism might involve in the OSCC progression. The PPI network contains 260 nodes and 655 edges. Then, two important clusters were picked out, and these two cluster are mainly enriched in cell cycle, DNA replication, cellular senescence, IL-17 signaling pathway, and influenza A. 10 genes were identified as hub genes conforming to the degree of connectivity. We then validated the levels of hub genes in ONCOMINE, TCGA and The Human Protein Atlas database. Six genes (SPP1, FN1, CXCL8, BIRC5, PLAUR, and AURKA) were markedly elevated in OSCC and related with poor prognosis. Additionally, genetic analysis demonstrated that hub genes were changed in about 44% OSCC patients and these genetic alternations include amplification, missense mutation and so on.
Importantly, we also focused on finding the therapeutic targets in OSCC. As shown in Fig. 5c, only AURKA, BIRC5, and PLAUR were served as the chemotherapy targets for cancer treatment currently. Therefore, more investigations and clinical trials are required to explore whether the other seven genes (FN1, CXCL8, CXCL10, SPP1, FOXM1, ISG15, and TPX2) could serve as novel therapeutic targets for OSCC patients. According to the results in TIMER website, FOXM1 and TPX2 were regarded as the potential immunotherapeutic targets with future clinical significance. Thus, our analysis may provide valuable clues for the targeted therapy in OSCC.
From construction of the risk gene signature, we found that the nine-gene based risk model could effectively discriminate the OSCC patients with different outcome(P=1.059e-09), and it presents a good performance in prognosis judgement. In addition, the Cox regression analysis further demonstrated that the risk score based on the gene signature is an independent prognostic indicator with the highest HR value than other factors.
Driven-genes play a crucial role in cancer progression through complex pathways and networks. Similar to the findings of our study, it has been reported that secreted phosphoprotein 1 (SPP1) is a cancer-related gene, which presents clearly upregulated level in many cancers [42,43,44]. Additionally, SPP1 also plays a significant role in extracellular matrix binding . Consistently, SPP1 was enriched in ECM receptor interaction according to KEGG pathway analysis, which plays crucial role in cancer metastasis . Huang et al found that overexpressed SPP1 was linked to carcinogenesis and progression of OSCC . Besides, Fibronectin 1(FN1), predominantly overexpressed in many tumor tissues [47, 48], was also involved in this pathway. Cai et al confirmed that downregulated FN1 could inhibit colorectal carcinogenesis through suppressing proliferation, migration, and invasion . Therefore, SPP1 and FN1 might be potential therapeutic targets for inhibiting OSCC metastasis.
Chemokine (C-X-C motif) ligand 8 (CXCL8) and CXCL10 belong to the chemokine family . Since CXCL8 and CXCL10 integrates with multiple intracellular signaling pathways associated with pro-inflammatory and pro-oncogenic processes, upregulated CXCL8 and CXCL10 are related to carcinogenesis and might predict prognosis of patients [51, 52]. For example, the level of CXCL8 was upregulated in endothelial cells co-cultured with HNSCC, showing that CXCL8 might play a pro-oncogenic role in the pathobiology of tumor cells . Moreover, CXCL8 and CXCL10 are two important modulators in immune response, and they might provide new opportunities for improving immune therapies and enhancing the effectiveness of existing chemotherapies [52, 54].
Plasminogen activator urokinase receptor (PLAUR) is one of glycosyl-phosphatidylinositol (GPI)-anchored membrane proteins . Downregulation of PLAUR could inhibit cancer proliferation and metastasis in several cancers . Consistently, our results showed that PLAUR was considered as a cancer therapeutic target and high expression of PLAUR can predicted poor OS in OSCC patients. Forkhead box protein M1 (FOXM1) is widely participated in the carcinogenesis of several malignances . For example, FOXM1-induced epigenetic signature may serve as ideal biomarkers for early cancer screening in head and neck carcinoma . FOXM1 may also act as a therapeutic target against head and neck carcinoma [59, 60]. Recent study showed that upregulated basal FOXM1 activity predisposes HPV positive HNSCC to WEE1i-induced toxicity . Consistently, we found that FOXM1 is differentially expressed in OSCC patient with different tumor stages, genders, races, and molecular subtypes. Moreover, we were surprised to find that FOXM1 may act as an immune-related therapeutic targets according to the TIMER website. FOXM1 was markedly correlated with PDCD1 (Fig 8b). However, little research has been conducted in this field. Previous study showed that FOXM1 could maintain the dynamic balance between cell apoptosis and proliferation through regulating the essential genes . Thus, FOXM1 may also participate in the PDCD1-regulated cell death. Moreover, the association of FOXM1 with immune cell infiltration may enhance our understanding of the correlationship between FOXM1 and PDCD1. The functions of regulatory T cells (Tregs) in immune regulation is well known, and much is being made of their potential for PDCD1 based therapy . Recently, researchers found that the expression of FOXM1 was positive correlated with FOXP3 (the specific molecular marker of Tregs) and FOXM1 may induce immune suppression via recruiting FOXP3 positive Tregs in cancer therapy . These results highlight the therapeutic potential of FOXM1 in OSCC. Previously, Li et al also found that there was a direct link between PLAUR and FOXM1 . Their results showed that FOXM1 can regulate the level of PLAUR via binding to its promoter. Moreover, the FOXM1-PLAUR axis contributed to colon carcinomas . We therefore supposed that FOXM1-PLAUR signaling might also be promising for designing novel therapeutic drugs for OSCC.
Targeting protein for xenopus kinesin-like protein 2 (TPX2) is involved in the mitotic spindle assembly and cell-cycle progression . Recent studies reported that TPX2 dysregulation was related to the progression of esophageal cancer , hepatocellular carcinoma , and colorectal cancer . Consistently, increased expression of TPX2 was also observed in OSCC tissues in this study. Moreover, TPX2 is differentially expressed in patients with different tumor stages, genders, races, and molecular subtypes according to UALCAN website. We also analyzed role of TPX2 in tumor immunity and found that TPX2 is involved in the infiltration of B cells, CD4+ T cells, CD8+ T cells, Macrophage cells and Neutrophils cells in OSCC. Importantly, TPX2 has a significant co-expression relationship with CD274, PDCD1LG2, and TOX. Since TOX is a significant regulator for T cell differentiation , we believe that TPX2 may participate in OSCC development by regulating TOX molecules.
There were still some limitations in our study. Firstly, our study mainly focused on the biological function of 10 hub genes and did not deeply explore the other DEGs. Therefore, more investigations are required in future. Secondly, because we just utilized ONCOMINE database, TCGA database, and The Human Protein Atlas database to verify the expression of hub genes, the experimental assays are also needed to demonstrate the above results. Finally, the prognostic model should be further validated in large clinical cohort.
Through conducting in silico analyses, 261 overlapping DEGs were identified in OSCC, which were mainly associated with PI3K-Akt signaling pathway, ECM-receptor interaction, focal adhesion, metabolic pathways, serotonergic synapse, and tyrosine metabolism. We identified 10 hub genes (FN1, CXCL8, CXCL10, SPP1, FOXM1, AURKA, ISG15, PLAUR, TPX2, and BIRC5), which can act as available targets for OSCC treatment. Additionally, we constructed a nine-gene prediction model that can be utilized as the prognostic tool in OSCC. Our study is helpful to explore the carcinogenesis of OSCC. Moreover, these findings can improve the prognosis judgement and targeted therapy of OSCC.
Availability of data and materials
All analyzed data related to this paper are included in this paper.
Differentially expressed genes
Search tool for the retrieval of interacting genes
Oral squamous cell carcinoma
Head and neck squamous cell carcinoma
Forkhead box protein M1
Targeting protein for xenopus kinesin-like protein 2
Gene Expression Omnibus
The Cancer Genome Atlas
Kyoto Encyclopedia of Genes and Genomes
Receiver operating characteristic
The area under the ROC curve
Miller KD, Siegel RL, Lin CC, Mariotto AB, Kramer JL, Rowland JH, Stein KD, Alteri R, Jemal A. Cancer treatment and survivorship statistics, 2016. CA Cancer J Clin. 2016;66(4):271–89.
Jiang S, Dong Y. Human papillomavirus and oral squamous cell carcinoma: A review of HPV-positive oral squamous cell carcinoma and possible strategies for future. CurrProbl Cancer. 2017;41(5):323–7.
Eleftheriadou A, Chalastras T, Ferekidou E, Yiotakis I, Kyriou L, Tzagarakis M, Ferekidis E, Kandiloros D. Association between squamous cell carcinoma of the head and neck and serum folate and homocysteine. Anticancer Res. 2006;26(3b):2345–8.
Lo AK, Lo KW, Tsao SW, Wong HL, Hui JW, To KF, Hayward DS, Chui YL, Lau YL, Takada K, et al. Epstein-Barr virus infection alters cellular signal cascades in human nasopharyngeal epithelial cells. Neoplasia. 2006;8(3):173–80.
Zhang X, Feng H, Li D, Liu S, Amizuka N, Li M. Identification of Differentially Expressed Genes Induced by Aberrant Methylation in Oral Squamous Cell Carcinomas Using Integrated Bioinformatic Analysis. Int J Mol Sci. 2018;19(6):1698.
Kamps R, Brandao RD, Bosch BJ, Paulussen AD, Xanthoulea S, Blok MJ, Romano A. Next-Generation Sequencing in Oncology: Genetic Diagnosis, Risk Prediction and Cancer Classification. Int J Mol Sci. 2017;18(2):308.
Sipari H, Rantala-Ylinen A, Jokela J, Oksanen I, Sivonen K. Development of a chip assay and quantitative PCR for detecting microcystin synthetase E gene expression. Appl Environ Microbiol. 2010;76(12):3797–805.
Kulasingam V, Diamandis EP. Strategies for discovering novel cancer biomarkers through utilization of emerging technologies. Nat ClinPractOncol. 2008;5(10):588–99.
Ishigami T, Uzawa K, Higo M, Nomura H, Saito K, Kato Y, Nakashima D, Shiiba M, Bukawa H, Yokoe H, et al. Genes and molecular pathways related to radioresistance of oral squamous cell carcinoma cells. Int J Cancer. 2007;120(10):2262–70.
Guo Y, Bao Y, Ma M, Yang W. Identification of Key Candidate Genes and Pathways in Colorectal Cancer by Integrated Bioinformatical Analysis. Int J Mol Sci. 2017;18(4):722.
Clough E, Barrett T. The Gene Expression Omnibus Database. Methods MolBiol. 2016;1418:93–110.
Yang W, Mai J, Zhou W, Li Z, Zhou X, Cao B, Zhang Y, Liu J, Yang Z, Zhang H, et al. Identification of hub genes and outcome in colon cancer based on bioinformatics analysis. Cancer Manage Res. 2019;11:323–38.
Ambatipudi S, Gerstung M, Pandey M, Samant T, Patil A, Kane S, Desai RS, Schaffer AA, Beerenwinkel N, Mahimkar MB. Genome-wide expression and copy number analysis identifies driver genes in gingivobuccal cancers. Genes Chromosomes Cancer. 2012;51(2):161–73.
Chen C, Mendez E, Houck J, Fan W, Lohavanichbutr P, Doody D, Yueh B, Futran ND, Upton M, Farwell DG, et al. Gene expression profiling identifies genes predictive of oral squamous cell carcinoma. Cancer Epidemiol Biomarkers Prev. 2008;17(8):2152–62.
Oghumu S, Knobloch TJ, Terrazas C, Varikuti S, Ahn-Jarvis J, Bollinger CE, Iwenofu H, Weghorst CM, Satoskar AR. Deletion of macrophage migration inhibitory factor inhibits murine oral carcinogenesis: Potential role for chronic pro-inflammatory immune mediators. Int J Cancer. 2016;139(6):1379–90.
Lee CH, Wong TS, Chan JY, Lu SC, Lin P, Cheng AJ, Chen YJ, Chang JS, Hsiao SH, Leu YW, et al. Epigenetic regulation of the X-linked tumour suppressors BEX1 and LDOC1 in oral squamous cell carcinoma. J Pathol. 2013;230(3):298–309.
Guo T, Ma H, Zhou Y. Bioinformatics analysis of microarray data to identify the candidate biomarkers of lung adenocarcinoma. PeerJ. 2019;7:e7313.
Giulietti M, Occhipinti G, Principato G, Piva F. Weighted gene co-expression network analysis reveals key genes involved in pancreatic ductal adenocarcinoma development. Cell Oncol (Dordrecht). 2016;39(4):379–88.
Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics (Oxford, England). 2003;19(2):185–93.
Chen J, Yu K, Zhong G, Shen W. Identification of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of clear cell renal carcinoma. Cancer Cell Int. 2020;20:157.
Du F, Sun Z, Jia J, Yang Y, Yu J, Shi Y, Jia B, Zhao J, Zhang X. Development and Validation of an Individualized Nomogram for Predicting Survival in Patients with Esophageal Carcinoma after Resection. J Cancer. 2020;11(14):4023–9.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45(D1):D353-d361.
Huang DW, Sherman BT, Tan Q, Kir J, Liu D, Bryant D, Guo Y, Stephens R, Baseler MW, Lane HC, et al. DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 2007;35(Web Server issue):W169-175.
Pathan M, Keerthikumar S, Ang CS, Gangoda L, Quek CY, Williamson NA, Mouradov D, Sieber OM, Simpson RJ, Salim A, et al. FunRich: An open access standalone functional enrichment and interaction network analysis tool. Proteomics. 2015;15(15):2597–601.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.
Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, Lin J, Minguez P, Bork P, von Mering C, et al. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013;41(Database issue):D808-815.
Yang W, Zhao X, Han Y, Duan L, Lu X, Wang X, Zhang Y, Zhou W, Liu J, Zhang H, et al. Identification of hub genes and therapeutic drugs in esophageal squamous cell carcinoma based on integrated bioinformatics strategy. Cancer Cell Int. 2019;19:42
Ni M, Liu X, Wu J, Zhang D, Tian J, Wang T, Liu S, Meng Z, Wang K, Duan X, et al. Identification of Candidate Biomarkers Correlated With the Pathogenesis and Prognosis of Non-small Cell Lung Cancer via Integrated Bioinformatics Analysis. Front Genet. 2018;9:469.
Rhodes DR, Kalyana-Sundaram S, Mahavisno V, Varambally R, Yu J, Briggs BB, Barrette TR, Anstet MJ, Kincead-Beal C, Kulkarni P. et al. Oncomine 3.0: genes, pathways, and networks in a collection of 18,000 cancer gene expression profiles. Neoplasia. 2007;9(2):166–80.
Song X, Du R, Gui H, Zhou M, Zhong W, Mao C, Reports MJJO. Identification of potential hub genes related to the progression and prognosis of hepatocellular carcinoma through integrated bioinformatics analysis. Oncol Rep. 2020;43(1):133–46.
Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017;45(W1):W98-w102.
Ponten F, Schwenk JM, Asplund A, Edqvist PH. The Human Protein Atlas as a proteomic resource for biomarker discovery. J Intern Med. 2011;270(5):428–46.
Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, et al. ThecBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2(5):401–4.
Gyorffy B, Surowiak P, Budczies J, Lanczky A. Online survival analysis software to assess the prognostic value of biomarkers using transcriptomic data in non-small-cell lung cancer. PLoS One. 2013;8(12):e82241.
Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017;77(21):e108–10.
Chandrashekar DS, Bashel B, Balasubramanya SAH, Creighton CJ, Ponce-Rodriguez I, Chakravarthi B, Varambally S. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia. 2017;19(8):649–58.
Liu Y, Wu L, Ao H, Zhao M, Leng X, Liu M, Ma J, Zhu J. Prognostic implications of autophagy-associated gene signatures in non-small cell lung cancer. Aging (Albany NY). 2019;11(23):11440–62.
Lakshminarayana S, Augustine D, Rao RS, Patil S, Awan KH, Venkatesiah SS, Haragannavar VC, Nambiar S, Prasad K. Molecular pathways of oral cancer that predict prognosis and survival: A systematic review. J Carcinog. 2018;17:7.
Weyandt JD, Thompson CB, Giaccia AJ, Rathmell WK. Metabolic Alterations in Cancer and Their Potential as Therapeutic Targets. Am SocClinOncolEduc Book Am SocClinOncol Meet. 2017;37:825–32.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424.
Zhang W, Fan J, Chen Q, Lei C, Qiao B, Liu Q. SPP1 and AGER as potential prognostic biomarkers for lung adenocarcinoma. OncolLett. 2018;15(5):7028–36.
Choe EK, Yi JW, Chai YJ, Park KJ. Upregulation of the adipokine genes ADIPOR1 and SPP1 is related to poor survival outcomes in colorectal cancer. J SurgOncol. 2018;117(8):1833–40.
Huang CF, Yu GT, Wang WM, Liu B, Sun ZJ. Prognostic and predictive values of SPP1, PAI and caveolin-1 in patients with oral squamous cell carcinoma. Int J ClinExpPathol. 2014;7(9):6032–9.
Fudge NJ, Mearow KM. Extracellular matrix-associated gene expression in adult sensory neuron populations cultured on a laminin substrate. BMC Neurosci. 2013;14:15.
Xu C, Sun L, Jiang C, Zhou H, Gu L, Liu Y, Xu Q. SPP1, analyzed by bioinformatics methods, promotes the metastasis in colorectal cancer by activating EMT pathway. Biomed Pharmacother. 2017;91:1167–77.
Ouchi K, Miyachi M, Tsuma Y, Tsuchiya K, Iehara T, Konishi E, Yanagisawa A, Hosoi H. FN1: a novel fusion partner of ALK in an inflammatory myofibroblastic tumor. Pediatr Blood Cancer. 2015;62(5):909–11.
Zhang H, Sun Z, Li Y, Fan D, Jiang H. MicroRNA-200c binding to FN1 suppresses the proliferation, migration and invasion of gastric cancer cells. Biomed Pharmacother. 2017;88:285–92.
Cai X, Liu C, Zhang TN, Zhu YW, Dong X, Xue P. Down-regulation of FN1 inhibits colorectal carcinogenesis by suppressing proliferation, migration, and invasion. J Cell Biochem. 2018;119(6):4717–28.
Lee NH, Nikfarjam M, He H. Functions of the CXC ligand family in the pancreatic tumor microenvironment. Pancreatology. 2018;18(7):705–16.
Liu Q, Li A, Tian Y, Wu JD, Liu Y, Li T, Chen Y, Han X, Wu K. The CXCL8-CXCR1/2 pathways in cancer. Cytokine Growth Factor Rev. 2016;31:61–71.
Tokunaga R, Zhang W, Naseem M, Puccini A, Berger MD, Soni S, McSkane M, Baba H, Lenz HJ. CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation - A target for novel cancer therapy. Cancer Treat Rev. 2018;63:40–7.
Neiva KG, Zhang Z, Miyazawa M, Warner KA, Karl E, Nor JE. Cross talk initiated by endothelial cells enhances migration and inhibits anoikis of squamous cell carcinoma cells through STAT3/Akt/ERK signaling. Neoplasia. 2009;11(6):583–93.
Russo RC, Garcia CC, Teixeira MM, Amaral FA. The CXCL8/IL-8 chemokine family and its receptors in inflammatory diseases. Expert Rev ClinImmunol. 2014;10(5):593–619.
Ploug M, Ronne E, Behrendt N, Jensen AL, Blasi F, Dano K. Cellular receptor for urokinase plasminogen activator. Carboxyl-terminal processing and membrane anchoring by glycosyl-phosphatidylinositol. J BiolChem. 1991;266(3):1926–33.
Zhou J, Kwak KJ, Wu Z, Yang D, Li J, Chang M, Song Y, Zeng H, Lee LJ, Hu J, et al. PLAUR Confers Resistance to Gefitinib Through EGFR/P-AKT/Survivin Signaling Pathway. Cell PhysiolBiochem. 2018;47(5):1909–24.
Gartel AL. FOXM1 in Cancer: Interactions and Vulnerabilities. Cancer Res. 2017;77(12):3135–9.
Teh MT, Gemenetzidis E, Patel D, Tariq R, Nadir A, Bahta AW, Waseem A, Hutchison IL. FOXM1 induces a global methylation signature that mimics the cancer epigenome in head and neck squamous cell carcinoma. PLoS One. 2012;7(3):e34329.
Chen Y, Liu Y, Ni H, Ding C, Zhang X, Zhang Z. FoxM1 overexpression promotes cell proliferation and migration and inhibits apoptosis in hypopharyngeal squamous cell carcinoma resulting in poor clinical prognosis. Int J Oncol. 2017;51(4):1045–54.
Luo W, Gao F, Li S, Liu L. FoxM1 Promotes Cell Proliferation, Invasion, and Stem Cell Properties in Nasopharyngeal Carcinoma. Front Oncol. 2018;8:483.
Diab A, Gem H. FOXM1 drives HPV+ HNSCC sensitivity to WEE1 inhibition. Proc Natl Acad Sci U S A. 2020;117(45)28287–96.
Liao GB, Li XZ, Zeng S, Liu C, Yang SM, Yang L, Hu CJ, Bai JY. Regulation of the master regulator FOXM1 in cancer. Cell Commun Signal. 2018;16(1):57.
Wang D, Liu J, Liu S, Li W. Identification of Crucial Genes Associated With Immune Cell Infiltration in Hepatocellular Carcinoma by Weighted Gene Co-expression Network Analysis. Front Genet. 2020;11:342.
González-Navajas JM, Fan DD, Yang S, Yang FM, Lozano-Ruiz B, Shen L, Lee J. The Impact of Tregs on the Anticancer Immunity and the Efficacy of Immune Checkpoint Inhibitor Therapies. Front Immunol. 2021;12:625783.
Li X, Ma K, Song S, Shen F, Kuang T, Zhu Y, Liu Z. Tight correlation between FoxM1 and FoxP3+ Tregs in gastric cancer and their clinical significance. ClinExp Med. 2018;18(3):413–20.
Li D, Wei P, Peng Z, Huang C, Tang H, Jia Z, Cui J, Le X, Huang S, Xie K. The critical role of dysregulated FOXM1-PLAUR signaling in human colon cancer progression and metastasis. Clin Cancer Res. 2013;19(1):62–72.
Warner SL, Stephens BJ, Nwokenkwo S, Hostetter G, Sugeng A, Hidalgo M, Trent JM, Han H, Von Hoff DD. Validation of TPX2 as a potential therapeutic target in pancreatic cancer cells. Clin Cancer Res. 2009;15(21):6519–28.
Hou Y, Liu H, Pan W. Knockdown of circ_0003340 induces cell apoptosis, inhibits invasion and proliferation through miR-564/TPX2 in esophageal cancer cells. Exp Cell Res. 2020;394(2):112142.
Ouyang G, Yi B, Pan G, Chen X. A robust twelve-gene signature for prognosis prediction of hepatocellular carcinoma.Cancer Cell Int. 2020;20:207.
Taherdangkoo K, KazemiNezhad SR, Hajjari MR, TahmasebiBirgani M. miR-485-3p suppresses colorectal cancer via targeting TPX2. BratislLekListy. 2020;121(4):302–7.
Scott AC, Dündar F, Zumbo P, Chandran SS, Klebanoff CA, Shakiba M, Trivedi P, Menocal L, Appleby H, Camara S, et al. TOX is a critical regulator of tumour-specific T cell differentiation. Nature. 2019;571(7764):270–4.
We thank the GEO, DAVID, STRING, GEPIA2, cBioPortal, TCGA and TIMER databases for providing their platforms.
This study was supported in part by grant from the Scientific Foundation of Shaanxi Province (S2019ZDLSF01-02-01; 2018SF-240) and grant from the National Clinical Research Center for Digestive Diseases (2015BAI13B07) and Boost Program of Xijing Hospital (XJZT18ML06).
Ethics approval and consent to participate
Consent for publication
All the authors have consented for the publication.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
GO and KEGG pathway terms of DEGs in OSCC. Figure S2. PPI network of 10 hub genes. Figure S3. Validation of hub genes in ONCOMINE database.
The enriched pathways in Cluster 1. Table S2. The enriched pathways in Cluster 2. Table S3. Univariate Cox regression analysis of OS in OSCC patients. Table S4. Multivariate Cox regression analysis of OS in OSCC patients.
About this article
Cite this article
Yang, W., Zhou, W., Zhao, X. et al. Prognostic biomarkers and therapeutic targets in oral squamous cell carcinoma: a study based on cross-database analysis. Hereditas 158, 15 (2021). https://doi.org/10.1186/s41065-021-00181-1
- Oral squamous cell carcinoma
- Bioinformatics analysis
- Differentially expressed genes