Abstract

   There are a few biological functions or phenomenon which are
   universally associated with majority of the cancers and hypoxia and
   immune systems are among them. Hypoxia often occurs in most of the
   cancers which helps the cells in adapting different responses with
   respect to the normal cells which may be the activation of signaling
   pathways which regulate proliferation, angiogenesis, and cell death.
   Similar to it, immune signaling pathways are known to play critical
   roles in cancers. Moreover, there are a number of genes which are known
   to be associated with these hypoxia and immune system and appear to
   direct affect the tumor growth and propagations. Cancer is among the
   leading cause of death and oral cancer is the tenth-leading cause due
   to cancer death. In this study, we were mainly interested to understand
   the impact of alteration in the expression of hypoxia and immune
   system-related genes and their contribution to head and neck squamous
   cell carcinoma. Moreover, we have collected the genes associated with
   hypoxia and immune system from the literatures. In this work, we have
   performed meta-analysis of the gene and microRNA expression and
   mutational datasets obtained from public database for different grades
   of tumor in case of oral cancer. Based on our results, we conclude that
   the critical pathways which dominantly enriched are associated with
   metabolism, cell cycle, immune system and based on the survival
   analysis of the hypoxic genes, we observe that the potential genes
   associated with head and neck squamous cell carcinoma and its
   progression are STC2, PGK1, P4HA1, HK1, SPIB, ANXA5, SERPINE1, HGF,
   PFKM, TGFB1, L1CAM, ELK4, EHF, and CDK2.

   Subject terms: Cancer genomics, Oral cancer, Computational biology and
   bioinformatics, Immunology

Introduction

   Hypoxia often occurs in cancer and helps the cells in adapting
   different responses than the normal cells such as the triggering of
   signaling pathways regulating critical biological processes
   (proliferation, angiogenesis, and cell death or apoptosis)^[34]1. So
   far from the previous work, a number of genes associated with these
   processes and functions have been explored and investigated. Similarly,
   different forms of cancer require different immune systems and
   associated signalling pathways, and the immune signalling network (ISN)
   may be a major component in cancer genesis and progression. Although it
   has been established that cancers, including head and neck cancer, are
   immunogenic tumours for which immunotherapy is aggressively sought by
   targeting immunological checkpoints, an immune-based prognostic
   signature remains a viable option^[35]2–[36]4. Several prior
   works^[37]2–[38]4 propose pathway-level knowledge and analytic
   methodologies, as well as Hansen and Iyengar’s^[39]4 computational
   strategy to bridge the gap between precision medicine and systems
   treatments. Comprehension and unravelling comprehensive and minute
   understanding of cell phenotypes and disease pathophysiology remains a
   basic problem, as do the molecular mechanisms that lead to disease
   initiation and oral cancer progression.

   Oral cancer is the 10th most prevalence cancer globally and in general
   classified as head and neck squamous cell carcinoma (HNSCC). It is a
   malignant neoplasia which arises in oral cavity of lip, tongue,
   gingiva, mouth floor and glands^[40]5–[41]7 which may originate by a
   number of factors such as genetic alterations, gene expression
   alterations, and mutations^[42]5,[43]8–[44]13. Furthermore, more
   factors may act as the potential cause of such cancer and one of the
   cause is HPV (Human Pappiloma Virus) infection and this virus is
   non-enveloped icosahedral capsid with circular double standard DNA
   which majorly cause cervical cancer in human^[45]12,[46]14,[47]15. In
   addition, there exit potential difference between HPV-induced oral
   cancers than that of HPV-negative (oral) tumors in terms of the
   clinical response and finally the overall survival rates^[48]16,[49]17.
   In the previous study, there are lots of work which have been performed
   for the study of oral cancer such as mutational, gene, and miRNA
   expression profiling, epigenetic changes, and proteomics for
   HNSCC^[50]5,[51]12,[52]18,[53]19.

   When we’re looking for a profound understanding of something, from a
   computational study to a therapeutic method^[54]20–[55]22. It provides
   a ray of hope for a revolutionary diagnostic technique. Thus, in HNSCC,
   the hypoxia and immune-based prognostic signatures maintain a
   diagnostic potential that can be further explored and examined. We
   chose a publically available gene expression dataset for this purpose
   and evaluated the data with the goal of understanding how signaling
   networks and their components are relevant to the immune system. In
   this study, our goal was to understand the impact of alteration in the
   expression of hypoxia and immune system-related genes and their
   contribution to head and neck cancer^[56]23,[57]24. For this purpose,
   we have collected the hypoxia-associated genes based on the literature
   related to diverse biological processes and functions and have also
   collected the mutational and expression (both microRNA and gene)
   datasets from Gene Expression Omnibus (GEO) freely accessible public
   database ([58]http://www.ncbi.nlm.nih.gov/geo/) for which we have
   performed comparative analysis and the clinical relevance i.e.,
   survival analysis^[59]25–[60]27.

   Based on our work, we observe that there are certain sets of genes
   which are always differentially expressed irrespective of the stages
   and similar to it there are a number of pathways which are potentially
   altered in result to the differential gene expression patterns. Based
   on our results, we conclude that the critical pathways which are
   dominantly enriched are associated with metabolism, cell cycle, immune
   system.

Materials and methods

   In the first step, we have selected the data of our interest for
   mutation, gene, and miRNA expression analysis. For the gene expression
   and miRNA expression datasets, the samples have been analyzed by using
   the inbuilt tool GEO2R^[61]28,[62]29 and the mutational dataset which
   have been obtained from TCGA database have been analyzed from in-house
   MATLAB code. For pathway enrichment analysis, the similar protocols
   have been followed as per DAVID and panther databases^[63]30–[64]33.
   For differential gene expression analysis, we have compared the tumor
   samples with normal samples to generate differentially expressed genes
   and miRNAs lists and to generate the list of mutated genes the
   threshold has been set i.e., 5% of the samples showing the mutation for
   specific genes. After preparing the lists (overexpressed genes and
   microRNAs) and the mutated genes, we proceed for our goal which is to
   understand the expression and mutational patterns^[65]10,[66]34 and its
   inferred functions^[67]33,[68]34. All the pathways have p-values less
   than 0.05, blue color means highest p-values and the yellow color for
   lowest p-value. For generating DEGs network, FunCoup2.0^[69]35 has been
   used for all the networks throughout the work and cytoscape^[70]36 has
   been used for network visualization. For most of our coding and
   calculations MATLAB has been used. FunCoup predicts four different
   classes of functional coupling or associations such as protein
   complexes, protein–protein physical interactions, metabolic, and
   signaling pathways^[71]35.

Results

Gene expression profiling reveals that there are selected sets of pathways
are mostly affected as a result of alteration in gene expression

   In the first step, the major goal was to understand the gene expression
   patterns and profiles between different tumor stages for which we have
   used the oral cancer expression dataset from GEO [72]GSE84846^[73]7.
   The dataset contains the samples which have tongue squamous cells
   carcinoma cells of both male and female of different ages and different
   stages from stage I to IV. We have performed comparative analysis of
   grade I with grade II, III, and IV to investigate the evolved DEGs and
   altered functions from grade I to IV and the total number of samples
   were 99. Here, we observe that the number of DEGs is comparatively low
   for grades I and II i.e., 24, for grades I and III is 40 and for grades
   I and IV is 175 and number of shared genes are quite low (Fig. [74]1a)
   and alterations in gene expression pattern increases exponentially from
   grade I to IV. Irrespective of the tumor location and the number of
   DEGs there are 59 pathways which are commonly enriched even in
   combination with the dataset [75]GSE31056 (Fig. [76]1b,c). All these 59
   enriched pathways have been displayed with their respective p-values
   under different conditions (Fig. [77]1d) and more details have been
   presented as Supplementary Data [78]1. Further, we have also compared
   this data with another dataset ([79]GSE37991)^[80]37–[81]41 in
   Supplementary Data [82]2. Here, we also observe that both the datasets
   share a large number of vital DEGs and the pathways.

Figure 1.

   [83]Figure 1
   [84]Open in a new tab

   Differential gene expression profiling of different grades of oral
   cancer. (a) Venn diagram for differentially expressed genes, inferred,
   and enriched pathways; (b) Plot to show the overall number of DEGs,
   inferred, and enriched pathways.

   We have compared gene expression pattern of overlapped genes in HPV
   infected samples from grade I–IV (Fig. [85]2a–c). Here, we tried to
   present a list of genes (table/supplementary) irrespective of tumor
   grades, at the same time it also shows dissimilarity in level of gene
   expression and functional impact. Figure [86]2d states that Grade I and
   IV are having similar pattern of the expression, grade II shares their
   pattern with all four grades of HPV infected tumor whereas grade III
   showed distinct from all other grades (Fig. [87]2a–c). It claims that
   grade II is important stage where all HPV mediated oncogenic components
   are expressed. In Fig. [88]2d, we have inserted a box with black line
   which is just to show the reverse behavior in terms of expression for
   the selected genes grade 3 versus all other conditions. Overlapped,
   DEGs of oral tumor also showed difference in functional effect.
   Figure [89]2e shows that there are 32 pathways which are mainly altered
   in HPV-infected oral cancer and those p-value were compared with
   altered pathways of non-infected oral cancer (GEO datasets) pathways,
   it revealed that few of the pathways such as circadian entrainment,
   arginine and proline metabolism, butanoate metabolism, PI3K-AKT, cell
   cycle, TGF-beta, cAMP, neuro active ligand receptor interaction are
   highly altered in HPV-infected tumor than non-infected oral cancer
   which implies enhanced vulnerability of HPV infection in oral cancer.
   Here, we also observe that functional effect and gene expression
   followed the same pattern (Fig. [90]2d,e). Thus, we conclude that
   different grades of oral tumor lead to diverse impact on gene
   expression pattern and their functions.

Figure 2.

   [91]Figure 2
   [92]Open in a new tab

   Differentially expressed genes and the enriched pathways for HPV
   infected oral cancer. Venn diagram for the different combinations of
   DEGs and the enriched pathways for (a,b) [93]GSE31056 and (c) NGS
   dataset. (d) Heatmap and cluster for the 152 commonly DEGs in case of
   NGS data. (e) Commonly enriched 32 pathways for NGS dataset and the
   p-values for these pathways including the array dataset.

Higher mutations leads to the potential change in critical biological
functions associated with oral cancer

   For mutational profiling, we have used the datasets from TCGA database
   which contains 530 samples (Head and Neck Squamous Cell Carcinoma,
   Firehose Legacy) and from here, we have selected the mutated genes
   which appear minimum in 5% of the samples and performed the pathway
   enrichment analysis where PI3K-Akt signaling, focal adhesion, thyroid
   hormone signaling, calcium signaling, cAMP signaling, FoxO signaling,
   phospholipase-d signaling, cell cycle, ubiquitin mediated proteolysis,
   apelin signaling, long-term potentiation, oxytocin signaling, Longevity
   regulating, ECM-receptor interaction, circadian entrainment, estrogen
   signaling, and melanogenesis (Fig. [94]3a) are among the enriched
   pathways for the selected genes which were mutated in minimum of the 5%
   of the samples. Furthermore, the top mutated (≥ 10%) genes and observe
   that TP53, TTN, FAT1, CDKN2A, FRG1BP, CSMD3, MUC16, PIK3CA, SYNE1,
   NOTCH1, LRP1B, KMT2D, PCLO, FLG, DNAH5, USH2A, NSD1, RYR2, PKHD1L1,
   XIRP2, CASP8, SI, and AHNAK (Fig. [95]3b) are among the highly mutated
   genes. Majority of these genes are well known to be associated with a
   number of cancers including the head and neck cancer and the similar
   case is with the enriched pathways. Furthermore, we have drawn a venn
   diagram to look over the commonly and specific altered pathways both
   because of altered expression or mutations (Fig. [96]3c). PI3K-Akt
   signaling, cAMP signaling, Focal adhesion, Calcium signaling, Oxytocin
   signaling, Apelin signaling, ECM-receptor interaction, and thyroid
   hormone signaling are those pathways which are commonly altered in
   terms of gene over expression and mutations which gives more
   significance to these pathways while circadian entrainment, cell cycle,
   phospholipase-d signaling, longevity regulating pathway, melanogenesis,
   FoxO signaling, estrogen signaling, long-term potentiation, and
   ubiquitin mediated proteolysis signaling pathways are specifically
   altered due to mutation. There are 16 signaling pathways which are
   exclusively altered due to over expression and some of them are insulin
   signaling, CAMs, TGF, TNF, tight junction, cGMP-PKG, phagosome
   signaling pathways and more.

Figure 3.

   [97]Figure 3
   [98]Open in a new tab

   Mutational profiling and functional impact in oral cancer. (a) We have
   performed pathway enrichment analysis for those genes which appear to
   have more than 5% mutation for the selected dataset from TCGA database;
   (b) Genes with mutations ≥ 10%; (c) Comparison of the altered functions
   with respect to mutations and differential expression; (d) p-value for
   Kaplan-Meyer plots after survival analysis of the hypoxic genes in case
   of head and neck cancer.

   Furthermore, we have also investigated the clinical relevance of the
   hypoxic genes by performing survival analysis (Kaplan-Meyer plot) and
   plot the heatmap of the p-values of all those genes (Fig. [99]3d) which
   appear significant (p-value < 0.05). Overall, 95 genes appear
   significant in case of head and neck cancer and STC2, PGK1, P4HA1, HK1,
   SPIB, ANXA5, SERPINE1, HGF, PFKM, TGFB1, L1CAM, ELK4, EHF, and CDK2
   appear to be highly significant which have p-values even lower than or
   equal to 0.0013. After analyzing the clinical relevance, we have
   fetched the inferred pathways for these top-ranked genes and we observe
   that these genes not only relate with the hypoxic condition but also a
   number of fully those pathways which directly affect the tumor
   initiation, propagation, and growth as shown in Table [100]1. Moreover,
   these top-ranked genes have been processed for detailed clinical
   relevance for which their overexpression have been checked in patients
   samples in TCGA database (Fig. [101]4) for all the four grades and the
   percentage of the patients with overexpressed genes have been shown.

Table 1.

   Top hypoxic genes with p-values ≤ 0.0013 and the associated pathways.
   Genes    Pathways
   HGF      Cytokine-cytokine_receptor_interaction
   HGF      Focal_adhesion
   HGF      Pathways_in_cancer
   HGF      Renal_cell_carcinoma
   HGF      Melanoma
   PGK1     Glycolysis_/_gluconeogenesis
   PGK1     Carbon_fixation_in_photosynthetic_organisms
   TGFB1    MAPK_signaling_pathway
   TGFB1    Cytokine-cytokine_receptor_interaction
   TGFB1    Cell_cycle
   TGFB1    TGF-beta_signaling_pathway
   TGFB1    Leishmaniasis
   TGFB1    Chagas_disease
   TGFB1    Pathways_in_cancer
   TGFB1    Colorectal_cancer
   TGFB1    Renal_cell_carcinoma
   TGFB1    Pancreatic_cancer
   TGFB1    Chronic_myeloid_leukemia
   SERPINE1 p53_signaling_pathway
   SERPINE1 Complement_and_coagulation_cascades
   P4HA1    Arginine_and_proline_metabolism
   CDK2     Cell_cycle
   CDK2     Oocyte_meiosis
   CDK2     p53_signaling_pathway
   CDK2     Progesterone-mediated_oocyte_maturation
   CDK2     Prostate_cancer
   CDK2     Small_cell_lung_cancer
   PFKM     Glycolysis_/_gluconeogenesis
   PFKM     Pentose_phosphate_pathway
   PFKM     Fructose_and_mannose_metabolism
   PFKM     Galactose_metabolism
   PFKM     Insulin_signaling_pathway
   HK1      Glycolysis_/_gluconeogenesis
   HK1      Fructose_and_mannose_metabolism
   HK1      Galactose_metabolism
   HK1      Starch_and_sucrose_metabolism
   HK1      Amino_sugar_and_nucleotide_sugar_metabolism
   HK1      Streptomycin_biosynthesis
   HK1      Insulin_signaling_pathway
   HK1      Type_II_diabetes_mellitus
   ELK4     MAPK_signaling_pathway
   L1CAM    Axon_guidance
   L1CAM    Cell_adhesion_molecules_(CAMs)
   TGFB1    Hippo_Signaling_Pathway
   SERPINE1 Hippo_Signaling_Pathway
   HGF      Ras_signaling_pathway
   HGF      Rap1_signaling_pathway
   SERPINE1 Apelin_signaling_pathway
   HK1      HIF-1_signaling_pathway
   SERPINE1 HIF-1_signaling_pathway
   PGK1     HIF-1_signaling_pathway
   CDK2     FoxO_signaling_pathway
   TGFB1    FoxO_signaling_pathway
   CDK2     PI3K-Akt_signaling_pathway
   HGF      PI3K-Akt_signaling_pathway
   PFKM     AMPK_signaling_pathway
   TGFB1    Osteoclast_differentiation
   [102]Open in a new tab

Figure 4.

   [103]Figure 4
   [104]Open in a new tab

   Clinical relevance. Clinical Relevance for the top-ranked genes (based
   on connectivity of the genes within the network generated through
   network database) and respective inferred pathways. p-value represents
   the clinical significance in terms of survival analysis and the TCGA
   database and cBioPortal have been used.

Differential microRNAs expression also potentially impact the cancer
associated functions

   After analyzing the gene expression and mutational profiling, we have
   performed miRNA expression profiling and for this purpose, the dataset
   was collected from GEO database i.e., [105]GSE31227 and the platform
   was [106]GPL9770^[107]19. This dataset contains 15 patient surgical
   margin as controls and 15 patient Oral Squamous Cell Carcinoma (OSCC)
   and we found that there are 46 miRNAs out of 739 miRNAs which are
   overexpressed in case of head and neck cancer (Fig. [108]5a) and in
   terms of functions the most affected biological pathways are thyroid
   cancer, pathways in cancer, pancreatic cancer, Foxo signaling, chronic
   myeloid leukemia, HIF1 signaling and more and most of these pathways
   appears to strongly associated with head and neck cancer and oral
   cancer (Fig. [109]5b). Here, it can be clearly seen that the pathways
   enriched in DEGs and the miRNA pathways list have a number of common
   pathways such as HIF-1 signaling, Ras, MAPK, immune system associated
   pathways, and directly cancer-associated pathways.

Figure 5.

   [110]Figure 5
   [111]Open in a new tab

   miRNA expression profiling and their functional impact in head and neck
   cancer. (a) Differentially expressed miRNAs and (b) the enriched
   pathways with the respective miRNAs.

Discussion

   Hypoxia and ISN often help the cells in adapting different responses
   than the normal cells such as the triggering of signaling pathways
   regulating critical biological processes (proliferation, angiogenesis,
   and cell death or apoptosis)^[112]15,[113]20,[114]42–[115]49. These two
   processes have been explored in a number of human diseases and
   different from the previous works, here these two systems have been
   explored simultaneously in HNSCC and also the role of microRNAs have
   been analyzed with their effect on biological processes and functions
   and furthermore the survival analysis^[116]12,[117]50–[118]52 has been
   analyzed for the hypoxic genes. The reason to focus on hypoxia and
   mainly the ISN was that the associated signaling pathways with these
   two systems are considered as a master regulator for cancer initiation
   and progression and has also been proven that HNSCC is an immunogenic
   tumor and immunotherapy is strongly pursued through targeting on the
   immune checkpoints the immune based prognostic signature remains a
   potential that can be applied^[119]2–[120]4. In the previous
   works^[121]2–[122]4, pathway-level understanding and analysis
   approaches has also been presented which present computational approach
   to bridge between precision medicine and systems therapeutics.

   In this study, the main focus of the study was to understand the
   expression pattern of both the genes and the microRNAs and the
   mutational profiling followed by the survival analysis for HNSCC for
   which the datasets have been utilized from the GEO and the TCGA
   database. After the expression and mutational profiling, we performed
   comparative analysis for the functions in both the cases. PI3K-Akt
   signaling, cAMP signaling, focal adhesion, calcium signaling, oxytocin
   signaling, apelin signaling, ECM-receptor interaction, and thyroid
   hormone signaling are those pathways which are commonly enriched for
   both the cases differential expression and mutation in HNSCC which
   gives higher significance to these pathways for the selected disease
   while there are specific pathways for DEGs and mutated genes lists
   which means there are pathways which may be altered only because of
   overexpression of the genes or higher mutations rate. Mutation-specific
   altered pathways are circadian entrainment, cell cycle, phospholipase-d
   signaling, longevity regulating pathway, melanogenesis, FoxO signaling,
   estrogen signaling, long-term potentiation, and ubiquitin mediated
   proteolysis signaling pathways while insulin signaling, CAMs, TGF, TNF,
   tight junction, cGMP-PKG, phagosome signaling pathways are altered gene
   expression specific. From OSCC microRNAs analysis, 46 miRNAs appear
   overexpressed and the most affected functions are thyroid cancer,
   pathways in cancer, pancreatic cancer, Foxo signaling, chronic myeloid
   leukemia, HIF1 signaling and more and most of these pathways appears to
   strongly associated with head and neck cancer and oral cancer. Based on
   the clinical relevance of the hypoxic genes, there are a large number
   of genes which are highly significant and STC2, PGK1, P4HA1, HK1, SPIB,
   ANXA5, SERPINE1, HGF, PFKM, TGFB1, L1CAM, ELK4, EHF, and CDK2 are
   highly significant which have p-values even lower than or equal to
   0.0013. Similar to the expression and mutational profiling, the
   inferred pathways of the top-ranked genes are direct components of
   those pathways which directly affect the tumor initiation, propagation,
   and growth. Moreover, the RNA and miRNA expression analysis shows that
   there are common functions in the RNA and the miRNA pathways list such
   as HIF-1 signaling, Ras, MAPK, immune system associated pathways, and
   directly cancer-associated pathways.

Conclusions

   As mentioned that the major goal of this study was to understand the
   role of expression profiling of genes and the microRNAs and the
   mutational profiling of the genes and also the clinical relevance in
   case of HNSCC and based on results and the analysis, it leads to the
   conclusion that the critical pathways which could be dominantly
   enriched or altered in case of HNSCC are associated with metabolism,
   cell cycle, immune system, and hypoxia and the three different datasets
   of gene expression and microRNA expression, and the mutational data
   also leads to the conclusion that the pathways and the pathways
   components potentially associated with HNSCC and its progression.

Supplementary Information

   [123]Supplementary Information 1.^ (5.1MB, xlsx)
   [124]Supplementary Information 2.^ (288.6KB, pdf)

Acknowledgements