Abstract Background This study was performed to identify disease-related genes and analyze prognostic values in nonsmoking females with non-small cell lung carcinoma (NSCLC). Materials and methods Gene expression profile [35]GSE19804 was downloaded from the Gene Expression Omnibus (GEO) database and analyzed by using GEO2R. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes were used for the functional and pathway enrichment analysis. Then, the Search Tool for the Retrieval of Interacting Genes, Cytoscape, and Molecular Complex Detection were used to construct the protein–protein interaction (PPI) network and identify hub genes. Finally, the Kaplan–Meier plotter online tool was used for the overall survival analysis of hub genes. Results A cohort of 699 differentially expressed genes was screened, and they were mainly enriched in the terms of ECM–receptor interaction, focal adhesion, and cell adhesion molecules. A PPI network was constructed, and 15 hub genes were identified base on the subset of PPI network. Then, two significant modules were detected and several genes were found to be associated with the cell cycle pathway. Finally, nine hub genes’ (UBE2C, DLGAP5, TPX2, CCNB2, BIRC5, KIF20A, TOP2A, GNG11, and ANXA1) expressions were found to be associated with the prognosis of the patients. Conclusion Overall, we propose that the cell cycle pathway may play an important role in nonsmoking females with NSCLC and the nine hub genes may be further explored as potential targets for NSCLC diagnosis and treatment. Keywords: non-small cell lung carcinoma, NSCLC, nonsmoking females, GEO2R, prognostic biomarkers, Kaplan–Meier plotter Introduction Lung cancer is one of the leading causes of cancer mortality worldwide.[36]^1 Non-small-cell lung cancer (NSCLC), the main type of lung cancer, accounts for 80%–85% of all cases, and small cell lung cancer (SCLC) accounts for around 20% of all cases. Although smoking is the major risk factor for lung cancer,[37]^2 only 7% of female patients with lung cancer have a history of cigarette smoking in Taiwan.[38]^3 In addition, other factors, such as environmental exposure[39]^4^,[40]^5 or hereditary factors,[41]^6 have been reported in association with nonsmoking lung cancer patients. However, the molecular mechanisms of NSCLC among nonsmoking women are unclear, although some genes, such as EML4-ALK,[42]^7 EGFR,[43]^8 TP53,[44]^9 and PIK3CA,[45]^10 have been found to be associated with lung cancer in never smokers. Thus, it is important to explore the molecular mechanisms involved in NSCLC’s onset and progression among nonsmoking females and identify effective biomarkers. Recently, gene microarray and bioinformatics analysis were widely used to identify potential biomarkers of cancer. Interestingly, studies were performed to identify disease-related genes, which were associated with prognosis of breast cancer, by using integrated bioinformatics methods.[46]^11^,[47]^12 Similarly, some studies were performed to find important key genes in lung cancer.[48]^13^–[49]^15 But there are limited studies on nonsmoking females with NSCLC. In 2010, Lu et al[50]^3 obtained a panel of differentially expressed genes (DEGs) from nonsmoking female NSCLC patients and normal samples. Using the same data, the present study was performed to further screen the DEGs and predict their underlying function by functional and pathway enrichment analyses. Then, protein–protein interaction network (PPI) networks and modules of PPI network were performed to identify hub genes. More importantly, the prognostic values of the hub genes were further confirmed by survival analysis in nonsmoking females with NSCLC. Materials and methods Microarray data Microarray expression profiles of [51]GSE19804 were obtained from the Gene Expression Omnibus (GEO) database. [52]GSE19804, which was based on the platform of the [53]GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array, Santa Clara, CA, USA), consisted of 60 samples from nonsmoking female NSCLC patients and 60 normal samples.[54]^3 Identification of DEGs The DEGs between NSCLC samples and normal controls were identified using GEO2R ([55]https://www.ncbi.nlm.nih.gov/geo/geo2r/), which is a web tool that is applied to screen DEGs by comparing two groups of samples. |log FC| >1.5 and P<0.01 were selected as the cutoff criterion. Functional and pathway enrichment analysis Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs were performed using Database for Annotation, Visualization, and Integrated Discovery (DAVID),[56]^16 which is a comprehensive set of functional annotation tools. P<0.05 was set as the cutoff criterion. PPI network analysis of DEGs and modules selection The Search Tool for the Retrieval of Interacting Genes (STRING) database[57]^17 was used to construct a PPI network for DEGs. Then, the modules of PPI network with significant gene pairs (combined score >0.8) was obtained by the Molecular Complex Detection (MCODE) plugin[58]^18 in Cytoscape.[59]^19 MCODE scores >14 and the number of nodes >15 were used as the cutoff values. Moreover, the function and pathway enrichment analysis of genes in each module was performed using DAVID. Survival analysis of the hub genes Kaplan–Meier plotter is a web tool that predicts the prognostic values of genes in some cancer patients, including breast, ovarian, lung, and gastric cancer patients ([60]http://kmplot.com/analysis/index.php?p=background). The nonsmoking female patients with NSCLC were divided into two groups according to the particular gene expression level (high vs low expression).[61]^20 Based on these categories, overall survival (OS) analysis of the two patient groups was then compared by the tool. The HR with 95% CIs and log-rank P-value were calculated and showed. Results Identification of DEGs The total number of samples included 60 NSCLC samples and 60 normal samples. A total of 699 genes were identified, including 161 upregulated and 538 downregulated genes. GO and KEGG pathway enrichment analysis We used the DAVID for GO and KEGG pathway enrichment analysis based upon DEGs. The top five GO terms of DEGs are shown in [62]Table 1. As to biological process, the DEGs were significantly enriched in the regulation of cell proliferation, cell adhesion, biological adhesion, vasculature development, and blood vessel development. About cellular component, the DEGs were significantly enriched in the extracellular region part, extracellular region, extracellular space, cell surface, and proteinaceous extracellular matrix. In addition, the most enriched GO terms in molecular function were growth factor binding, carbohydrate binding, calcium ion binding, pattern binding, and polysaccharide binding. On the other hand, the most enriched KEGG pathway terms were as follows: ECM–receptor interaction, focal adhesion, cell adhesion molecules, neuroactive ligand-receptor interaction, and complement and coagulation cascades. Table 1. Functional and pathway enrichment analysis of DEGs in nonsmoking females with NSCLC ID Description Count P-value __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ GO_function GO-BP:0042127 Regulation of cell proliferation 55 2.29E-10 GO-BP:0007155 Cell adhesion 51 2.94E-10 GO-BP:0022610 Biological adhesion 51 3.04E-10 GO-BP:0001944 Vasculature development 28 1.03E-09 GO-BP:0001568 Blood vessel development 27 2.83E-09 GO-CC:0044421 Extracellular region part 87 5.32E-22 GO-CC:0005576 Extracellular region 130 2.84E-20 GO-CC:0005615 Extracellular space 61 8.88E-15 GO-CC:0009986 Cell surface 36 1.30E-10 GO-CC:0005578 Proteinaceous extracellular matrix 33 9.70E-10 GO-MF:0019838 Growth factor binding 15 6.47E-07 GO-MF:0030246 Carbohydrate binding 27 2.57E-06 GO-MF:0005509 Calcium ion binding 48 8.91E-06 GO-MF:0001871 Pattern binding 16 1.39E-05 GO-MF:0030247 Polysaccharide binding 16 1.39E-05 KEGG_ PATHWAY hsa04512 Ecm-receptor interaction 16 3.74E-08 hsa04510 Focal adhesion 16 0.001596 hsa04514 CAMs 11 0.008625 hsa04080 Neuroactive ligand-receptor interaction 16 0.014903 hsa04610 Complement and coagulation cascades 7 0.021259 [63]Open in a new tab Note: Top five terms were selected according to P-value. Abbreviations: BP, biological process; CAMs, cell adhesion molecules; CC, cellular component; DEGs, differentially expressed genes; ECM, extracellular matrix; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; NSCLC, non-small-cell lung cancer. Construction of PPI network and module identification To predict the interactions of the identified DEGs from the protein level, a PPI network from STRING was constructed. In total, 439 nodes and 1,056 edges were shown in the PPI network ([64]Figure 1). The subset of PPI network for the DEGs with a combined score >0.8 was performed to determine the hub genes. The top 15 genes were selected as hub genes, which were GNG11, UBE2C, CCNB2, NMU, ANXA1, FPR2, KIF20A, MELK, NUSAP1, DLGAP5, CENPF, BIRC5, TOP2A, TPX2, and HMMR ([65]Table 2). Subsequently, two modules with MCODE scores >14 and the number of nodes >15 were selected ([66]Figure 2), and enrichment analysis showed that the genes in the modules were mainly associated with cell cycle, chemokine signaling pathway, cytokine-cytokine receptor interaction, and neuroactive ligand-receptor interaction ([67]Table 3). Figure 1. [68]Figure 1 [69]Open in a new tab PPI network of DEGs. Abbreviatins: DEGs, differentially expressed genes; PPI, protein–protein interaction. Table 2. The hub genes that had a degree >20 in PPI network Gene Regulation Degree __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ GNG11 Down 29 UBE2C Up 28 CCNB2 Up 23 NMU Up 22 ANXA1 Down 22 FPR2 Down 22 KIF20A Up 22 MELK Up 21 NUSAP1 Up 21 DLGAP5 Up 21 CENPF Up 21 BIRC5 Up 21 TOP2A Up 21 TPX2 Up 21 HMMR Up 21 [70]Open in a new tab Abbreviation: PPI, protein–protein interaction. Figure 2. Figure 2 [71]Open in a new tab The two modules identified in the PPI network of the DEGs. Note: (A) Module 1 and (B) module 2. Blue represents upregulated DEGs; yellow represents downregulated DEGs. Abbreviatins: DEGs, differentially expressed genes; PPI, protein–protein interaction. Table 3. Functional and pathway enrichment analysis of the DEGs in modules ID Description Count P-value __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ Module 1 GO_BP GO:0000280 Nuclear division 13 4.29E-18 GO_BP GO:0007067 Mitosis 13 4.29E-18 GO_BP GO:0000087 M phase of mitotic cell cycle 13 5.34E-18 GO_BP GO:0000279 M phase 14 6.27E-18 GO_BP GO:0048285 Organelle fission 13 6.99E-18 GO_CC GO:0015630 Microtubule cytoskeleton 14 3.15E-14 GO_CC GO:0005819 Spindle 10 2.30E-13 GO_CC GO:0044430 Cytoskeletal part 14 3.57E-11 GO_CC GO:0005856 Cytoskeleton 15 1.91E-10 GO_CC GO:0043232 Intracellular non-membrane- bounded organelle 16 6.69E-08 GO_MF GO:0005524 ATP binding 9 2.60E-04 GO_MF GO:0032559 Adenyl ribonucleotide binding 9 2.85E-04 GO_MF GO:0030554 Adenyl nucleotide binding 9 4.10E-04 GO_MF GO:0001883 Purine nucleoside binding 9 4.56E-04 GO_MF GO:0001882 Nucleoside binding 9 4.78E-04 KEGG_PATHWAY hsa04110 Cell cycle 3 0.008428591 Module 2 GO_BP GO:0006935 Chemotaxis 8 4.59E-11 GO_BP GO:0007186 G-protein coupled receptor protein signaling pathway 12 8.22E-11 GO_BP GO:0007166 Cell surface receptor linked signal transduction 13 4.90E-10 GO_BP GO:0006954 Inflammatory response 6 8.52E-06 GO_BP GO:0006955 Immune response 7 2.17E-05 GO_CC GO:0005615 Extracellular space 5 0.003963123 GO_CC GO:0044421 Extracellular region part 5 0.013058837 GO_CC GO:0005576 Extracellular region 6 0.040928448 GO_MF GO:0008009 Chemokine activity 5 9.63E-08 GO_MF GO:0042379 Chemokine receptor binding 5 1.25E-07 GO_MF GO:0005125 Cytokine activity 5 3.17E-05 GO_MF GO:0001653 Peptide receptor activity 4 1.77E-04 GO_MF GO:0008528 Peptide receptor activity, G-protein coupled 4 1.77E-04 KEGG_PATHWAY hsa04062 Chemokine signaling pathway 7 9.05E-07 KEGG_PATHWAY hsa04060 Cytokine-cytokine receptor interaction 6 1.25E-04 KEGG_PATHWAY hsa04080 Neuroactive ligand-receptor interaction 5 0.001563563 [72]Open in a new tab Abbreviations: BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function. Survival analysis Kaplan–Meier plotter was used to predict the prognostic value of 15 identified hub genes. Our results showed that high expression of UBE2C was associated with worse OS for NSCLC patients, as well as DLGAP5, TPX2, CCNB2, BIRC5, KIF20A, and TOP2A (P<0.05) ([73]Figure 3). Additionally, low expression of ANXA1 was associated with poorer OS for NSCLC patients, as well as GNG11 (P<0.05) ([74]Figure 3). Figure 3. [75]Figure 3 [76]Open in a new tab Kaplan–Meier curves depicting OS in nonsmoking females with NSCLC with high and low expression of ANXA1 (A), GNG11 (B), BIRC5 (C), CCNB2 (D), DLGAP5 (E), KIF20A (F), TOP2A (G), TPX2 (H), and UBE2C (I). Abbreviations: NSCLC, non-small cell lung carcinoma; OS, overall survival. Discussion In this study, a total of 699 DEGs were screened, including 161 upregulated genes and 538 downregulated genes. Moreover, we selected two significant modules with several key DEGs in the regulatory network of nonsmoking female patients with NSCLC. Then, survival analysis of these genes determined that seven overexpressed genes (UBE2C DLGAP5, TPX2, CCNB2, BIRC5, KIF20A, and TOP2A) and two downregulated genes (GNG11 and ANXA1) were significantly correlated with worse OS of nonsmoking female patients with NSCLC. The data showed that the seven overexpressed genes are involved in “module 1” of the subnetwork, which is enriched in the cell cycle pathway. As we know that, ubiquitin-conjugating enzyme E2C (UBE2C) is a key regulator of cell cycle progression. High expression of UBE2C is associated with aggressive progression and poor outcome of malignant glioma,[77]^21 and UBE2C has been identified as a prognostic protein marker in bladder cancer.[78]^22 Moreover, UBE2C, DLGAP5, and TPX2 are associated with the progression and prognosis of pancreatic carcinoma.[79]^23 Additionally, both DLGAP5 and TPX2 are mitosis-associated genes and correlated with poor prognosis in NSCLC patients.[80]^24 TPX2 overexpression is associated with poor survival in gastric cancer.[81]^25 A study also found that TPX2 promotes glioma cell proliferation and invasion via activation of the AKT signaling pathway in glioblastoma multiforme.[82]^26 Elevated DLGAP5 expression was negatively correlated with both OS and relapse-free survival of lung cancer.[83]^27 A study found that silencing of DLGAP5 by siRNA significantly inhibits the proliferation and invasion of hepatocellular carcinoma cells.[84]^28 Cyclin B2 (CCNB2) is a member of cyclin family proteins. Using the same data (GEO:[85]GSE19804), Qian et al further demonstrated that both mRNA and protein expressions of CCNB2 were higher in NSCLC than in normal lung tissues and CCNB2 overexpression is a poor prognostic biomarker in Chinese NSCLC patients.[86]^29 BIRC5, which codes for survivin, can participate in cellular survival functions, such as cell cycle progression.[87]^30 Moreover, BIRC5 (survivin) is a pejorative prognostic marker in stage II/III breast cancer.[88]^31 Similarly, survivin is upregulated in lung cancer tissues and high expression of BIRC5 is associated with poor survival in lung cancer.[89]^32^,[90]^33 Kinesin family member 20A (KIF20A), also known as RAB6KIFL, plays an important role in cytokinesis.[91]^34 In addition, it has been reported that high expression of KIF20A is associated with poor OS in cervical squamous cell carcinoma.[92]^35 Similarly, a study has found that positive expression of KIF20A indicates poor prognosis of glioma patients.[93]^36 High expression of TOP2A, which is involved in DNA synthesis and transcription, is associated with recurrence and metastasis of prostate cancer.[94]^37 In addition, high mRNA levels of TOP2A are independent predictors of poor outcome in renal cell carcinoma patients.[95]^38 Moreover, another study has found that TOP2A is associated with worse prognosis in NSCLC patients.[96]^39 Taken together, these data suggest that UBE2C, DLGAP5, TPX2, CCNB2, BIRC5, KIF20A, and TOP2A are involved in the cell cycle pathway and play a significant role in cancer development, which supports our findings. On the other hand, GNG11 and ANXA1, which are found in “module 2”, were associated with chemokine signaling pathway, cytokine-cytokine receptor interaction, and neuroactive ligand-receptor interaction. GNG11, a lipid-anchored protein, acts as a tumor suppressor in lung adenocarcinoma.[97]^40 Annexin A1 (ANXA1) is a calcium- and phospholipid-binding protein. Some studies found that low ANXA1 expression was associated with better OS in gastric cancer[98]^41 and esophageal squamous cell carcinoma.[99]^42 However, low expression of ANXA1 was correlated with worse OS in breast cancer patients.[100]^43 Downregulation of ANXA1 is correlated with radioresistance in nasopharyngeal carcinoma.[101]^44 Moreover, decreased expression of ANXA1 was correlated with poor survival of pancreatic ductal adenocarcinoma (PDAC) patients, and ANXA1 knockdown inhibited cell proliferation, induced G1 phase cell cycle arrest, and increased PDAC cell migration and invasion capacity compared with controls.[102]^45 Similarly, ANXA1 knockdown suppressed the proliferation, migration, and invasion of NSCLC cells.[103]^46 Together, we speculate that GNG11 and ANXA1 may play a crucial role in NSCLC. ANXA1 could regulate the gastric cancer cell invasion through the formyl peptide receptor (FPR)/extracellular signal-regulated kinase/integrin beta-1-binding protein pathway, and all three FPRs (FPR1 through FPR3) were involved in the regulation process.[104]^41 FPR2 promotes invasion and metastasis of gastric cancer cells and predicts the prognosis of patients.[105]^47 A study also found that FPR2 was overexpressed in epithelial ovarian cancer (EOC) tissues and FPR2 overexpression indicated poor prognosis of EOC patients.[106]^48 However, in the present study, we found that FPR2 was downregulated and low expression indicated better prognosis of NSCLC patients. Similarly, we found that Neuromedin U (Nmu), a secreted neuropeptide, was upregulated and high expression indicated better prognosis of NSCLC patients. But in HNSCC, the expression levels of Nmu in primary tumors with regional metastasis were higher, compared with those without metastasis, and overexpression of Nmu may be involved in the process of regional metastasis of HNSCC.[107]^49 So, the roles of FPR2 and Nmu need to be further investigated. In summary, the present study was intended to identify the potential biomarkers and analyze the prognostic values in nonsmoking females with NSCLC by bioinformatics analysis. We found that hub genes of complex networks, such as UBE2C, DLGAP5, TPX2, CCNB2, BIRC5, KIF20A, TOP2A, GNG11, and ANXA1, may act as potential biomarkers for nonsmoking females with NSCLC. However, the current study was performed by bioinformatics analysis and the conclusions remain to be confirmed by corresponding experiments. Therefore, further investigation is required to verify our findings and determine the potential clinical value of these as biomarkers. Footnotes Disclosure The authors report no conflicts of interest in this work. References