Identification of key differentially expressed mRNAs and microRNAs in non-small cell lung cancer using bioinformatics analysis

Non-small cell lung cancer (NSCLC) is a leading cause of mortality worldwide. However, the pathogenesis of NSCLC remains to be fully elucidated. Therefore, the present study aimed to explore the differential expression of mRNAs and microRNAs (miRNAs/miRs) in NSCLC and to determine how these RNA molecules interact with one another to affect disease progression. Differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) were identified from the GSE18842, GSE32863 and GSE29250 datasets downloaded from the Gene Expression Omnibus (GEO database). Functional and pathway enrichment analysis were performed based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. STRING, Cytoscape and MCODE were applied to construct a protein-protein interaction (PPI) network and to screen hub genes. The interactions between miRNAs and mRNAs were predicted using miRWalk 3.0 and a miRNA-mRNA regulatory network was constructed. The prognostic value of the identified hub genes was then evaluated via Kaplan-Meier survival analyses using datasets from The Cancer Genome Atlas. A total of 782 DEGs and 46 DEMs were identified from the 3 GEO datasets. The enriched pathways and functions of the DEGs and target genes of the DEMs included osteoclast differentiation, cell adhesion, response to a drug, plasma membrane, extracellular exosome and protein binding. A subnetwork composed of 11 genes was extracted from the PPI network and the genes in this subnetwork were mainly involved in the cell cycle, cell division and DNA replication. A miRNA-gene regulatory network was constructed with 247 miRNA-gene pairs based on 6 DEMs and 210 DEGs. Kaplan-Meier survival analysis indicated that the expression of ubiquitin E2 ligase C, cell division cycle protein 20, DNA topoisomerase IIα, aurora kinase A and B, cyclin B2, maternal embryonic leucine zipper kinase, slit guidance ligand 3, phosphoglucomutase 5, endomucin, cysteine dioxygenase type 1, dihydropyrimidinase-like 2, miR-130b, miR-1181 and miR-127 was significantly associated with overall survival of patients with lung adenocarcinoma. In the present study, a miRNA-mRNA regulatory network in NSCLC was established, which may provide future avenues for scientific exploration and therapeutic targeting of NSCLC.


Introduction
Lung cancer is one of the most common types of cancer and is the leading cause of cancer-associated mortality worldwide, accounting for >25% of all deaths due to cancer (1). Approximately 85% of lung cancers are non-small cell lung cancer (NSCLC) (2). In turn, NSCLCs are comprised of lung adenocarcinoma (LUAD, 40-50% of lung cancers), lung squamous cell carcinoma (LUSC, 20-30% of lung cancers) and large cell carcinoma (9% of lung cancers) (2). Despite advances in diagnostics, surgery and medication in recent decades, the average 5-year survival rate of patients with NSCLC remains as low as 15% (1). This poor prognosis is a consequence of high rates of tumor metastasis and recurrence, and numerous signaling pathways having been identified to be involved in these processes (1). Thus, an enhanced understanding of the molecular mechanisms controlling NSCLC progression is required to improve the low survival rate. The development of high-throughput sequencing has allowed for comprehensive comparisons of gene expression profiles, thereby identifying differentially expressed genes (DEGs) between tumor and normal tissues. Changes in expression levels usually indicate pathological states, as proteins encoded by DEGs may be involved in tumorigenesis and tumor progression (3). microRNAs (miRNAs/miRs) are short non-coding RNA molecules that mediate the post-transcriptional regulation of mRNAs via binding to complementary sequences in the 3 ' -untranslated region (UTR) of mRNAs and suppressing their translation or mediating their degradation (4). An individual miRNA may regulate hundreds of different mRNA molecules, highlighting the existence of miRNA-mRNA regulatory networks. Depending on how it is expressed, a specific miRNA may therefore act to suppress or promote oncogenesis via specific effects on relevant target mRNAs (5). Indeed, miRNA profiling efforts have been used to identify specific miRNA signatures associated with particular tumor subtypes, thereby allowing for cancer diagnosis, treatment planning and prediction of patient prognosis (6). Bioinformatics analysis of gene expression microarray data provides a useful tool for revealing numerous previously unrecognized mRNAs and miRNAs that may be implicated in the pathogenesis of cancer or other diseases (7).
In the present study, 3 gene expression datasets were analyzed using an integrated bioinformatics approach in order to identify DEGs and differentially expressed miRNAs (DEMs) between NSCLC tumors and healthy control tissues. Functional enrichment and protein-protein interaction (PPI) network analyses were performed to better establish the functions of these mRNAs and miRNA, and these approaches were combined with an analysis of mRNA-miRNA interactions to screen hub genes and miRNAs in this regulatory network. Through this approach, the present study aimed to further elucidate the molecular mechanisms of NSCLC to identify potentially novel therapeutic strategies for its treatment.

Materials and methods
Microarray data information. The datasets used in the present study were obtained from the Gene Expression Omnibus (GEO) database (https://www. ncbi.nlm.nih.gov/geo/) (8). The original gene expression profiles from the datasets GSE18842 (9), GSE32863 (10) and GSE29250 (11) were used and the clinical information of the patients was obtained from the original research articles. The GSE18842 dataset included 91 samples (46 tumors and 45 controls) and all samples were paired except 2 tumors and 1 control. The platform used for the GSE18842 dataset was GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. The GSE32863 dataset included 116 samples (58 tumors and 58 controls) and all samples were paired. The platform used for the GSE32863 dataset was the GPL6884 Illumina Human WG-6 v.3.0 expression beadchip. The GSE29250 dataset included 12 samples (6 tumors and 6 controls) and all samples were paired. The platform used for the GSE29250 dataset was the GPL8179 Illumina Human v.2 MicroRNA expression beadchip.

Identification of DEGs.
To compare gene expression profiles, the GEO2R tool (http://www.ncbi.nlm.nih.gov/geo/geo2r/; accessed March 2019) (12), which is based on the limma package in R, was used to individually identify DEMs and DEGs in each dataset. To control for type I error as a result of multiple comparisons within each dataset, the false discovery rate (FDR) determination feature automatically included in the GEO2R tool was employed. Significant DEGs were those that remained significant after FDR correction when tested via multiple-comparisons t-tests, and fold change (FC)>2 and P<0.05 were set as the cut-off criteria. Any probes that did not correspond to a specific gene symbol were then filtered from the resultant data.
Gene function analysis. For gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of DEMs and DEGs, the database for annotation, visualization and integrated discovery (DAVID) database (v6.8; http://david.abcc.ncifcrf.gov/) was used, focusing specifically on humans and using all genes as an enrichment background. Significant enrichment in these analyses was determined based on an adjusted P-value of 0.05 as established via the Benjamini-Hochberg method (13). These P-values were determined on the basis of a cumulative hypergeometric distribution, calculating q-values based on the Benjamini-Hochberg procedure as a means of controlling for multiple testing (13). For comparisons of hierarchical clustering of enriched terms, clusters were designated as groupings that had a similarity score >0.3, with the most significant term within a given cluster being selected to represent the cluster as a whole.
PPI network construction and analysis. The online Search Tool for the Retrieval of Interacting Genes and proteins (STRING) database (v.11.0; https://string-db.org/) was used for PPI network construction. PPI pairs with a combined score ≥0.7 were used to generate the network. Cytoscape (v.3.4.0; https://cytoscape.org/) was used to visualize the regulatory interactions between these genes and CentiScaPe (v.2.2; Center for Biomedical Computing, University of Verona, Italy) (14) was used to analyze network distributions based on topological properties. The Molecular Complex Detection application (MCODE; v.1.6) was used to identify and extract subnetworks from the global PPI network based on the k-core algorithm (15). The genes with a degree ≥30 in this regulatory network were identified as hub genes, as described previously (16).
Prediction of the mRNA-miRNA interactions. An online tool called miRWalk 3.0 (http://mirwalk.umm.uni-heidelberg.de/), which integrates predictive outputs of TargetScan (17) and miRDB (18), was used to predict DEG and DEM interactions. A score ≥0.95 was considered as the critical criterion for the miRWalk predictive analysis. Only the target mRNAs predicted by all 3 tools (miRWalk, TargetScan and miRDB) were used for further analysis. By overlaying identified DEGs and these predicted mRNA targets, a miRNA-mRNA regulatory network was constructed and then visualized using Cytoscape.

Analysis of datasets from the cancer genome atlas (TCGA).
TCGA is an online database that may be used to research and explore publicly available datasets (https://cancergenome.nih. gov/; accessed March 2019) (19), including RNA sequencing (RNA-seq) data from TCGA samples of 31 different types of cancer. In the present study, the tumor types were limited to LUAD and LUSC. RNA-seq data and clinical data from 478 LUAD and 482 LUSC samples from TCGA datasets were used. Based on the approach previously outlined by Li and Dewey (20), a PERL program was used to multiply the 'scaled estimate' by 10 6 , yielding transcripts per million (TPM) values for all gene expression, as TPM values were thought to be a more reliable means of comparing gene expression than fragments per kilobase of TPM-mapped reads or reads per kilobase of TPM-mapped reads values (21). In the present study, to improve the reliability of the analysis, the expression of hub genes was validated in TCGA datasets using Gene Expression Profiling Interactive Analysis (GEPIA; v1.0; http:// gepia.cancer-pku.cn/). For each of the hub genes, patients were stratified into 2 groups based on expression levels of each gene and differences in patient survival were analyzed, generating hazard ratio (HR) and 95% CI values, as well as log-rank P-values for each comparison.

Identification of DEGs and DEMs.
In the GSE18842 and GSE32863 datasets, 3,167 and 1,270 DEGs were identified, respectively, of which 1,395 and 514 were upregulated, and 1,772 and 756 were downregulated ( Fig. 1A and B). A total of 782 DEGs were shared between these 2 datasets (232 upregulated and 550 downregulated). The GSE32863 dataset yielded a list of 46 DEMs, of which 26 were upregulated and 20 were downregulated (Fig. 1C).
Functional enrichment analysis of overlapped DEGs and target genes of DEMs. To assess the biological roles of these DEGs and target genes of DEMs, KEGG and GO enrichment analyses were performed. The top 10 enriched terms for each analysis were compiled in Fig. 2. KEGG pathway enrichment analysis indicated that the DEGs and target genes of DEMs were mainly enriched in osteoclast differentiation, complement and coagulation cascades, Staphylococcus aureus infection and pertussis ( Fig. 2A). GO analysis in the category biological process suggested these DEGs and target genes of DEMs were primarily enriched in 'cell adhesion,' 'response to drugs' and 'extracellular matrix organization' (Fig. 2B). GO analysis in the category cellular component suggested that the DEGs and target genes of DEMs were mainly enriched in 'plasma membrane,' 'extracellular exosome' and 'extracellular localization' (Fig. 2C). In the category molecular function, the DEGs and target genes of DEMs were mostly enriched in 'protein binding,' 'identical protein binding' and 'calcium ion binding' (Fig. 2D).
PPI network construction and analysis of modules. The 782 overlapping DEGs which were shared between the 2 datasets (GSE18842 and GSE32863) indicated a distinct set of interactions and networks. PPI pairs with a combined score ≥0.7 were used to generate the network. A PPI network was constructed using 445 of the 782 DEGs and the resultant network had 445 nodes and 1,490 edges (Fig. 3A). There were 137 upregulated and 308 downregulated genes among the 445 DEGs. A total of 11 nodes had a degree of >30 and were designated as hub genes, including interleukin 6, Jun proto-oncogene (JUN),  Table I. Furthermore, a subnetwork containing 25 nodes and 284 edges was extracted from the global PPI network (Fig. 3B). The results of the KEGG and GO analyses for the genes in the subnetwork are presented in Table II. The most significantly enriched terms in this network were 'cell cycle, division' and 'DNA replication' associated with cancer, confirming the relevance of the present analysis to NSCLC progression and prognosis.
Survival analysis of hub genes. Survival analysis was performed for the 11 hub genes based on TCGA data. Increased expression of 7 of the hub genes (UBE2C, CDC20, TOP2A, AURKA, AURKB, CCNB2 and MELK) was significantly associated with poorer overall survival (OS) in patients with LUAD (Fig. 4).
miRNA-gene network. Using the miRWalk application, putative DEM targets as established by the TargetScan and miRDB databases were identified. A total of 210 putative target mRNAs overlapped with the DEG dataset, yielding 247 miRNA-gene pairs based on 6 DEMs and the 210 DEGs that were putative targets. These were used to construct an overlapping regulatory network (Fig. 5). A total of 4 upregulated DEMs were predicted to downregulate 185 DEGs, whereas decreased expression of 2 DEMs was predicted to be associated with increased expression of 25 DEGs. TCGA-based survival analysis suggested that none of the hub genes were associated with OS in patients with NSCLC. However, specifically among patients with LUAD, those with elevated expression of slit guidance ligand 3 (SLIT3), phosphoglucomutase 5 (PGM5), endomucin (EMCN), cysteine dioxygenase type 1 (CDO1), dihydropyrimidinase-like 2 (DPYSL2), miR-130b and miR-1181 had a significantly higher OS compared with patients who had low expression of these genes. By contrast, elevated miR-130b and miR-127 expression was associated with poorer OS in patients with LUAD (Fig. 6). Furthermore, the differences of the DEGs and DEMs in miRNA-gene networks of NSCLCs with or without KRAS mutation were examined. The results suggested that the expression levels of TOP2A (P=0.0357; Fig. 7A), AURKA (P=0.0409; Fig. 7B) and MELK (P=0.0190; Fig. 7C) were significantly lower in KRAS mutation groups compared to KRAS wild-type groups (Fig. 7).

Discussion
Cancer is a genetic disease wherein cumulative mutations drive the multi-step progression towards oncogenesis, eventually culminating in unrestrained cancer growth. NSCLC remains one of the most common and deadliest forms of cancer, making the elucidation of the molecular mechanisms governing this disease paramount (22). In the present study, DEMs and DEGs associated with NSCLC were identified via a bioinformatics analysis, yielding 782 DEGs and 46 DEMs    based on overlapping hits in the GSE18842, GSE32863 and GSE29250 datasets. These hits included 232 upregulated and 550 downregulated genes, as well as 26 upregulated and 20 downregulated miRNAs. Through functional enrichment analyses, it was determined that these DEGs were primarily associated with processes including 'osteoclast differentiation', 'complement and coagulation cascades', 'cell adhesion, drug responses', 'plasma membrane', 'extracellular exosome' and 'protein binding'. In addition, a DEG PPI network was generated and a significant subnetwork module was identified that contained genes associated with the cell cycle, DNA replication and oocyte meiosis, with the GO terms enrichment for 'mitotic nuclear division', 'cell division', 'G2/M transition of mitotic cell cycle', 'spindle', 'midbody', 'nucleoplasm', 'ATP binding', 'protein binding' and 'microtubule binding'. Cell cycle dysregulation is known to be a key factor linked to tumor development and progression (23,24). Recent studies indicated that microtubule binding is linked to tumor metastasis and drug resistance (25,26). Complement activation and coagulation cascade activation are similarly able to promote tumor development as a consequence of their ability to mediate the recruitment of myeloid cells that support tumor growth (27). To summarize, the identified DEGs may regulate the proliferation, invasion, migration and drug-resistance of cancer cells through these pathways, thus affecting the occurrence and development of NSCLC. The investigation of these DEGs may pave a way towards novel targeted therapies for NSCLC. Based on the PPI network, 11 hub genes with high degrees of interaction (degree ≥30) were extracted. A prognostic analysis of these 11 hub genes was performed using the online tool GEPIA. The results revealed that patients with LUAD who had upregulation of UBE2C, CDC20, TOP2A, AURKA, AURKB, CCNB2 and MELK had a worse prognosis.
The expression of UBE2C has been previously indicated to be upregulated in lung cancer (28), and the results of the present study suggested that it was associated with poor survival. Similar observations have previously been made in ovarian cancer (29), breast cancer (30) and gastric cancer (31). UBE2C is involved in the progression of the cell cycle and transcription, and upregulation of UBE2C may induce an enhanced growth and colony formation of tumors (32), as well as decreased autophagy in cancer cells (33). Furthermore, UBE2C reduces the sensitivity of cells to common chemotherapy drugs for lung cancer, including cisplatin (34) and docetaxel (35). UBE2C may be used as a therapeutic target for NSCLC.
Another oncogene in several types of tumor and a hub gene identified in the present study was CDC20 (36). CDC20 is an important regulator of the cell cycle and altered expression or functional impairment may induce mitotic arrest to prevent activation of adenomatous polyposis coli and hence, increase premature anaphase manifesting as aneuploidy in daughter cells (37). CDC20 was observed to be upregulated at mRNA and protein levels in NSCLC, and was significantly correlated with tumor size, pleural invasion and histological classification (38). Of note, knockdown of CDC20 caused inhibition of growth, migration ability and formation of colonies in lung cancer cells, as well as cell cycle arrest in G2/M phase and induction of apoptosis (39), making this oncogene a potential target molecule to address NSCLC therapy. Upregulation of CDC20 has been associated with shorter OS in patients with LUAD, but not in patients with LUSC (40), which is consistent with the results of the present study.
TOP2A encodes for a DNA topoisomerase involved in torsional dynamics during replication and transcription (41), which is also associated with cell proliferation (42). TOP2A has been indicated to be upregulated in numerous types of tumor, including breast, nasopharyngeal and renal cell carcinomas, and is associated with poor prognosis; therefore, TOP2A has important roles in cancer (43)(44)(45). The ability of NSCLC cells to proliferate and invade tissues is associated with elevated TOP2A expression. Several anti-cancer agents have been developed to target this gene (46), and the development of drug resistance has been associated with mutation of TOP2A (47). AURKA and AURKB are highly conserved serine/threonine kinases, the former of which is associated with regulating centrosome duplication and spindle formation (48), and the latter of which is important for regulating chromatin modifications and suppressing cytokinesis (49). NSCLC prognosis is known to be associated with elevated expression of these 2 genes (50,51). In the present study, TCGA dataset analysis suggested that the prognosis of patients with LUAD was associated with AURKA and AURKB, which require further clinical trial validation.
CCNB2 is a cyclin gene that activates cyclin-dependent kinase 1 to drive the G2/M cell cycle transition, and inhibition of CCNB2 leads to cell cycle arrest (52). It has been previously confirmed that CCNB2 is upregulated in tissue and serum samples from patients with NSCLC (53,54). Elevated CCNB2 mRNA levels are known to be closely associated with tumor differentiation grade and histological type, and upregulation of CCNB2 at the protein level has been significantly associated with the degree of differentiation, tumor size, lymph node metastasis, distant metastasis and clinical stage (54,55). Previous studies of NSCLC have suggested that there was no statistically significant correlation between the levels of CCNB2 protein and mRNA in NSCLC (55). The results of the present study indicated that upregulation of CCNB2 mRNA was a poor prognostic biomarker in patients with LUAD, while a previous study suggested that the protein levels of CCNB2 may serve as an independent prognostic marker in NSCLC (54). Therefore, the role of CCNB2 in NSCLC should be further elucidated.
MELK is a serine/threonine kinase that has been indicated to be highly expressed in several human cancer types (prostate, breast, brain, colorectal and gastric cancer) and glioblastoma multiforme stem cells (56). Elevated expression of MELK is associated with the degree of tumor malignancy and with poor survival in cervical cancer (57), breast cancer (58) and gastric cancer (59). Furthermore, the present study suggested that upregulation of MELK is associated with the progression of NSCLC. miRNAs regulate a wide array of target mRNAs via 3 ' -UTR binding and subsequent translational repression. As a  result, complex miRNA-mRNA networks may govern a wide range of biological pathways, making miRNAs critical for the progression of numerous types of cancer (60). In the present study, a total of 46 DEMs were identified, and miRWalk-mediated predictive analyses were performed to identify those DEMs that were predicted to interact with DEGs, yielding 6 hub miRNAs and associated mRNAs, including miR-127-5p, miR134-5p, miR-130b-3p, miR-1181, miR-145-3p, miR-153-3p, CDO1, SLIT3 and PGM5. Survival analyses revealed that dysregulation of miR-127-5p, miR-130b-3p, miR-1181, CDO1, SLIT3, PGM5, EMCN and DPYSL2 were significantly associated with the prognosis of patients with LUAD.
miR-127 has previously been indicated to function as either a promoter or suppressor of cancer development depending on the specific context (61,62). Based on the NSCLC network established in the present study, miR-127 was among the most prominent regulatory miRNAs, suggesting it serves complex regulatory functions in the context of NSCLC. In a previous study, miR-127 expression was indicated to be elevated in LUAD and associated with poor prognosis (63), consistent with the results of the present study. High levels of miR-127 induce epithelial-to-mesenchymal transition, rendering tumor cells with stem cell-like properties, and propagate tumor resistance to epidermal growth factor receptor inhibitor (63). The aggressiveness of the cancer was associated with a circuit involving miR-127, NF-κB and tumor necrosis factor α-induced protein 3, which are markers of inflammation (63).
miR-130b has also been documented in several other types of tumor, with upregulation observed in prostate cancer (64), while downregulation was identified in thyroid carcinomas (65). In the present study, miR-130b-3p expression was determined to be significantly increased in NSCLC and associated with poor survival, although this was specifically restricted to patients with LUAD in the TCGA dataset, warranting further investigation.
miR-134 has been indicated to be differentially regulated in lung cancer and other types of cancer (including gastric cancer, breast cancer and oral cancer), with certain studies reporting increased expression in lung cancer (66,67), while other studies observed that it was downregulated (68,69). miR-134 may function to either promote or suppress tumor progression (69,70), highlighting complex mechanisms warranting further investigation.
miR-1181 has been observed to be downregulated in nasopharyngeal carcinoma, ovarian cancer and pancreatic cancer (71)(72)(73). In addition, miR-1181 inhibited invasion and proliferation via STAT3 in pancreatic cancer cells and inhibited metastasis by modulating the WNT/β-catenin pathway in nasopharyngeal carcinoma (71,73). Thus, the role of miR-1181 in NSCLC requires to be further investigated in the future.
In summary, the present study identified 782 DEGs and 46 DEMs between NSCLC tumor and normal tissues, and a miRNA-mRNA regulatory network was established. Certain hub genes were screened out from the PPI network and miRNA-mRNA regulatory network, including UBE2C, CDC20, TOP2A, AURKA, AURKB, CCNB2, MELK, SLIT3, PGM5, EMCN, CDO1, DPYSL2, miR-130b, miR-1181 and miR-127. These analyses suggested a comprehensive overview of the mechanistic basis of NSCLC, potentially highlighting future avenues for treatment. However, analysis of TCGA datasets indicated that the expression of certain hub genes was only associated with the prognosis in patients with LUAD, which requires further validation. The present results remain to be verified by further clinical investigation in the future.