Abstract

Background. Lung cancer is the most common cancer and the most common cause of cancer-related death worldwide. However, the molecular mechanism of its development is unclear. It is imperative to identify more novel biomarkers. Methods. Two datasets (GSE70880 and GSE113852) were downloaded from the Gene Expression Omnibus (GEO) database and used to identify the differentially expressed genes (DEGs) between lung cancer tissues and normal tissues. Then, we constructed a competing endogenous RNA (ceRNA) network and a protein-protein interaction (PPI) network and performed gene ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and survival analyses to identify potential biomarkers that are related to the diagnosis and prognosis of lung cancer. Results. A total of 41 lncRNAs and 805 mRNAs were differentially expressed in lung cancer. The ceRNA network contained four lncRNAs (CLDN10-AS1, SFTA1P, SRGAP3-AS2, and ADAMTS9-AS2), 21 miRNAs, and 48 mRNAs. Functional analyses revealed that the genes in the ceRNA network were mainly enriched in cell migration, transmembrane receptor, and protein kinase activity. mRNAs DLGAP5, E2F7, MCM7, RACGAP1, and RRM2 had the highest connectivity in the PPI network. Immunohistochemistry (IHC) demonstrated that mRNAs DLGAP5, MCM7, RACGAP1, and RRM2 were upregulated in lung adenocarcinoma (LUAD). Survival analyses showed that lncRNAs CLDN10-AS1, SFTA1P, and ADAMTS9-AS2 were associated with the prognosis of LUAD. Conclusion. lncRNAs CLDN10-AS1, SFTA1P, and ADAMTS9-AS2 might be the biomarkers of LUAD. For the first time, we confirmed the important role of lncRNA CLDN10-AS1 in LUAD.

1. Introduction

Lung cancer has the highest incidence of all cancers (11.6% of the total cases) and the highest death rate (18.4% of the total cancer deaths) [1]. It contains two main types, small cell (SCLC, 15% cases) and non-small-cell lung cancer (NSCLC, 85% cases). Lung adenocarcinoma (LUAD) and squamous cell carcinoma (LUSC) are the main histological subtypes of NSCLC [2]. The underlying biological and molecular mechanisms of lung cancer are gradually being understood over the past decades [3]. The molecular biomarkers may improve the early lung cancer detection [4]. Furthermore, traditional chemotherapy that is based on histopathology is replaced by the individualized precision treatment which is based on carcinogenic factors [5].

Protein-coding genes are widely studied in lung cancer over the past decades, but the human genome transcribes less than 2% of the protein-coding genes. 85% are noncoding RNAs, including long noncoding RNAs (lncRNAs) [6]. lncRNAs account for a large part of the human genome [7], and they were once considered to be the insignificant “noise” in genome’s repertoire of non-protein-coding transcripts. However, recent studies have revealed the roles of lncRNAs in many biological processes, including transcriptional regulation and cell differentiation [8,9]. It is well recognized that dysregulation of lncRNAs plays an important role in cancer, including lung cancer [10,11]. Nevertheless, a large number of lncRNAs are still unexplored [12]. Therefore, it is imperative to recognize more lncRNAs as biomarkers of lung cancer for better diagnosis, therapy, and prediction of the prognosis.

This study aimed to explore more biomarkers of lung cancer via integrated bioinformatics analysis. Gene Expression Omnibus (GEO) is a NCBI’s publicly available genomics database (https://www.ncbi.nlm.nih.gov/gds/), which provides us a large amount of genomic data about lung cancer. Two datasets (GSE70880 and GSE113852) were downloaded from GEO and used to identify the differentially expressed genes (DEGs) between lung cancer tissues and normal tissues. Then, we constructed a competing endogenous RNA (ceRNA) network and a protein-protein interaction (PPI) network and performed gene ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and survival analyses to identify potential biomarkers that are related to the diagnosis and prognosis of lung cancer. We verified the expression differences and measured the diagnostic roles of lncRNAs through The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/).

2. Materials and Methods

2.1. Gene Expression Microarray Datasets

We downloaded gene expression microarray datasets GSE70880 and GSE113852 from the GEO database. The criteria they met were as follows: (1) the studies were about lung cancer, (2) tissue samples included the tumor and corresponding adjacent tissues or normal tissues, (3) the number of total samples was not less than 40, and (4) the datasets must include lncRNAs and mRNAs.

2.2. Integrated Analysis of Microarray Datasets

We merged the two datasets not only to increase the sample size but also to facilitate subsequent analyses. Based on different platforms, days, environments, and people, the samples had heterogeneity and potential variables, which might lead to a bias. Therefore, we batch-normalized the merged dataset by suing limma and sva packages in software R (v3.6.1). Next, we performed gene differential analysis between tumor and normal tissues by limma package. |LogFC| > 1 and adjusted value <0.05 (the correction for value was done by Benjamini–Hochberg method) were considered statistically significant for the DEGs.

2.3. Construction of ceRNA Network

In order to find out the competing endogenous regulating network mediated by lncRNAs and miRNAs, the ceRNA network was constructed. CeRNA networks link the functions of protein-coding mRNAs with that of noncoding RNAs such as microRNA, long noncoding RNA, pseudogenic RNA, and circular RNA [12]. The construction of the ceRNA network included two steps. (1) Interactions between differentially expressed lncRNAs (DElncRNAs) and miRNAs were predicted by the miRcode database (http://www.mircode.org/) [13]. lncRNAs could competitively bind to the shared miRNAs. The upregulation and downregulation of one lncRNA result in more and less sequestrated copies of shared miRNAs, respectively [14]. (2) miRNA-mRNA interactions were predicted by miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/) [15], TargetScan (http://www.targetscan.org/) [16], and miRDB (http://www.mirdb.org/miRDB/) [17]. The interactions should match all the three databases. We selected the mRNAs that were differentially expressed in our study from the target mRNAs of miRNAs for the construction of the ceRNA network. We used the software Cytoscape (v3.7.1) to visualize the network.

2.4. Functional Enrichment Analysis and PPI Network

Functional enrichment analysis includes GO and KEGG analysis. Genes in the ceRNA network were analyzed by GO and KEGG analysis by package clusterprofiler in software R. In GO analysis, functions were divided into biological processes (BP), molecular functions (MF), and cellular component (CC). PPI network was constructed by using STRING (http://string‐db.org/cgi/input.pl).

2.5. Online Survival Analyses in LUAD and LUSC

To confirm whether lncRNA is associated with the survival of LUAD and LUSC was one of our research contents. Due to lack of clinical information, we performed Kaplan–Meier (KM) analyses of lncRNAs in the ceRNA network and genes in the PPI network through the Kaplan–Meier Plotter (http://kmplot.com/analysis/), respectively [18]. We only performed the survival analyses of these lncRNAs and mRNAs in LUAD and LUSC. Log rank was considered statistically significant. Genes related to the prognosis of LUAD or LUSC would be considered as the hub genes.

2.6. The Differences of Expression Levels and Diagnostic Roles of lncRNAs Selected from ceRNA Network

To confirm whether the expression levels of lncRNAs in the ceRNA network were different between normal and tumor tissues, we downloaded the gene expression profile from TCGA database and performed gene differential analysis of lncRNAs selected from the ceRNA network. means there is a statistical difference. In order to measure the diagnostic values of these lncRNAs, we used the data from TCGA database to perform receiver operating characteristic (ROC) curves by GraphPad Prism 7.0, and the values of cutoff, specificity, sensitivity, and area under the curve (AUC) were also calculated. The cutoff value is the best expression of a lncRNA for the diagnosis of lung cancer. indicated significant difference.

3. Results

3.1. Gene Expression Profile Data

There were two datasets and a total of 47 lung cancer tissue samples and 47 normal tissue samples in this study (Table 1). We merged the two datasets into one dataset and then batch-normalized it to reduce variability. 846 genes (321 genes were upregulated and 525 genes were downregulated) were confirmed as DEGs. We divided the DEGs into mRNA group (a total of 805 mRNAs, 307 mRNAs were upregulated and 498 mRNAs were downregulated) and lncRNA group (a total of 41 lncRNAs, 14 lncRNAs were upregulated and 27 lncRNAs were downregulated) (Table 2). Volcano plot showed the results of gene differential analysis (Figure 1). Heatmaps showed gene expression changes (Figure 2).

3.2. ceRNA Network

A total of 4 lncRNAs (CLDN10-AS1, SFTA1P, SRGAP3-AS2, and ADAMTS9-AS2), 21 miRNAs, and 48 mRNAs were involved in the ceRNA network (Figure 3). In tumor tissues, CLDN10-AS1 was upregulated, while SFTA1P, SRGAP3-AS2, and ADAMTS9-AS2 were downregulated. The number of nodes and lines in the network was 73 and 85, respectively. Among the 4 lncRNAs, the number of miRNAs that connected with ADAMTS9-AS2 was larger than that of other lncRNAs, which indicated the important role of ADAMTS9-AS2. miRNAs miR-125a-5p and miR-125b-5p were only connected with lncRNA CLDN10-AS1, which showed the special role of CLDN10-AS1 among the four lncRNAs.

3.3. GO and KEGG Pathway Analysis

GO terms and KEGG pathway analyses were performed by package clusterProfiler in software R. The results showed that biological processes were mainly associated with epithelial cell migration, regulation of epithelial cell migration, ameboidal-type cell migration, endothelial cell migration, endothelium development, and so on, while molecular functions were associated with carbohydrate binding, transmembrane receptor protein serine/threonine kinase activity, transmembrane receptor protein tyrosine phosphatase activity, transmembrane receptor protein phosphatase activity, carboxylic acid binding, and organic acid binding. Cellular component gathered in the basal part of cell, MCM complex, pore complex, cell-cell adherens junction, and basal plasma membrane (Figure 4). KEGG pathway analysis revealed the genes were only enriched in cell adhesion molecules (CAMs) and axon guidance (Figure 5(a)).

3.4. Protein-Protein Interaction

We uploaded the mRNAs in the ceRNA network to STRING to construct a PPI network (Figure 5(b)). There were a total of 25 nodes and 42 edges in the network. PPI enrichment value was . We calculated the number of genes connected to each gene in the PPI network. Results revealed that five genes had the highest connectivity (DLGAP5, E2F7, MCM7, RACGAP1, and RRM2, and the number of genes connected with each of the five genes was seven). All of them were upregulated. Immunohistochemistry (IHC) in The Human Protein Atlas (THPA) (https://www.proteinatlas.org/) database demonstrated that DLGAP5, MCM7, RACGAP1, and RRM2 were upregulated in LUAD (Figure 6).

3.5. Survival Analyses in LUAD and LUSC

We uploaded the four lncRNAs and five mRNAs to the Kaplan–Meier Plotter (http://kmplot.com/analysis/) to perform survival analyses in LUAD and LUSC, respectively (Figures 7(a) and 7(b)). The results showed that lncRNAs CLDN10-AS1 (logrank ), SFTA1P (logrank ), and ADAMTS9-AS2 (logrank ) were associated with prognosis of LUAD. Upregulation of CLDN10-AS1 predicted poor prognosis, while upregulation of the SFTA1P and ADAMTS9-AS2 was related to good prognosis. mRNAs DLGAP5 (logrank ), E2F7 (logrank ), MCM7 (logrank ), RACGAP1 (logrank ), and RRM2 (logrank ) were associated with prognosis of LUAD. Their high expression related to poor prognosis. No lncRNAs and mRNAs were associated with the prognosis of LUSC. DLGAP5, E2F7, MCM7, RACGAP1, and RRM2 were considered as hub genes. The heatmap showed the expression changes of the three lncRNAs and five mRNAs in the merged dataset (Figure 8).

Furthermore, we performed survival analyses of the three lncRNAs and five mRNAs in the cohorts of stage I LUAD patients, stage II LUAD patients, LUAD patients with smoking history, and LUAD patients without smoking history. The results showed that lncRNAs SFTA1P (logrank ) and ADAMTS9-AS2 (logrank ) and mRNAs RACGAP1 (logrank ) and RRM2 (logrank ) were associated with the prognosis of stage I LUAD patients (Figure 7(c)); only mRNA E2F7 (logrank ) was associated with the prognosis of stage II LUAD patients (Figure 7(d)); lncRNAs SFTA1P (logrank ), ADAMTS9-AS2 (logrank ), and mRNA RACGAP1 (logrank ) were associated with the prognosis of LUAD patients with smoking history (Figure 7(e)); lncRNA CLDN10-AS1 (logrank ), ADAMTS9-AS2 (logrank ), and mRNA DLGAP5 (logrank ) were associated with the prognosis of LUAD patients without smoking history (Figure 7(f)).

3.6. Expression Levels and Diagnostic Values of Four lncRNAs in LUAD Confirmed through TCGA

A total of 551 samples (54 normal samples and 497 tumor samples) of LUAD were downloaded from the TCGA database. After gene differential analysis, we observed that, in tumor samples, CLDN10-AS1 was upregulated, while SFTA1P, SRGAP3-AS2, and ADAMTS9-AS2 were downregulated (Figure 9). The results were consistent with our previous conclusions. Furthermore, the ROC curves showed that all of the four lncRNAs had diagnostic values in LUAD (Table 3, Figure 10).

4. Discussion

In the present study, we investigated the gene expression patterns of lncRNAs and mRNAs in lung cancer cells and corresponding adjacent tissues or normal tissues. We found that 525 genes were downregulated and 321 genes were upregulated in lung cancer cells. The ceRNA network revealed the correlation among lncRNAs, miRNAs, and mRNAs. Four lncRNAs (CLDN10-AS1, SFTA1P, SRGAP3-AS2, and ADAMTS9-AS2) were involved in ceRNA network and survival analyses showed that CLDN10-AS1, SFTA1P, and ADAMTS9-AS2 were associated with prognosis of LUAD. The PPI network showed the interaction among these mRNAs in ceRNA network and five mRNAs (DLGAP5, E2F7, MCM7, RACGAP1, and RRM2) had the highest connectivity. Survival analyses also revealed the relationship between the five mRNAs and the prognosis of LUAD, and they were confirmed as hub genes of LUAD finally. However, survival analyses indicated that these lncRNAs and mRNAs were not related to the prognosis of LUSC. To increase the credibility of our results, we verified the expression levels and measured the prognostic values of the four lncRNAs in LUAD through TCGA database. The results confirmed our conclusions and revealed that all the four lncRNAs were associated with the diagnosis of LUAD. lncRNAs CLDN10-AS1, SFTA1P, and ADAMTS9-AS2, and their related mRNAs DLGAP5, E2F7, MCM7, RACGAP1, and RRM2 might be biomarkers of LUAD.

Previous studies have explored the functions of lncRNAs as diagnostic and prognostic biomarkers [6]. HOTAIR is a famous cancer-related lncRNA and it is highly expressed in NSCLC and SCLC [10]. MALAT1 is another important lncRNA, and in patients with non-small-cell lung cancer, it is significantly related to metastasis potential and poor prognosis [19, 20]. Schmidt et al. indicated that MALAT1 could be considered as an independent prognostic parameter for both LUAD and LUSC [19, 20, 22]. Due to the high expression in LUAD, CCAT2 is a potential diagnostic biomarker for LUAD [22]. lncRNAs have also been studied as potential drug targets [6]. HOTAIR might be a therapeutic target in NSCLC because of its role in the chemoresistance to cisplatin [6]. Liu et al. indicated that MEG3 might be a potential therapeutic target in lung cancer, for tumor cells would be sensitive to cisplatin when MEG3 was overexpressed in A549 cells [23].

lncRNA SFTA1P had been reported previously. Huang et al. demonstrated that SFTA1P and CASC2 were associated with the regulation and development of LUSC and could be used as prognostic and predictive indicators of LUSC via integrated bioinformatics analysis [24]. Zhang et al. reported that the downregulation of SFTA1P affected LUAD patients’ survival time but had no influence on LUSC patients. SFTA1P exerted tumor inhibition in LUAD [25]. There were no experiments on the mechanism of SFTA1P in lung cancer so far. However, this opens up many avenues of study to pursue on this topic.

The number of miRNAs that connected with lncRNA ADAMTS9-AS2 was the largest among the four lncRNAs, suggesting it was a critical lncRNA. Studies about lncRNA ADAMTS9-AS2 in lung cancer were rare. Liu et al. demonstrated that lncRNA ADAMTS9‐AS2 was lowly expressed in lung cancer tissues by qRT-PCR. In their study, high expression of lncRNA ADAMTS9-AS2 reduced proliferation ability and inhibited migration and elevated their apoptosis rate. They also verified the relationship among lncRNA ADAMTS9-AS2, mRNA TGFBR3, and miRNA miR-223-3p. ADAMTS9-AS2 increased TGFBR3 expression, but miR-223-3p decreased both of them. miR-223-3p targeted TGFBR3 to enhance the ability of proliferation, migration, and invasion of lung cancer. They came to a conclusion that DAMTS9-AS2, TGFBR3, and miR-223-3p might provide potential therapeutic targets in lung cancer [26].

miRNAs miR-125a-5p and miR-125b-5p only connected with lncRNA CLDN10-AS1, indicating lncRNA CLDN10-AS1 had a different binding trend among the four lncRNAs. CLDN10-AS1 might play a role in the development and progression of lung cancer, by regulating the G1/S transition of mitotic cell cycle through miR-125b-5p/PPAT and by regulating endothelium development, angiogenesis, and cell-cell adherens junction through miR-125b-5p/CDH5, respectively. Furthermore, CLDN10-AS1 was associated with the prognosis of LUAD (logrank ), indicating the potential functions in the prognosis of LUAD patients. It is noteworthy that CLDN10-AS1 has not been reported in lung cancer. For the first time, we confirmed the crucial role of lncRNA CLDN10-AS1 in LUAD. It might be related to the diagnosis, ability of migration, and prognosis of LUAD.

Limitations existed in our study. (1) We only included two datasets. Although there were some other datasets, their sample size was too small. In order to increase the quality of our study, we excluded datasets with a sample size less than 40. (2) The two datasets are noncoding mRNA and lncRNA microarray, so we used the probe reannotation method to annotate gene symbol, which might drop some genes due to the failure of matching the probes. (3) Neither dataset contained clinical information. We could only perform survival analyses through the Kaplan–Meier Plotter (http://kmplot.com/analysis/). (4) We did not validate our results by experiment, giving us a direction for future research.

5. Conclusion

In our study, we identified the differentially expressed lncRNAs CLDN10-AS1, SFTA1P, SRGAP3-AS2, and ADAMTS9-AS2 by analyzing gene expression profiles from GEO. Among them, CLDN10-AS1, SFTA1P, and ADAMTS9-AS2 and their related mRNAs DLGAP5, E2F7, MCM7, RACGAP1, and RRM2 were associated with the prognosis of LUAD, suggesting they were more critical in LUAD. What is more, the four lncRNAs had diagnostic values in LUAD. lncRNAs CLDN10-AS1, SFTA1P, and ADAMTS9-AS2 and their related mRNAs DLGAP5, E2F7, MCM7, RACGAP1, and RRM2 might be biomarkers of LUAD. For the first time, we confirmed the important role of lncRNA CLDN10-AS1 in LUAD. It might be related to the diagnosis, ability of migration, and prognosis of LUAD.

Data Availability

All data can be obtained from GEO (https://www.ncbi.nlm.nih.gov/gds/), TCGA (https://portal.gdc.cancer.gov/), and THPA (https://www.proteinatlas.org/).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Donghui Jin designed the study, analyzed the data, and wrote the manuscript. Yuxuan Song designed the study and analyzed the data. Yuan Chen revised the paper critically for important intellectual content. Peng Zhang designed the study and revised the paper critically for important intellectual content.

Acknowledgments

The authors gratefully acknowledge contributions from the GEO, TCGA, and THPA database.