Identification of genes and analysis of prognostic values in nonsmoking females with non-small cell lung carcinoma by bioinformatics analyses

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 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.¹ 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,² only 7% of female patients with lung cancer have a history of cigarette smoking in Taiwan.³ In addition, other factors, such as environmental exposure^4,5 or hereditary factors,⁶ 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,⁷ EGFR,⁸ TP53,⁹ and PIK3CA,¹⁰ 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.^11,12 Similarly, some studies were performed to find important key genes in lung cancer.^13–15 But there are limited studies on nonsmoking females with NSCLC.

In 2010, Lu et al³ 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 GSE19804 were obtained from the Gene Expression Omnibus (GEO) database. GSE19804, which was based on the platform of the 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.³

Identification of DEGs

The DEGs between NSCLC samples and normal controls were identified using GEO2R (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),¹⁶ 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¹⁷ 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¹⁸ in Cytoscape.¹⁹ 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 (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).²⁰ 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 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

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 (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 (Table 2). Subsequently, two modules with MCODE scores >14 and the number of nodes >15 were selected (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 (Table 3).

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

Abbreviation: 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

Abbreviations: BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Figure 1 PPI network of DEGs.

Abbreviatins: DEGs, differentially expressed genes; PPI, protein–protein interaction.

Figure 2 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.

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) (Figure 3). Additionally, low expression of ANXA1 was associated with poorer OS for NSCLC patients, as well as GNG11 (P<0.05) (Figure 3).

Figure 3 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,²¹ and UBE2C has been identified as a prognostic protein marker in bladder cancer.²² Moreover, UBE2C, DLGAP5, and TPX2 are associated with the progression and prognosis of pancreatic carcinoma.²³ Additionally, both DLGAP5 and TPX2 are mitosis-associated genes and correlated with poor prognosis in NSCLC patients.²⁴ TPX2 overexpression is associated with poor survival in gastric cancer.²⁵ A study also found that TPX2 promotes glioma cell proliferation and invasion via activation of the AKT signaling pathway in glioblastoma multiforme.²⁶ Elevated DLGAP5 expression was negatively correlated with both OS and relapse-free survival of lung cancer.²⁷ A study found that silencing of DLGAP5 by siRNA significantly inhibits the proliferation and invasion of hepatocellular carcinoma cells.²⁸ Cyclin B2 (CCNB2) is a member of cyclin family proteins. Using the same data (GEO: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.²⁹ BIRC5, which codes for survivin, can participate in cellular survival functions, such as cell cycle progression.³⁰ Moreover, BIRC5 (survivin) is a pejorative prognostic marker in stage II/III breast cancer.³¹ Similarly, survivin is upregulated in lung cancer tissues and high expression of BIRC5 is associated with poor survival in lung cancer.^32,33 Kinesin family member 20A (KIF20A), also known as RAB6KIFL, plays an important role in cytokinesis.³⁴ In addition, it has been reported that high expression of KIF20A is associated with poor OS in cervical squamous cell carcinoma.³⁵ Similarly, a study has found that positive expression of KIF20A indicates poor prognosis of glioma patients.³⁶ High expression of TOP2A, which is involved in DNA synthesis and transcription, is associated with recurrence and metastasis of prostate cancer.³⁷ In addition, high mRNA levels of TOP2A are independent predictors of poor outcome in renal cell carcinoma patients.³⁸ Moreover, another study has found that TOP2A is associated with worse prognosis in NSCLC patients.³⁹ 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.⁴⁰ 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⁴¹ and esophageal squamous cell carcinoma.⁴² However, low expression of ANXA1 was correlated with worse OS in breast cancer patients.⁴³ Downregulation of ANXA1 is correlated with radioresistance in nasopharyngeal carcinoma.⁴⁴ 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.⁴⁵ Similarly, ANXA1 knockdown suppressed the proliferation, migration, and invasion of NSCLC cells.⁴⁶ 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.⁴¹ FPR2 promotes invasion and metastasis of gastric cancer cells and predicts the prognosis of patients.⁴⁷ A study also found that FPR2 was overexpressed in epithelial ovarian cancer (EOC) tissues and FPR2 overexpression indicated poor prognosis of EOC patients.⁴⁸ 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.⁴⁹ 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.

Disclosure

The authors report no conflicts of interest in this work.