The Gene Expression Omnibus database (GEO) is a public functional database for high- throughput screening of gene expression data, microarray data, and gene chips. In this study, we recovered genome expression datasets from GEO (Affymetrix Human Genome U133A Plus 2.0 Array, Affymetrix Human Genome U95A Array, and Affymetrix Human Genome U133 Plus 2.0 Array) [GSE1987, GSE17073, GSE 54495, GSE118370]. The GSE1987, GSE17073, GSE 54495, and GSE118370 datasets contain 28, two, 17, and six NSCLC tissue samples and 9, 10, 13, and six non-cancer tissue samples, respectively ( Table 1 ). The software tools used in this study are listed in Supplementary Table 1 .

The differentially expressed genes were analyzed using GEO2R (http://www.ncbi.nlm.gov/geo2r), an interactive network tool for screening between NSCLC and non-cancer samples. The tool helps to compare and analyze two or more datasets under experimental conditions. The significant DEGs were detected using the adjusted p-values and the Benjamini- Hochberg false discovery rate. The probe sets with no gene names or genes with multiple probe sets were removed from the study. The upregulated genes were identified using cut- off values of P<0.05 and logFC >1, while the downregulated genes were identified with logFC <1.

The c uration of significant DEGs was performed with the help of the Comparative Toxicogenomics Database (CTD), the Online Mendelian Inheritance in Man (OMIM), PubMed, and MeSH databases in order to curate their role in NSCLC. The DAVID database was used to retrieve the gene symbol, name, and UniProt ID of the identified DEGs.

To identify the functional annotation and gene ontology of the shortlisted DEGs, the DAVID tool was used for functional enrichment analysis for their p- and FDR values. The biological function, transcription factors, and clinical phenotypes of NSCLC-related DEGs were analyzed using the FunRich tool.

The genotype-phenotype association helps in decoding the genetic variations, which aids in understanding the mutations arising from cancer and the inherited disease- related processes. Several single nucleotide variants (SNVs) are part of the human genome and are involved in the progression of the disease. The missense SNVs present at the PTM (post-translation modification) protein sites are associated with disease progression due to the substitution of approximately 21% of amino acids. This chemical modification of the amino acid ultimately alters the function of the protein (17). The online web tool ActiveDriverDB database was used to analyze the mutations associated with differentially expressed genes (17). The overall visual summary of the position, frequency, and functional significance of the mutations in the DEGs was provided by the needle plot mutation analysis. We observed the predicted disordered region and the PTM sites for all mutations in the protein sequence. The position of the gene sequence was observed by placing the pin along the protein. The figure legends explain the related mutational effects and PTM sites.

The protein-protein network analysis was performed to study the interaction of each protein with one or more genes associated with its molecular functions (18). The network reveals the altered activity of these genes in normal or pathological conditions. The potential NSCLC- related gene signatures associated with other genes whose dysfunction results in the disease state were identified through this network. The STRING database was used to analyze the protein-protein interactions of the cDNA dataset DEGs (19). The high confidence score interaction was used in this network analysis, having a score of 0.9-1. The target genes identified by this network were then studied for their role in NSCLC using Cancer GeneticsWeb, the National Cancer Institute, and the OMIM databases. Cytoscape (version 3.9.1) (20) was used to visualize the molecular network and Network Analyzer was used to calculate the topological properties of the networks. The nodes in the network categorized the degree of annotation between genes and diseases.

T o validate the expression pattern of identified DEGs in normal and tumor tissues, the GEPIA database was utilized (21). GEPIA (Gene Expression Profiling Interactive Analysis) is helpful in differential expression analysis, patient survival analysis, correlation and profiling plotting, and several other key interactive and customizable features. G ene expression plots in GEPIA are based on TCGA clinical annotations. The Human Protein Atlas database was used to validate the translational levels of the identified oncogenes. V alidation was performed using immunohistopathologic sections between the normal and OSCC samples (https://www.proteinatlas.org/).

The comparative toxicogenomic database (CTD) was used to perform the toxicogenomic analysis. The CTD helps in retrieving the exposome data, which explores the chemical-genome to phenome association. The analysis investigates the mechanism of functional pathway signaling toward the progression of the disease. The chemical-gene and disease interactions are inferred, revealing the particular expression of the gene and its association with disease (22).

microRNAs (miRNAs) are small non- coding RNAs that act as post-translational regulators that influence the genes involved in biological signaling pathways. In order to study the gene etiology, the expression and functional role of miRNA play a significant role (23). The disease- specific functional and molecular annotation of DEGs can be investigated by predicting the miRNA target (24). The miRDB online database was used to predict the NSCLC- related DEG miRNA targets for functional target prediction. The prediction data includes several important descriptions, such as 3′-UTR region length, miRNA-candidate target pairs along with target prediction scores, 3′-UTR sequences, miRNA seed binding sites, miRNA target sequences, etc. A score >80 was considered reliable, and the miRNA target was ranked (25, 26). The miRNA targets were analyzed for their site of expression and the biological pathways in which they are commonly used, with the help of the FunRich tool.

The drug -gene network analysis was performed to correlate the shortlisted DEGs with the chemicals/drugs that affect the activity of these genes. The CTD database was used to screen the chemical and disease relationships with default parameters. A direct link between the DEGs and anticancer drugs was developed. The FDA approval status of the identified drugs was verified with the Drug Bank database.

The four datasets GSE17073, GSE1987, GSE54495, and GSE118370 were used to screen and identify the DEGs after obtaining the microarray normalization results. The dataset GSE17073 contained 30 significant DEGs, GSE1987 contained 508 significant DEGs, GSE54459 contained five significant DEGs, and GSE118370 contained 707 significant DEGs. The volcano plot and the mean difference plot of DEGs from different datasets revealed the significant genes that were upregulated and downregulated. The blue color shows the down regulated genes, while the red indicates the upregulated genes ( Figure 1 ). In GSE1987, GSE17073, GSE54495, and GSE118370, 403, zero, one, and 354 genes are upregulated, respectively, while 405, 30, four, and 353 DEGs are downregulated. In the Venn diagram, the overlap between the datasets was observed, which showed an overlap of 10 significant DEGs with a p-value of <0.05 and logFC <1 for downregulated genes and logFC >1 for upregulated genes ( Figure 2 ). Table 2 shows the expression profiles of the microarray datasets.

From four datasets, we found 10 common DEGs whose gene symbols and biological annotations were retrieved using the David tool. D isease-gene curation was performed by text mining with the help of PubMed, PMC, PMIM, and MeSH ( Supplementary Table 1 ). It was observed that the DEGs DOCK4, ID2, GATA3, SASH1, GJA4, TBX2, HBEF, NPR1, CD24, and DDR1 were the most curated terms in the databases. The mapping of these genes by cancer genetics was performed to further analyze their role in carcinogenesis ( Figure 3 ).

The enrichment analysis of DEGs showed that these genes were significantly linked to the regulation of epithelial cell differentiation, ear development, mammary gland alveolus development, regulation of cardiac contraction, cellular senescence, negative regulation of transcription, DNA-templated synthesis, positive regulation of endothelial cell migration, positive regulation of smooth muscle cell proliferation, response to estrogen, and cell chemotaxis ( Table 3 ). The clinical phenotypes associated with the GATA3 gene were n ephrosis and renal agenesis. Vaginal agenesis, septate vagina, uterine agenesis, and uterus didelphys are rarely associated with this gene ( Figure 4 ). The transcriptome analysis revealed expressive transcription factors encoded by these genes, such as LMO2, ELF2, TBX5, PPARG, and FOXC1 ( Figure 4 ).

Dock4 consists of 30 post-translational modification (PTM) sites having 305 recurrent cancer mutations on the 7- chromosome number negative strand encoding 1966 protein residues with 22.69% predicted disordered regions. The mutation visualization plot shows the DOCK4 isoform mutation at position 378 with the reference amino acid residue V compared to the mutated amino acid residue L, which is enriched for phosphorylation-type mutations, showing the distal mutation PTM impact with the affected site. At position 1770, the L amino acid residue that was enriched with the phosphorylation-type mutation showed proximal mutation. The PTM impact with affected sites were 1769S. SASH1 showed a 68.48% predicted disordered region with 239 mutations observed on chromosome number 6 on the positive strand. 247 protein residues were encoded and 63 PTM sites were observed. At positions 370 and 421, the reference G and R amino acid residues were compared with W and C mutated amino acid residues showing distal PTM mutation impact with affected sites 374S and 417S, respectively. The affected sites at positions 925, 1231, and 1008 have P, L, and S residue sites enriched with a phosphorylation-type mutation affecting PTMs showing a network-rewiring motif gain mutation with mutated amino acid residue R, respectively. A proximal PTM impact at the 837S and 839S affected sites at position 841 was observed with reference amino acid residue D compared to the mutated Y amino acid residue. Similarly, the mutational analysis of GJA4 showed 48.95% of the sequence predicted for disease pathophysiology. A total of 64 mutations were found on the positive strand of chromosome number 1 for GJA4. The number of PTM sites was one, with 333 amino acid residues. In protein HBEGF, 21 mutations were observed on the negative strand of chromosome number 5 encoding 208 protein residues with 47.12% predicted disordered regions ( Figure 5 ).

The mutational analysis of NPR1 showed 15 PTM-affected sites with 151 mutations on the positive strand of chromosome number 1 encoding 1061 protein residues, representing 5.84% of predicted disordered regions. GATA3 showed 76.98% disordered regions with 104 mutations on the positive strand of chromosome number 10 encoding 443 protein residues. In position 157, reference amino acid P showed network rewiring motif loss PTM impacts at affected sites 156S and 162S compared with mutated amino acid T. In positions 278, 268, and 200, reference amino acids G, A, and H and mutated amino acids S, E, and Q showed distal PTM impact enriched with phosphorylation, methylation, and ubiquitination-type mutations, respectively. Similarly, the number of PTM sites for DDR1 was 20, with 913 protein residues having 22.34% of the predicted sequence for the disordered region. A total of 155 mutations were found on the positive strand of chromosome number 6 for DDR1 ( Figure 5 ). Eight isoforms of DDR1 were found ( Supplementary Table 2 ).

The STRING database was used to retrieve the related nodes and edges of all NSCLC- associated DEGs ( Figure 6 ). The PPI enrichment p-value was 0.012 with 15 nodes and 17 edges. Figure 6B shows the upregulated and downregulated oncogenes and their associated genes, allowing us to evaluate their biological functions.

The Cancer Genome Atlas (TCGA) and the GEPIA databases were employed to further validate our findings. The GEPIA NSCLC data showed that the expressions of these identified DEGs were significantly different between the normal and tumor tissues ( Figure 7 ). The trend was the same as observed in our data, which is consistent with the GEO analysis. Moreover, the Human Protein Atlas database, which showed deregulation of the expression of these seven genes, was used to obtain their immunohistochemistry staining data ( Figure 8 ). E xpression of the hub genes HBEGF, GJA4, DDR1, CD24, TBX2, GATA3, and SASH1 did not appear to be prognostic in lung cancer. No pathological data were found for DOCK4, NPR1, or ID2 genes.

The chemical-genotype-phenotype exposome data that may lead to disease progression were explored with the help of toxicogenomic analysis. The effect of different environmental chemicals on the activity and expression of NSCLC- associated DEGs was curated. The data revealed the effects of several chemicals on the increase or decrease of the expression of these NSCLC DEGs on gene activity at different levels ( Figure 9 ). It was also revealed that the co-treatment expression leading to disease occurrence and the same chemical exposure showed different reactivity for different genes. For example, arsenic trioxide decreases the expression of GJA4, but it affects the binding of GATA3. The details of the effect of these chemicals on NSCLC genes are shown in Table 4 .

The miRNA targets were predicted with the help of the miRDB database based on the algorithms. The miRNA targets such as hsa-miR-302c-5p, hsa-miR-4531, hsa-miR-11181-5p, hsa-miR-4476, hsa-miR-338-5p, hsa-miR-194-5p, hsa-miR-4279, hsa-miR-4742-3p, hsa-miR-4530, and hsa-miR-199a-5p were predicted for DOCK4, SASH1, ID2, GJA4, TBX2, HBEGF, NPR1, GATA3, and CD24, respectively. The onset and progression of the disease may be caused by the deregulation of these genes. Table 5 shows the predicted scores, the total number of miRNA hits, and the seed location of the DEGs. The functional enrichment analysis of these miRNA targets was performed using the FunRich tool. The analysis revealed some important biological pathways associated with these targets, such as the PDGF receptor signaling network, the ErbB receptor signaling network, the glypican signaling pathway, the TRAIL signaling pathway, and the plasma membrane estrogen receptor signaling pathway ( Figure 10 ). The sites of expression for these miRNA targets analyzed were brain (51.1%), placenta (68.6%), ovary (56.2%), kidney (65.1%), heart (38.5%), skeletal muscles (58.5%), and lung (62.9%), with p<0.001.

In order to explore the available treatment options for NSCLC, a toxicogenomic analysis was performed. The anti-cancer drugs gefitinib, doxorubicin, lapatinib, cisplatin, mitomycin, cyclophosphamide, and thalidomide showed interactions with the hub genes using the CTD database. Several other FDA- approved drugs also showed interactions with the seed genes, such as theophylline, rosiglitazone, aspirin, estradiol, phenylephrine, atorvastatin, and valproic acid. These drugs may serve as novel targets for these genes in the treatment of NSCLC ( Figure 9 ).