The gene expression profiles of GSE19804 and GSE10072 were obtained from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo). These two gene expression datasets were analyzed using the Affymetrix platform (Affymetrix Human Genome U133 Plus 2.0 Array; Thermo Fisher Scientific, Inc., Waltham, MA, USA). The GSE19804 gene expression profile submitted by Lu et al (8) included 60 pairs of clinical NSCLC samples, which consisted of 56 adenocarcinoma, 3 bronchioloalveolar carcinoma, and 1 squamous carcinoma, and corresponding adjacent normal tissue samples. The GSE10072 gene expression profile consisted of 58 adenocarcinoma samples (16 non-smokers, 18 former smokers and 24 current smokers) and 49 non-tumor samples (15 non-smokers, 18 former smokers and 16 current smokers) (9).

Raw microarray data files of the two datasets were downloaded from the GEO database. GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), an online tool that compares two or more groups of samples in the same experimental setting, was used to analyze the raw data (10). False Discovery Rate (FDR) adjusted P-value of 0.05 and |logFC|>1 were set as the cut-off criteria.

Gene ontology (GO) analysis was processed by the Database for Annotation, Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/) to elucidate the biological function of genes in NSCLC. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed to identify DEGs using the DAVID database. P<0.05 was set as the threshold.

The functional interaction of proteins can shed light on the molecular mechanism underlying NSCLC. The online database STRING (version 10.0, http://string.embl.de/) can be used in the evaluation of PPIs (11). The STRING database includes 9,643,763 proteins from 2,031 organisms. In order to evaluate the PPIs among the DEGs, DEGs were mapped to the STRING database. A confidence score >0.7 was selected as significant. In addition, the degree of the nodes in PPI network was calculated, and the nodes with a higher degree were selected as hub proteins. Furthermore, Cytoscape software (version 3.4.0, http://cytoscape.org/) was employed to construct PPI networks. The plug-in Molecular Complex Detection (MCODE) was performed to screen modules of the PPI network with the threshold set as follows: MCODE scores >10. The GO and KEGG analysis of genes in the module was performed using the DAVID online tool as aforementioned.

The cell lines, human bronchial epithelial (HBE1), A549 and H322, were gifted from Professor Zeyao Tang (Dalian Medical University, Dalian, China) (12). The cells were maintained in high-glucose Dulbecco's modified Eagles medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (GE Healthcare Life Sciences, Logan, UT, USA), 100 U/ml penicillin and 100 µg/ml streptomycin (Gibco; Thermo Fisher Scientific, Inc.). The cells were incubated at 37°C in a humidified chamber with 5% CO₂.

Total RNA from cells lines, including HBE1, A549, and H322, were extracted by using the TRIzol^® regent (Invitrogen; Thermo Fisher Scientific, Inc.). The cDNA of mRNA was synthesized using the PrimeScript^™ RT reagent kit (Takara Bio Inc., Otsu, Japan). RT-qPCR was carried out using the 7500 Real-time PCR system (Thermo Fisher Scientific, Inc.) at 95°C for initial denaturation for 10 min, followed by 40 cycles at 95°C for 15 sec, and 60°C for 1 min with the SYBR^® Green mix (Takara Bio Inc., Japan). Data were analyzed by using the comparative Cq (ΔΔCq) to determine the relative gene expression, and GAPDH was used as an endogenous control (13). The primers were synthesized by Shanghai GenePharma Co., Ltd., (Shanghai, China). The following primer pairs was used to measure the amount of GAPDH: Forward, 5′-GGAGCGAGATCCCTCCAAAAT-3′ and reverse, 5′-GGCTGTTGTCATACTTCTCATGG-3′.

ROC curve analysis was performed using the MedCalc software packages (version 16.8.4; MedCalc Software bvba, Ostend, Belgium). The area under the curve (AUC) values with 95% confidence interval (CI) were calculated to evaluate the overall performance of the diagnostic tests.

Kaplan-Meier plotter (www.kmplot.com), an online survival analysis tool, was used to evaluate the prognostic value of biomarkers of breast, ovarian, lung and gastric cancer (14). Patients with NSCLC were divided into high and low expression groups using the median level, which was included in the low expression group, as the cutoff value. To analyze the association between gene expression and clinical outcomes, Kaplan-Meier plots was employed to compare the overall survival ratio between the two groups, and the log rank P-value and hazard ratio (HR) with 95% confidence intervals (CI) were calculated and displayed.

The data are expressed as the mean ± standard deviation of three replicates. Statistical differences were assessed using one-way analysis of variance test and Tukey's multiple comparisons test. SPSS software (version 17.0; SPSS, Inc., Chicago, IL, USA) was used to analyze the data. P<0.05 was considered to indicate a statistically significant difference.

By using the threshold (adjusted P-vale of 0.05 and fold change >2), a total of 1,412 DEGs were identified in the GSE19804 dataset. Among these genes, 453 genes were upregulated, and 959 genes were downregulated. A heat-map illustrating the expression of the top 50 up and downregulated DEGs is shown in Fig. 1.

To further elucidate the functions of the identified DEGs in NSCLC, GO and KEGG pathway enrichment analyses were employed. As shown in Table I, GO analysis of upregulated DEGs in NSCLC indicated that these genes were associated with ‘mitotic cell cycle’, ‘mitotic nuclear division’ and the ‘cell cycle process’. KEGG pathway enrichment analysis of DEGs revealed that upregulated DEGs were largely enriched in cell cycle and extracellular matrix (ECM)-receptor interaction pathways, while downregulated DEGs were enriched in ‘malaria’ and ‘tumor necrosis factor (TNF) signaling pathways’ (Table I). These results suggest that upregulated DEGs in NSCLC may be largely involved in the progression of the cell cycle.

Based on the analysis of DEGs in the STRING database, a PPI network of DEGs containing 1,291 nodes and 2,854 edges was constructed. By using the plug-in MCODE in Cytoscape, the top 3 modules in the PPI network was obtained (Fig. 2A-C), and KEGG analysis of genes in the corresponding modules was also performed (Fig. 2D-F). Consistent with the KEGG analysis of DEGs, function enrichment analysis of genes in the top 3 modules indicated that these hub genes were also enriched in ‘cell cycle progression’ (Fig. 2E). Therefore, the present study focused on the 5 hub genes associated with cell cycle progression including CDC20, CENPF, KIF2C, BUB1 and ZWINT.

Although 5 hub genes were selected by KEGG analysis the genes in these 3 modules, these 5 selected genes may be limited to the diagnosis or prognosis for non-smoking female patients with NSCLC. In order to elucidate whether these genes can be non-selectively applied to patients with NSCLC, as previously reported (4,15), an additional dataset and RT-qPCR were employed to validate the mRNA level of these genes in NSCLC samples and cell lines. Since the GSE19804 dataset included 56 non-smoking female adenocarcinoma samples (8), the present study searched for a dataset that included adenocarcinoma and simultaneously excluded the effects of sex and smoking. Based on the aforementioned criterion, the GSE10072 database was identified as suitable. Using the GSE10072 dataset, it was detected that the mRNA level of these 5 genes were also overexpressed in NSCLC samples (Fig. 3A-E). In addition, the RT-qPCR results also validated that the mRNA level of these genes were overexpressed in NSCLC cell lines including A549 and H322 (16), when compared with the control cell line HBE1 (Fig. 3F). H322 may be identical to another uncommonly used NSCLC cell line H322M (https://web.expasy.org/cellosaurus/CVCL_1556). Taken together, these results suggest that these 5 hub genes may be novel gene signatures for patients with NSCLC.

To evaluate the diagnostic value of these 5 hub genes, ROC analysis was conducted based on these 2 datasets. The present study demonstrated that the sensitivity and specificity of these 5 genes was relatively high. As shown in Fig. 4A, the AUC values for CDC20, CENPF, KIF2C, BUB1 and ZWINT were 0.927, 0.906, 0.887, 0.876, and 0.937, respectively in the GSE19804 dataset, while the values were 0.958, 0.944, 0.923, 0.897, and 0.942, respectively in the GSE10072 dataset (Fig. 4B). These results indicate that these 5 hub genes may be sensitive and specific in distinguishing NSCLC tissues from normal tissues.

The prognostic value of these 5 genes in PPI network was evaluated using the Kaplan-Meier plotter as previous described (14). Based on the low and high expression of each hub gene, the overall survival of patients with NSCLC was obtained for each gene. As shown in Fig. 5, the high mRNA expression of CDC20 (HR, 1.82; CI, 1.6–2.07) was associated with a poorer overall survival for patients with NSCLC. Similar associations were detected for: CENPF (HR, 1.57, CI, 1.38–1.78), KIF2C (HR, 1.78; CI, 1.57–2.03), BUB1 (HR, 1.83; CI, 1.61–2.09) and ZWINT (HR, 1.5; CI, 1.32–1.71). These results indicate that these 5 hub genes may serve as potential prognostic biomarkers for patients with NSCLC.