Unraveling the therapeutic potential of Ocimum sanctum phytoconstituents in the mechanism of oral squamous cell carcinoma: An in silico approach

2.1. Identification of phytochemicals in OS

The phytochemicals found in OS were discerned through an extensive review of scholarly literature, specifically published research articles that elaborate on the chemical constituents of the plant [19].

2.2. Analysis of the structural composition of phytoconstituents in OS

The molecular weight, molecular formula, and PubChem Compound Identifications (CIDs) of the phytoconstituents derived from the extract of OS were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Compounds that were not present in the PubChem database were omitted from subsequent analysis.

2.3. Identification of target genes associated with the phytoconstituents of OS

The PubChem CIDs of phytoconstituents were submitted to the Bioinformatics Annotation Tool for Molecular Mechanism of TCM (BATMAN-TCM) database (http://bionet.ncpsb.org.cn/batman-tcm/) for the purpose of target protein (gene) analysis. BATMAN-TCM is an online platform designed to catalogue both established and anticipated interactions between various compounds and their corresponding target proteins.

2.4. Genes implicated in OSCC

GeneCards (https://www.genecards.org/) serves as an extensive repository that systematically compiles and organizes data regarding gene-disease associations from a multitude of sources. By organizing the data systematically, it offers significant insights into the genetic underpinnings of human diseases. This database was employed to discern genes linked to OSCC.

2.5. Network development

The networks were constructed utilizing Cytoscape 3.10.3 (https://cytoscape.org/), an open-source software platform developed in Java. Cytoscape serves as a prominent tool for the visualization of intricate networks, facilitating the integration of diverse forms of attribute data. The networks are composed of nodes, which symbolize compounds and/or genes, edges that denote the relationships between these compounds and genes, and communication points that link various interconnected elements, including sources, bioactives, and targets. In the examination of network and biological interactions, a variety of Cytoscape tools were utilized. Sub-networks were developed to investigate particular relationships, including Source-Bioactive and Bioactive-Target networks. Utilizing the bioactive-target network, we have identified potential compounds that may exhibit therapeutic efficacy against OSCC. The prominent active constituents and prevalent genes of OS in relation to OSCC were discerned through network analysis.

2.6. Anticipating the toxicological profiles of phytoconstituents

ProTox-3 represents a sophisticated web-based instrument engineered to forecast the toxicity of chemical compounds, thereby facilitating the processes of drug discovery and development. Available through the ProTox-3 Web Server (https://tox.charite.de/protox3/index.php?site=compound_input), this platform utilizes a similarity-based methodology, examining compounds with established median lethal doses (LD50) and discerning toxic structural fragments. ProTox-3 offers an extensive evaluation of toxicity, encompassing acute toxicity, hepatotoxicity, carcinogenicity, immunotoxicity, mutagenicity, cytotoxicity, and LD50 values (articulated in mg/kg body weight). Furthermore, it categorizes compounds into distinct toxicity classes (I–VI) according to their toxic dose thresholds. Employing the degree layout within Cytoscape, three principal bioactive compounds—Quercetin, Apigenin, and Ursolic Acid—were discerned for the purpose of toxicity analysis. Their corresponding SMILES representations are as follows: C1=CC(=C(C=C1C2=C(C(=O)C3=C(C=C(C=C3O2)O)O)O)O)O for Quercetin, C1=CC(=CC=C1C2=CC(=O)C3=C(C=C(C=C3O2)O)O)O for Apigenin, and C[C@@H]. 1CC[C@@] 2(CC[C@@]) 3(C(=CC[C@H] 4[C@] 3(CC[C@H] 5[C@@] 4(CCC@@HO)C)C)[C@@H] 2[C@H] 1C)C)C(=O)O (Ursolic Acid) was input into the ProTox-3 platform for the purpose of predicting toxicity. The examination yielded a comprehensive toxicity profile for each compound, assessing their potential hazards across various toxicological parameters.

2.7. Analysis of protein–protein interactions (PPIs)

PPIs are fundamental to biological processes (BPs) and are pivotal in comprehending the intricate systems that operate within living cells. The assembly of target genes underwent examination through the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) to develop the PPI network. The ultimate collection of shared genes was presented to STRING, leading to the establishment of a PPI network along with two distinct clusters. The resultant data were subsequently visualized and subjected to analysis utilizing Cytoscape 3.10.3. This platform synthesizes experimental data, computational forecasts, and meticulously curated knowledge to investigate molecular interactions and the functional relationships that exist between proteins. In order to maintain a high level of reliability, the analysis was confined to Homo sapiens.

2.8. Identification of central genes

The PPI network was subsequently constructed utilizing Cytoscape 3.10.3, a widely recognized bioinformatics application for the visualization and integration of biological data. The Cytoscape plugin cytoHubba was employed to examine clusters and discern highly interconnected areas within the PPI network. The proteins exhibiting the highest maximum neighborhood component (MNC) rankings were identified as the Hub Genes.

2.9. Gene ontology (GO) and pathway enrichment analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/, ver. 6.8) and the ShinyGO (ver.0.77) database (http://bioinformatics.sdstate.edu/) were employed for the analysis of GO and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment, respectively. The prevalent genes discerned from the Cytoscape database were subsequently imported into the DAVID database, wherein KEGG pathway enrichment analysis and GO annotation were performed utilizing the “Homo sapiens” configuration. DAVID serves as a bioinformatics instrument that aids in the annotation and interpretation of gene lists, whereas ShinyGO is predominantly centered on pathway enrichment analysis. KEGG constitutes an extensive repository of pathways, presenting visual depictions of BP, while GO functions as a repository for functional genomics, supplying definitions and classifications pertinent to gene functions. Pathways pertinent to oral cancer were meticulously chosen in accordance with clinical and pathological data. The ShinyGO database was employed to visualize the common pathways identified across various analyses. In order to conduct a more thorough analysis and effectively present the data, a bubble chart was created utilizing the SRplot database (http://www.bioinformatics.com.cn/en), an online platform dedicated to bioinformatics that excels in data processing and visualization. For DAVID, the Benjamini correction for multiple testing, the Functional Annotation Chart, and Homo sapiens as a background were utilized. The Benjamini–Hochberg method of false discovery rate (FDR) adjustment was used for ShinyGO. Adjusted p-value FDR < 0.05 was used to filter enriched KEGG pathways and GO keywords.

2.10. Evaluation of molecular docking interactions between central genes and compounds

Through an analysis of the degree layout in Cytoscape, we discerned three principal compounds exhibiting the highest degree: quercetin, apigenin, and ursolic acid. The application of molecular docking facilitated a more profound exploration of the interactions between the candidate proteins (hub genes) and these three compounds. The docking simulation was meticulously crafted to evaluate the interactions between the hub targets and the compounds involved. CB-Dock (https://cadd.labshare.cn/cb-dock2/php/blinddock.php) was employed for the molecular docking simulations, demonstrating its ability to autonomously identify active sites within a specified protein, estimate their centres and dimensions, and adjust the grid box size in accordance with the query ligands. The crystal structure of the target protein was sourced from the Protein Data Bank (PDB) (https://www.rcsb.org/). In a similar way, the three-dimensional architecture of the compounds was ascertained through the utilization of the PubChem compound database (https://pubchem.ncbi.nlm.nih.gov/). The protein and ligand structures served as inputs for CB-Dock, which conducted docking studies to assess the binding activities of the compounds and proteins. The Discovery Studio Visualizer software (https://discover.3ds.com/discovery-studio-visualizer-download) served as the tool for the visualization and analysis of the docked results.

2.11. Validation of protein structure through the Ramachandran plot

The Ramachandran plot is a fundamental tool for assessing protein structure quality by analyzing the backbone dihedral angles φ (phi) and ψ (psi). It distinguishes sterically favorable from unfavorable conformations, highlighting regions that correspond to stable secondary structures such as α-helices, β-sheets, and turns. Common in structural biology, protein modeling, and molecular dynamics, it serves as a key validation metric: well-folded proteins show most residues in allowed regions, while those in disallowed zones may indicate steric clashes, errors, or flexible functional sites. The present analysis and visualization were performed using Discovery Studio Visualizer.

2.12. Analysis of gene expression levels and overall survival concerning central genes

UALCAN, the University of Alabama at Birmingham Cancer Data Analysis Portal (https://ualcan.path.uab.edu/), is an extensive online resource tailored for the examination of cancer OMICS data. This work synthesizes publicly accessible datasets, such as The Cancer Genome Atlas (TCGA) and the Clinical Proteomic Tumor Analysis Consortium (CPTAC). This investigation employed TCGA data to analyze the expression of particular genes linked to OSCC. UALCAN provides a platform for researchers to investigate gene and protein expression, promoter methylation patterns, and survival correlations in OSCC, taking into account a range of clinical parameters including tumor stage, grade, race, and HPV status. Moreover, the platform enables comprehensive analysis across various cancer types, fosters the identification of biomarkers, and supports investigations into epigenetic regulation, thereby offering significant insights into the molecular changes associated with cancer. Moreover, UALCAN is interconnected with external databases, including GeneCards and DRUGBANK, thereby augmenting its value in the realm of molecular oncology research. Through the provision of an accessible interface and sophisticated bioinformatics functionalities, UALCAN emerges as an essential resource for elucidating the molecular mechanisms that underpin pathology. The platform facilitates the discernment of potential prognostic markers and therapeutic targets, thereby advancing the formulation of personalized treatment strategies for OSCC. The Benjamini–Hochberg method of FDR adjustment was used for Gene expression and methylation. TCGA—Head and Neck Squamous cell carcinoma (HNSCC) cohorts include cancers from multiple anatomical sites (oral cavity, oropharynx, larynx, and hypopharynx), and the use of HNSCC data without site-specific explanation may introduce this heterogeneity. Our TCGA/UALCAN analyses are based on the broader HNSC population, and we explicitly note that non-oral locations might have played a role in the findings.

3.1. Possible active compounds and bioactive substances derived from OS

A thorough examination of existing scholarly articles revealed a total of 82 bioactive substances present in OS. The analysis of these compounds was conducted through the PubChem database, allowing for the retrieval of their molecular weight, molecular formula, and PubChem CIDs. Compounds that were not present in the PubChem database were omitted from further analysis, leading to a final selection of 66 compounds for subsequent investigations.

3.2. Target identification of OS related bioactives and linked to OSSC

The BATMAN-TCM database was utilized to discern bioactive compounds associated with OS. Out of the 66 compounds examined, 55 were identified as possessing both established and anticipated target proteins. The screening resulted in the identification of 1,450 known and predicted target genes linked to these 55 compounds. After excluding the anticipated target genes, a total of 885 unique known target genes were preserved for subsequent analysis. The GeneCards database has catalogued a total of 6,409 genes associated with OSCC. Through the meticulous cross-referencing of the 55 bioactive-related genes of OS with the genes associated with OSCC, a total of 435 common genes were carefully curated for subsequent analysis (Fig. 1).

Figure 1. A Venn diagram illustrating the intersection of genes associated with OSCC (OSCC, yellow) and OS (blue). The examination uncovered 5,974 genes (92.4%) that are exclusive to OSCC, 435 genes (6.7%) that are shared between both entities, and 59 genes (0.9%) that are distinctive to OS. [Click here to view]

3.3. Principal compounds/bioactive elements, along with target genes of OS in relation to OSCC

Upon the establishment of the compound-gene interaction network (Fig. 2), it became evident that 12 compounds derived from OS demonstrated the most significant degree of association as determined by the degree layout analysis conducted in Cytoscape. The network illustrated the engagement of these 12 compounds with 89 shared target genes (Fig. 3).

Figure 2. The network illustrates the correlation of bioactive compounds (Green) found in OS (Pink) alongside the genes associated with OSCC (Blue). [Click here to view]

Figure 3. The network demonstrates the correlation of 12 bioactive compounds (Green) present in OS (Pink) with 89 associated genes related to OSCC (Blue). [Click here to view]

3.4. Anticipated toxicity of the bioactive compounds

The predicted toxicity levels for Quercetin, Apigenin, and Ursolic Acid—reflected in their respective LD50 values of 159, 2,500, and 2,000 mg/kg, and designated as toxicity classes 3, 5, and 4—highlight differences in acute toxicity, with Quercetin exhibiting the highest. Notably, assessments suggest possible risks for neurotoxicity, respiratory toxicity, and carcinogenicity (Table 1), which emphasizes the importance of tailored risk assessment and ongoing monitoring to ensure therapeutic benefit while limiting adverse effects during the use of these phytochemicals.

Table 1. Anticipated toxicity assessment of quercetin, apigenin, and ursolic acid utilizing protox-3.0, grounded in molecular similarity, pharmacophoric features, fragment propensities, and advanced machine-learning algorithms.

3.5. Analysis of the PPI network and identification of key targets

The STRING database was employed to examine interactions among the target genes, facilitating the development of a PPI network. Figure 4 illustrates the interaction network of 89 common target genes, providing significant insights into their potential collective impact on cellular functions and signaling pathways. The network analysis revealed two distinct clusters, indicating that the proteins within each cluster are functionally interconnected and probably engaged in common BP, including metabolic regulation, signal transduction, or pathways associated with disease.

Figure 4. Network of PPIs involving 89 genes: Clustering relies on the interrelations among proteins. Cluster 1 comprises 75 proteins, represented as nodes, along with their interconnections, totalling 1,753 edges. PPI enrichment p-value: less than 1.0e -16. Cluster 2 comprises 11 proteins, represented as nodes, and their interconnections, totaling 38 edges. PPI enrichment p-value: less than 1.0e-16 [Click here to view]

3.6. Analysis of principal hub genes

The analysis of PPI networks was performed utilizing Cytoscape, incorporating the cytoHubba plugin to discern hub genes—proteins characterized by their extensive interaction degrees and pivotal positions within the network. The proposed hub genes are believed to be pivotal in facilitating the therapeutic effects against OSCC. The identification of the top ten hub genes was achieved through the application of various algorithms, with the MNC method employed to ascertain the key hub genes: IL6, AKT1, TNF, TP53, IL1B, CASP3, PTGS2, BCL2, JUN, and MAPK3 (Table 2; Fig. 5).

Figure 5. The ten principal hub genes discerned from the 89 shared genes between OSCC and OS. [Click here to view]

Table 2. A compilation of diverse algorithms employed to identify the top 10 hub genes utilizing CytoHubba plugin.

3.7. GO enrichment analysis

The 89 common genes underwent a comprehensive analysis to elucidate their functional roles and associated pathways, employing the DAVID bioinformatics tool, which synthesizes data from various sources. The O enrichment analysis of the prevalent genes revealed around 30 distinct GO terms. The examination indicated that these targets play a role in the enhancement of gene expression, the reaction to foreign stimuli, signal transduction, and the apoptotic process (BP). The analysis of cellular components (CCs) revealed that these targets are situated within the cytoplasm, nucleus, mitochondrion, and extracellular space. Regarding molecular function (MF), the targets predominantly engage in enzyme binding, ubiquitin-protein ligase binding, protein homodimerization activity, protein kinase activity, and DNA-binding transcription factor activity (Fig. 6).

Figure 6. The GO enrichment analysis, highlighting 10 BP, 10 MFs, and 10 CC sourced from the DAVID Database, pertaining to the targets in O. sanctum for the treatment of OSCC. [Click here to view]

3.8. KEGG enrichment analysis

The analysis of KEGG pathway enrichment revealed a total of 176 signaling pathways linked to the 89 shared genes. The foremost 45 pathways were chosen for additional examination. The examination disclosed that the foremost ten hub genes were intricately connected to various pathways, with significant correlations noted in the PI3K-Akt signaling pathway, cancer-related pathways, and the IL-17 signaling pathway (Fig. 7). Within this context, BCL2, AKT1, and TP53 were persistently engaged in the foremost ten pathways (Table 3). These genes were markedly associated with the PI3K-Akt signaling pathway, which is fundamental in governing a range of cellular functions related to carcinogenesis, encompassing metabolic processes, cell survival, proliferation, gene expression, and protein synthesis. The PI3K-AKT signaling pathway serves as a vital regulator of cellular processes such as survival, proliferation, metabolism, and apoptosis, frequently exhibiting dysregulation in various cancers, including OSCC. The disruption of this pathway results in unrestrained tumor proliferation, evasion of apoptosis, and increased metastatic potential. In this pathway, AKT1 serves as a pivotal entity facilitating the survival of tumor cells and their resistance to apoptosis. The heightened expression of BCL2, BCL-xL, and Mcl-1 amplifies anti-apoptotic signaling, whereas the inactivation of TP53 (p53), a pivotal tumor suppressor, often results in unchecked tumor advancement. CASP3, a key player in the execution of apoptosis, frequently exhibits downregulation, which plays a significant role in the development of resistance to programmed cell death. The modifications in epigenetics and the dysregulation of pivotal genes within the PI3K-AKT pathway position them as promising candidates for therapeutic targeting and prognostic evaluation in OSCC, underscoring the necessity for continued investigation into targeted treatment strategies (Fig. 8).

Figure 7. A bar plot illustrating the enrichment of the leading 45 signaling pathways linked to OSCC. The horizontal axis delineates the quantity of genes, whereas the vertical axis illustrates the various pathways. [Click here to view]

Table 3. Top 10 signaling pathways indicating important values and genes present in them.

Figure 8. PI3K-AKT signalling pathway: The PI3K-AKT signalling pathway represents a fundamental mechanism that is prevalent across nearly all signalling pathways. The red boxes denote the Hub genes that are implicated in the pathway. This pathway ultimately influences cellular proliferation primarily through mechanisms such as apoptosis and cell cycle arrest. [Click here to view]

3.9. Assessment of Hub target via molecular docking analysis

Ten hub genes were selected as target proteins for molecular docking studies to evaluate the validity of drug-target interactions. CB-DOCK was employed to investigate the docking potential of compounds with JUN, IL6, CASP3, BCL2, IL1B, AKT1, TP53, TNF, MAPK3, and PTGS2. A diminished energy value signifies that the ligand adopts a more stable conformation during its interaction with the receptor, suggesting an increased probability of contact. The binding energies observed between the compounds and the core target proteins were determined to be less than −5.0 in this study, indicating a noteworthy level of binding activity between the compounds and the core targets. The binding energies of particular interest are delineated in Table 4. The docking schematic representations delineating the interactions between the target proteins and compounds are presented in Tables 5 and 6.

Table 4. Molecular docking evaluations of quercetin, apigenin, ursolic acid, and central target proteins showing PDB IDs and binding energy.

The interactions of quercetin, apigenin, and ursolic acid with various hub gene-encoded proteins (JUN, IL6, CASP3, BCL2, IL1B, AKT1, TP53, TNF, MAPK3, and PTGS2) were examined through molecular docking, utilizing structures from the PDB to emphasize their three-dimensional conformations. The structural representation encompasses alpha-helices depicted in red, beta-sheets illustrated in cyan, coil regions in grey signifying flexibility, and loop or turn regions represented in green. It further elucidates unique binding patterns that could impact their therapeutic efficacy. The structural analyses of these protein-ligand complexes reveal that quercetin, through its robust interactions, notably stabilizes proteins such as IL6, AKT1, CASP3, and MAPK3, which play pivotal roles in the regulation of apoptosis and inflammation. Apigenin, through its interactions with analogous proteins, demonstrates binding capabilities that could influence protein conformations via hydrogen bonds and hydrophobic interactions, especially impacting proteins characterized by mixed secondary structure elements such as AKT1, TNF, and TP53. Ursolic acid, on the other hand, exhibits a distinct affinity for proteins that play crucial roles in immune response and tumor progression, effectively stabilizing AKT1, PTGS2, and TNF, characterized by predominant α-helices, as well as TNF and TP53, which are enriched in β-sheet domains. In examining the three compounds, it becomes evident that the regulatory proteins TP53 and JUN demonstrate a notable degree of structural adaptability, indicating that ligand binding may play a significant role in modulating their transcriptional and apoptotic functions. The comparative analysis reveals that although all three bioactive compounds address inflammation, apoptosis, and cancer progression pathways, quercetin demonstrates the most stabilizing interactions, apigenin affects protein dynamics via structural modulation, and ursolic acid shows a strong affinity for immune and tumor regulatory proteins, underscoring their unique yet interconnected therapeutic potential (Table 5).

Table 5. A three-dimensional visual representation illustrating the molecular docking of quercetin, apigenin, and ursolic acid compounds alongside the top ten hub genes (target protein).

Hub genes	PDB ID	Quercetin	Apigenin	Ursolic acid
JUN	5T01
IL6	1ALU
CASP3	6BDV
BCL2	2YV6
IL1B	9ILB
AKT1	6CCY
TP53	5ZCJ
TNF	2AZ5
MAPK3	2ZOQ
PTGS2	1CX2

The two-dimensional molecular docking interactions of quercetin, apigenin, and ursolic acid with target proteins showed unique but interrelated binding patterns, suggesting therapeutic potential. Hydrogen bonding (green dashed lines), hydrophobic interactions (yellow/orange dashed lines), and electrostatic forces are also shown. Through a synergistic interaction of hydrogen bonds, hydrophobic contacts, and electrostatic forces, quercetin stabilizes essential amino acid residues in active or allosteric locations. Apigenin forms several hydrogen bonds and varies in ligand orientation among proteins, showing that binding pocket structure affects its interactions. It stabilizes proteins via hydrophobic and electrostatic interactions, notably with aromatic and aliphatic residues, increasing its affinity for catalytic or allosteric locations. Ursolic acid, like quercetin, has a strong binding affinity with many interaction patterns. Electrostatic interactions promote stability, while hydrogen bonding and hydrophobic forces dominate stabilization. Hydrogen bonding dominates all three molecules, but hydrophobic and electrostatic contributions greatly impact their specificity and functional consequences, highlighting their medicinal potential (Table 6).

Table 6. Two-dimensional representations of molecular docking quercetin, apigenin, ursolic acid compounds, and the ten principal hub genes (target proteins)

Hub genes	PDB ID	Quercetin	Apigenin	Ursolic acid
JUN	5T01
IL6	1ALU
CASP3	6BDV
BCL2	2YV6
IL1B	9ILB
AKT1	6CCY
TP53	5ZCJ
TNF	2AZ5
MAPK3	2ZOQ
PTGS2	1CX2

3.10. Assessment of the structural integrity of the target protein through Ramachandran plot analysis

To understand crucial target proteins’ drug design suitability, the Ramachandran plot provides valuable insights into their structural integrity and flexibility. IL6, AKT1, TNF, TP53, and IL1B have a high number of residues in the favored and allowed regions, indicating structural stability. Proteins with residue aggregation in α-helical and β-sheet domains show minimal steric hindrance, making them excellent candidates for molecular docking and drug development. Their small modifications improve computational precision and experimental confirmation, boosting their medicinal potential. CASP3, PTGS2, and BCL2 are structurally stable, yet numerous residues are in disallowed areas. CASP3, a key apoptosis enzyme, has modest structural changes that may be linked to flexible loops that activate it. PTGS2, a key enzyme involved in inflammation, and BCL2, which regulates apoptosis, also adapt to their roles in complex cellular processes. Conversely, JUN and MAPK3 have more outliers, indicating structural flexibility. JUN, a transcription factor, and MAPK3, a signaling kinase, interact with cofactors, DNA, and regulatory molecules, requiring dynamic conformational changes. Residues in restricted locations imply flexible loops or disordered domains that permit essential functional interactions. Their biological flexibility makes structure modeling and medication targeting difficult. The Ramachandran plot shows that IL6, AKT1, TNF, TP53, and IL1B are stable, making them good therapeutic targets. CASP3, PTGS2, and BCL2 have minor stability fluctuations, indicating functional dynamics. JUN and MAPK3’s structural flexibility is essential for their regulatory activities in cellular pathways, requiring customized docking techniques. The findings illuminate protein stability, which informs molecular docking and treatment design based on protein structure (Fig. 9).

Figure 9. Ramchandran plot illustrating compounds and the top 10 hub genes (target protein) relation. [Click here to view]

3.11. Validation of expression for hub genes within the TCGA dataset utilizing the UALCAN database

The expression levels of genes in samples of head and neck squamous cell carcinoma (HNSC) derived from TCGA. A comparative analysis is conducted on the expression levels of various genes between normal tissue (n = 44), represented in blue, and primary tumor tissue (n = 520), depicted in red. The x-axis delineates the various sample types, contrasting normal with primary tumor, whereas the y-axis illustrates the levels of gene expression. A number of inflammatory and oncogenic genes (IL6 (p = 6.393800E-02), TNF (p = 4.15190000158105E-08), IL1B (p = 1.11022302462516E-16), PTGS2 (p = .67987845856032E-12)) exhibit upregulation in tumors, indicating a significant connection between inflammation and HNSC. The upregulation of pro-survival genes such as AKT1 (p < 1E-12) and JUN (p = 1.358600E-01) indicates a notable increase in oncogenic signaling. The tumor suppressor TP53 (p = 1.28850000002956E-06) exhibits heightened expression, presumably as a consequence of mutations that result in its impairment. The elevation of the apoptotic marker CASP3 (p = 4.30139999951784E-07) suggests a nuanced interplay between cellular survival and the mechanisms of cell death. The absence of notable alterations in BCL2 (p = 1.073800E-01) and MAPK3 (p = 3.83700000000653E-05) indicates that their functions may be contingent upon the specific context of HNSC progression (Fig. 10).

Figure 10. Hub genes expression in HNSC categorized by sample types. The mRNA expression levels of hub genes in TCGA and Genotype-Tissue Expression (GTEx) databases. Figure depicts mRNA expression levels of the top 10 hub genes with statistically significant values in normal verses primary tumor (Box-Plot). The UALCAN database was used for the presentation of this data. [Click here to view]

The box plot delineates the expression levels of transcripts per million across various tumor grades in samples from TCGA. The x-axis delineates the samples into five distinct categories: The representation of Normal (n = 44) is depicted in blue, while Grade 1 (n = 27) is illustrated in orange, denoting well-differentiated, low-grade tumors. Grade 2 (n = 71) is shown in green, indicating moderately differentiated, intermediate-grade tumors. Grade 3 (n = 81) is presented in brown, corresponding to poorly differentiated, high-grade tumors, and Grade 4 (n = 264) is highlighted in pink, signifying undifferentiated, high-grade tumors. p < 0.05 was considered a statistically significant value. The y-axis measures gene expression levels, facilitating a comparative analysis among various tumor grades. The expression levels of genes including AKT1 (Normal-vs-Stage2; p = 9.16011710927478E-12), TNF (Normal-vs-Stage2; p = 8.02560000000563E-05), TP53 (Normal-vs-Stage4; p = 3.141100E-03), IL1B (Normal-vs-Stage4; p = 1.28597132942332E-12), CASP3 (Normal-vs-Stage2; p = 5.67910000000005E-05), PTGS2 (Normal-vs-Stage4; p=3.39379968572473E-10), and BCL2 (Normal-vs-Stage1; p= 9.616800E-01) rise with advancing tumor grades, indicating their potential involvement in the progression and aggressiveness of tumors. Conversely, IL6 (Normal-vs-Stage4; p = 4.432400E-02) and MAPK3 (Normal-vs-Stage3; p = 7.51570000000479E-05) exhibit reduced expression in tumors, suggesting possible tumor-suppressive roles or downregulation as a component of immune evasion strategies. JUN exhibits a variable expression pattern, demonstrating elevated levels in tumors; however, it lacks a robust correlation with tumor grade (Normal-vs-Stage1; p = 7.555900E-02), indicating a more intricate regulatory function in tumor development (Fig. 11).

Figure 11. The mRNA expression levels of hub genes in TCGA and Genotype-Tissue Expression (GTEx) databases. Figure depicts mRNA expression levels of top 10 hub genes with statistically significant values in normal as well as stages of oral cancer (Box-Plot). The UALCAN database was used for the presentation of this data. [Click here to view]

3.12. Verification of DNA methylation in hub genes within the TCGA dataset

The expression levels of genes in samples of head and neck squamous cell carcinoma (HNSC) derived from TCGA. A comparative analysis is conducted on DNA methylation levels of various genes between normal tissue (n = 50), represented in blue, and primary tumor tissue (n = 528), depicted in red. The Beta value indicates the level of DNA methylation ranging from 0 (unmethylated) to 1 (fully methylated). Different beta value cut-off has been considered to indicate hyper-methylation (Beta value: 0.7–0.5) or hypo-methylation (Beta-value: 0.3–0.25). p < 0.05 was considered a statistically significant value. The genes IL6 (p < 1E-12), CASP3 (p = 2.513600E-02), BCL2 (p = 6.16620000000179E-05), TNF (p = 1.63014046705712E-12), IL1B (p = <1E-12), TP53 (p = 9.283800E-01), AKT1 (p = 9.585900E-02), JUN (p = 6.827400E-01), MAPK3 (p = 3.909000E-01) and PTGS2 (p = 4.463600E-02) demonstrate notable hypomethylation, indicating a likelihood of overexpression and their roles in inflammation, cell survival, resistance to apoptosis, and oncogenic signaling in head and neck squamous cell carcinoma (HNSC). Among these, IL6, TNF, IL1B, and TP53 exhibit significant associations with the progression of HNSC, with IL6 and TNF assuming pivotal roles in the context of tumor-associated inflammation. Considering the epigenetic modifications present, focusing on these alterations may represent a potentially effective therapeutic approach for the treatment of HNSC (Fig. 12).

Figure 12. Promoter methylation level of genes in HNSC based on sample type. The Beta value indicates level of DNA methylation ranging from 0 (unmethylated) to 1 (fully methylated). Different beta value cut-off has been considered to indicate hyper-methylation (Beta value: 0.7–0.5) or hypo-methylation (Beta-value: 0.3–0.25). p < 0.05 was considered as statistically significant value. [Click here to view]

A number of genes, such as IL6 (Normal-vs-Stage2 and Normal-vs-Stage4; p < 1E-12), IL1B (Normal-vs-Stage2, Normal-vs-Stage3 and Normal-vs-Stage4; p = <1E-12), TNF (Normal-vs-Stage4; p = 1.7599255386358E-12), TP53 (Normal-vs-Stage3; p = 6.780000E-02), JUN (Normal-vs-Stage4; p = 8.212000E-01), CASP3 (Normal-vs-Stage2; p = 1.775610E-03)and BCL2 (Normal-vs-Stage1; p = 6.72310000005005E-07), demonstrate hypomethylation in tumor samples, suggesting their involvement in the advancement of cancer. Conversely, genes such as MAPK3 (Normal-vs-Stage4; p = 9.482400E-02) exhibit downregulation in advanced tumor grades, indicating a possible role in tumor suppression. The alterations in methylation patterns of critical genes may provide significant biomarkers for tumor progression and prospective therapeutic targets in HNSC (Fig. 13).

Figure 13. Analysis of promoter methylation levels of genes in head and neck squamous cell carcinoma in relation to tumor grade. The Beta value indicates level of DNA methylation ranging from 0 (unmethylated) to 1 (fully methylated). Different beta value cut-off has been considered to indicate hyper-methylation (Beta value: 0.7–0.5) or hypo-methylation (Beta-value: 0.3–0.25). p < 0.05 was considered as statistically significant value. [Click here to view]

3.13. Prognostic outcomes in individuals exhibiting protein expression

Survival curves utilizing the Kaplan–Meier method were constructed to assess the relationship between gene expression levels and overall survival in individuals diagnosed with head and neck cancer (n = 500). The x-axis delineates the duration of survival in days, whereas the y-axis signifies the likelihood of survival. Patients were categorized into groups of high expression (red) and low expression (blue) according to their gene expression levels. The evaluation of statistical significance was conducted through log-rank tests, accompanied by the relevant p-values. The cut off value was set as percentiles with median expression. IL6 [p value:0.01; FDR: over 50%; HR = 1.42 (1.09–1.86)] and BCL2 [p value:0.4229; FDR: 100%; HR = 1.13 (0.84–1.52)] demonstrate statistically significant correlations with survival, underscoring their potential utility as prognostic indicators or therapeutic targets. Conversely, genes including AKT1 [p value: 0.0098; FDR: over 50%; HR = 1.42 (1.09–1.86)], TNF [p value: 0.0491; FDR: over 50%; HR = 0.75 (0.56–1.0)], TP53 [p value: 0.1323; FDR:100%; HR = 0.78 (0.56–1.08)], IL1B, CASP3 [p value:0.1092; FDR:100%; HR = 1.29 (0.94–1.78)], TP53 [p value: 0.1323; ?FDR:100%; HR = 0.78 (0.56–1.08)], PTGS2 [p value: 0.1318; FDR: 100%; HR = 0.81 (0.62–1.06)], JUN [p value: 0.1957; FDR: 100%; HR = 1.23 (0.9–1.67)], and MAPK3 [p value: 0.2848; FDR: 100%; HR = 0.86 (0.65–1.13)], do not demonstrate notable survival disparities, indicating that their expression alone may not function as robust prognostic markers in HNSC. Future research may explore the synergistic impacts of these genes alongside their epigenetic alterations to enhance our understanding of their contributions to tumor progression (Fig. 14).?

Figure 14. The impact of gene expression levels on the survival of patients with HNSC. Kaplan–Meier overall survival analyses of patients with HNSC based on expression of the ten hub genes. (p < 0.05 was considered as statistically significant value, Hazard Ratio (95% CI), cut off as Percentiles with median expression and FDR) (“http://kmplot.com/analysis/index.php?p=service&cancer=pancancer_rnaseq”) [Click here to view]