Integrative Bioinformatics Analysis Reveals Potential Long Non-Coding RNA Biomarkers and Analysis of Function in Non-Smoking Females with Lung Cancer.

BACKGROUND Lung cancer is the most lethal cancer worldwide. The aim of this study was to identify the tumor-related lncRNAs and explore their functions in female non-smokers with lung cancer. MATERIAL AND METHODS The gene expression microarray datasets GSE19804, GSE31210, and GSE31548 were downloaded from the Gene Expression Omnibus database. The differentially-expressed lncRNAs between non-smoking female lung cancer samples and non-tumor lung tissues were identified using GEO2R. RESULTS In total, 25, 40, and 15 differentially-expressed lncRNAs were obtained from GSE19804, GSE31210, and GSE31548 datasets (|logFC| >1, adj. P<0.05), respectively. Eight lncRNAs were screened out in all 3 datasets. Of these, 5 lncRNAs were up-regulated and 3 lncRNAs were down-regulated in lung cancer tissues compared to non-tumor lung tissues. Then, the target miRNAs of aberrantly expressed lncRNAs and target mRNAs corresponding to miRNAs were predicted. Subsequently, the ceRNA network with 8 key lncRNAs, 20 miRNAs, and 38 mRNAs were constructed. Functional and pathway enrichment analysis showed these target genes were mainly enriched in biological processes associated with protein binding, nucleus, metal ion binding, regulation of transcription from RNA polymerase II promoter, nucleic acid binding, cell differentiation, microRNAs in cancer, and the hippo signaling pathway. Survival analysis of these lncRNAs revealed that low LINC00968 (P=0.0067) and TBX5-AS1 (P=0.0028) expression were associated with unfavorable prognosis in never-smoking female lung cancer patients. CONCLUSIONS The present study promotes understanding of the molecular mechanism of the pathogenesis of non-smoking female lung cancer and provides potential biomarkers for diagnosis and treatment.

Background

Lung cancer is the most common and lethal cancer worldwide, with an estimated 1.8 million new cases and 1.6 million deaths from this disease in 2012 [1]. In China, there were approximately 610 000 deaths due to lung cancer in 2015 [2]. Lung cancer has become a serious public health challenge and although research has focussed on developing new drug and treatment modalities for decades, the prognosis of lung cancer patients with recurrence or metastasis remains very poor.

Smoking is considered to the major risk factor for lung cancer [3]. However, only a small proportion of female lung adenocarcinoma patients are associated with smoking [4,5]. Molecular-level studies have indicated that some genes are correlated with lung adenocarcinoma in never-smokers, such as RET [6], EGFR [7], SEMA5A [8], PIK3CA [9], and KRAS [10]. In addition, the molecular mechanisms of tumorigenesis in female lung cancer patients remain unclear. lncRNAs, a class of RNA molecules, 200 to hundreds of thousands of nucleotides long, are deregulated in a variety of diseases and are involved in various biological processes [11-14]. To date, however, little is known about the roles of lncRNAs in non-smoking female lung cancer.

In this study, the microarray data (GSE19804, GSE31210, and GSE31548) were obtained from the Gene Expression Omnibus (GEO, ) database. Then, we determined the differentially expressed genes (DEGs) by using bioinformatics analysis. Furthermore, target miRNAs of differentially expressed lncRNAs and target mRNAs corresponding to miRNAs were predicted. Subsequently, the interaction network among lncRNAs, miRNAs and mRNAs were constructed in non-smoking female lung cancer. This study can promote our understanding of the role of lncRNAs and associated pathways in non-smoking female lung cancer.

Material and Methods

Microarray data

Three gene expression profiles (GSE19804, GSE31210, and GSE31548) were obtained from the Gene Expression Omnibus (GEO, ) database. GSE19804 and GSE31210 datasets were based on Affymetrix Human Genome U133 Plus 2.0 Array, and GSE31548 dataset was performed by Affymetrix Human Genome U133B Array. The array data of GSE19804 included 60 non-smoking female lung cancer tissue samples and 60 non-smoking healthy female lung tissues [8]. GSE31210 contained 246 samples, consisting of 98 never-smoking female lung cancer samples and 4 normal female lung tissues [15]. GSE31548 contained 50 samples, including 6 never-smoking female lung cancer samples and 5 normal female lung tissues.

Identification of differentially expressed genes

The GEO2R () was used to identify differentially expressed genes (DEGs) between lung cancer and normal samples. GEO2R is an interactive web analysis tool based on R that comes with the GEO database and can analyze almost all GEO data. The P value was adjusted (adj. P) with Benjamini and Hochberg method by default to obtain false discovery rate (FDR). The absolute fold change (|logFC|) >1 and adj. P<0.05 were set as the cut-off criterion. The gplots package was used to generate the heat maps under the R environment (version 3.3.4, ).

Construction of lncRNA-miRNA-mRNA ceRNA network

The lncRNA-miRNA-mRNA ceRNA network was constructed to explore the association among lncRNA, miRNA, and mRNA based on the hypothesis of ceRNA [16]. The target miRNAs of the differentially expressed lncRNAs were predicted using miRNA tools miRecords () and starBase v2.0 database () [17,18]. The mRNAs targeted by miRNAs were predicted by miRecords and starBase v2.0 database, and then the predicted results and DEGs of GSE19804, GSE31210, and GSE31548 datasets were combined to gain the intersection mRNAs. Subsequently, the ceRNA network was established and visualized using Cytoscape (v3.40) software. The intersection mRNAs were determined by pathway enrichment analysis using the DAVID database () [19].

Survival analysis of differentially expressed lncRNAs

Kaplan-Meier Plotter () is a public online database that can assess the relationship between many genes and breast or ovarian or lung or gastric cancer patient survival [20]. This database comprises clinical information and gene expression data for 2437 lung cancer patients. The overall survival (OS) of non-smoking female lung cancer was detected by a Kaplan-Meier plot. Hazard ratio (HR) and its 95% confidence intervals (CI) were computed and displayed on the page. P<0.05 was considered statistically significant.

Results

Differentially expressed lncRNAs in non-smoking female lung cancer patients

A total of 25, 40, and 15 differentially expressed lncRNAs were identified from GSE19804, GSE31210, and GSE31548 datasets (|logFC| >1, adj. P<0.05), respectively (Figure 1). Eight lncRNAs were screened out in all 3 datasets (Figure 2). Among them, 5 lncRNAs (AFAP1-AS1, LINC00467, LINC00511, LINC00673 and LINC01207) were up-regulated and 3 lncRNAs (LINC00472, LINC00968, and TBX5-AS1) were down-regulated in lung cancer tissues compared to non-tumor lung tissues (Figure 2, Table 1).

Figure 1

Heat map of differentially expressed lncRNAs in GSE19804 (A), GSE31210 (B), and GSE31548 (C) datasets. Red: up-regulation; green: down-regulation.

Figure 2

Identification of differentially expressed lncRNAs in datasets GSE19804, GSE31210, and GSE31548.

Table 1

Differentially expressed lncRNAs in datasets GSE19804, GSE31210 and GSE31548.

GEO accession	Number	lncRNAs
GSE19804GSE31210GSE31548	8	AFAP1-AS1, LINC00467, LINC00472, LINC00511, LINC00673, LINC00968, LINC01207, TBX5-AS1
GSE19804GSE31210	6	FENDRR, LINC00341, LINC00622, SYNE3, TINCR, WASIR2
GSE19804GSE31548	1	RNF157-AS1
GSE19804	10	BRE-AS1, DRAIC, EP300-AS1, HOXB-AS3, LCAL1, LINC00312, LINC00551, LINC01140, MAGI2-AS3, MBNL1-AS1
GSE31210	26	ADAMTS9-AS1, BDNF-AS, DEPDC1-AS1, FLVCR1-AS1, H1FX-AS1, LAMTOR5-AS1, LINC00115, LINC00638, LINC00689, LINC00702, LINC00842, LINC00857, LINC00944, LINC00961, LINC01106, LINC01123, LINC01314, LINC01569, MAFG-AS1, RALY-AS1, SBF2-AS1, STK24-AS1, STK4-AS1, SYNPR-AS1, THUMPD3-AS1, TMEM51-AS1
GSE31548	6	ATP13A4-AS1, HAGLR, LINC01088, LINC01279, PAX8-AS1, PAXIP1-AS1

Prediction of lncRNA targets and ceRNA network construction

To construct the ceRNA network, miRecords and starBase v2.0 were used to detect the potential target miRNAs [17,18]. The results indicated 10 aberrantly expressed lncRNAs targeted by 20 specific miRNAs (Table 2). In the next step, based on miRecords and starBase v2.0, the predicted targets of miRNAs described in Table 2 were gained. Subsequently, the predicted results and DEGs were combined to obtain the intersection mRNAs. We identified the targeted relationship between 19 specific miRNAs and the 38 intersection mRNAs (Table 3). Of these, some targets are cancer-related genes, such as RECK, EZH2, LATS2, CXCL12, NR4A2, TPPP3. Next, based on the data described in Tables 2 and 3, the ceRNA network among lncRNA, miRNA, and mRNA was constructed (Figure 3). A total of 8 lncRNAs, 19 miRNAs, and 38 mRNAs were involved in the ceRNA network.

Table 2

miRNAs that may interact with specific lncRNAs.

lncRNAs	miRNAs
AFAP1-AS1	hsa-miR-451a
LINC00511	hsa-miR-29c-3p, hsa-miR-29a-3p, hsa-miR-29b, hsa-miR-16-5p, hsa-miR-15b
LINC00467	hsa-miR-18a-5p, hsa-miR-18b-5p
TBX5-AS1	hsa-miR-92b-3p
LINC00472	hsa-miR-196a, hsa-miR-23a-3p, hsa-miR-302d-3p, hsa-miR-372-3p, hsa-miR-23b-3p,hsa-miR-204-5p
LINC00673	hsa-miR-150-5p, hsa-miR-1231
LINC00968	hsa-miR-26a-5p
LINC01207	hsa-miR-18

Table 3

miRNAs targeting specific mRNAs.

miRNAs	mRNAs
hsa-miR-15b	RECK, WEE1, CACUL1, IPO7
hsa-miR-16-5p	CADM1,RECK,TPPP3, SPDEF
hsa-miR-18	CTGF, TIMP3
hsa-miR-18a-5p	CTGF
hsa-miR-18b-5p	TIMP3
hsa-miR-196a	ANXA1,S100A9
hsa-miR-204-5p	SPDEF, AP1S2
hsa-miR-23a-3p	CXCL12, ZNF138
hsa-miR-23b-3p	PLAU
hsa-miR-26a-5p	EZH2, TGFBR2
hsa-miR-29a-3p	COL3A1, TET3
hsa-miR-29b	COL1A1, COL3A1
hsa-miR-29c-3p	COL1A1, TET3
hsa-miR-302d-3p	KLF13, NR4A2
hsa-miR-372-3p	KLF13, LATS2, NR4A2, RECK
hsa-miR-451a	ABCB1, OSR1
hsa-miR-92b-3p	CDKN1C
hsa-miR-150-5p	ZEB1, MYB, gag-pol, MUC4, EGR2, PFN2
hsa-miR-1231	AZIN1, ZAK, FAM127B, G3BP1, CTGF

Figure 3

The lncRNA-miRNA-mRNA ceRNA network. Red diamonds, miRNAs; Green rectangles, lncRNAs; Orange ovals, mRNAs.

Functional and pathway enrichment analysis

To learn more about the function of identified intersection DEGs in Table 3, functional and pathway analysis was carried out using DAVID [19]. These genes were mainly enriched in biological processes associated with protein binding, nucleus, metal ion binding, regulation of transcription from RNA polymerase II promoter, nucleic acid binding, cell differentiation, microRNAs in cancer, and hippo signaling pathway (Table 4).

Table 4

Functional and pathway enrichment analysis analysis of miRNA target genes in ceRNA network.

Category	Term	Description	Count	P value
GOTERM_BP_DIRECT	GO: 0030154	Cell differentiation	6	0.002
GOTERM_BP_DIRECT	GO: 0000122	Negative regulation of transcription from RNA polymerase II promoter	7	0.002
GOTERM_BP_DIRECT	GO: 0045944	Positive regulation of transcription from RNA polymerase II promoter	8	0.003
GOTERM_CC_DIRECT	GO: 0005615	Extracellular space	9	0.002
GOTERM_CC_DIRECT	GO: 0005634	Nucleus	19	0.002
GOTERM_CC_DIRECT	GO: 0005576	Extracellular region	8	0.022
GOTERM_MF_DIRECT	GO: 0005515	Protein binding	30	5.0 E-5
GOTERM_MF_DIRECT	GO: 0046872	Metal ion binding	10	0.018
GOTERM_MF_DIRECT	GO: 0003676	Nucleic acid binding	6	0.045
KEGG_PATHWAY	hsa05206	MicroRNAs in cancer	6	0.001
KEGG_PATHWAY	hsa04390	Hippo signaling pathway	3	0.027

Kaplan-Meier plotter identified potential prognostic lncRNAs for never-smoking female lung cancer

We used Kaplan-Meier plotter to assess the prognostic values of the 8 lncRNAs in never-smoking female lung cancer patients [20]. P value less than 0.05 was considered a significant statistical difference. The results showed that only low LINC00968 (HR=0.24, 95%CI: 0.08–0.74, P=0.0067) and TBX5-AS1 (HR=0.19, 95% CI: 0.05–0.64, P=0.0028) expression were related to unfavorable prognosis in never-smoking female lung cancer patients (Figure 4). However, there was no significant relationship among the expressions of other lncRNAs and prognosis of non-smoking female lung cancer.

Figure 4

Kaplan-Meier survival curves for LINC00968 and TBX5-AS1 expression in never-smoking female lung cancer patients.

Discussion

In this study, we first analyzed the microarray data of non-smoking female lung cancer from GEO database under the accession number GSE19804, GSE31210, and GSE31548 by GEO2R to obtain the DEGs. GEO2R is a powerful tool to process gene expression data and can analyze almost all GEO data. Then, the target miRNAs of these lncRNAs and the mRNAs targeted by miRNAs were predicted. Finally, the lncRNAs-miRNAs-mRNAs ceRNA network was constructed. This study provides important clues for exploring the key lncRNAs and associated regulatory network in pathogenesis of non-smoking female lung carcinoma.

According to our results, a total of 8 differentially expressed lncRNAs including 5 up-regulated lncRNAs and 3 down-regulated lncRNAs were identified in all 3 datasets. Among them, lncRNAs AFAP1-AS1, LINC00467, LINC00511, LINC00673, and LINC01207 were significantly up-regulated, and LINC00472, LINC00968, and TBX5-AS1 were significantly down-regulated in female lung cancer tissue in this study. Sui et al. analyzed the RNA sequencing data in 521 lung adenocarcinoma tissues and 49 non-tumor lung tissues from The Cancer Genome Atlas (TCGA) (), showing that AFAP1-AS1 was over-expressed and LINC00472 were under-expressed in lung cancer tissues compared to normal lung tissues [21]. Recently, some studies reported that AFAP1-AS1, LINC00511, LINC00673, and LINC01207 are up-regulated in lung cancer tissues, and promote cell proliferation, invasion, and metastasis [22-27]. Similarly, these lncRNAs were also involved in tumorigenesis and tumor progression and prognosis of other cancers [28-34]. Together, these results suggest that these aberrantly expressed lncRNAs play key roles in tumorigenesis and development of non-smoking female lung carcinoma. Furthermore, some other differentially expressed lncRNAs also related to the pathogenesis of cancers, such as enhanced LINC00467 expression, can promote neuroblastoma cell survival and reduced cell apoptosis [35]. Although we found LINC00968 and TBX5-AS1 were decreased in tumor tissue, the biological functions of LINC00968 and TBX5-AS1 in lung cancer or other carcinomas remain unclear. We reasonably surmise that these aberrantly expressed lncRNAs play important roles in initiation and development of female lung cancer.

The lncRNAs-miRNAs-mRNAs ceRNA network included 38 aberrantly expressed mRNAs. Functional and pathway enrichment analysis showed that the deregulated lncRNAs were associated with protein binding, metal ion binding, regulation of transcription from RNA polymerase II promoter, nucleic acid binding, microRNAs in cancer, cell differentiation, and hippo signaling pathway. Some cancer-related genes were found in the ceRNA network, such as LATS2 [36], EZH2[37], RECK [38] and NR4A2 [39], which were also correlated to the initiation and development of female lung cancer.

Conclusions

We identified aberrantly expressed key lncRNAs from the GEO database and constructed the lncRNA-miRNA-mRNA ceRNA network in non-smoking female lung cancer. Our results provide novel clues to understand the molecular mechanism of the pathogenesis of non-smoking female lung cancer and detect potential prognostic and diagnostic biomarkers. However, further molecular biological experiments are needed to confirm our findings.