Integrated Analysis of lncRNA–Mediated ceRNA Network in Lung Adenocarcinoma

Background A growing body of evidence indicates that long non-coding RNAs (lncRNAs) can act as competitive endogenous RNAs (ceRNAs) to bind to microRNAs (miRNAs), thereby affecting and regulating the expression of target genes. The lncRNA–miRNA–mRNA ceRNA network has been theorized to play an indispensable role in many types of tumors. However, the role of the lncRNA-related ceRNA regulatory network in lung adenocarcinoma (LUAD) remains unclear. Methods We downloaded the RNAseq and miRNAseq data of LUAD from The Cancer Genome Atlas (TCGA) data portal and identified differentially expressed lncRNAs (DElncRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) between LUAD and corresponding paracancerous tissues by using the edgeR package of R software. We constructed the lncRNA–miRNA–mRNA ceRNA network by using Cytoscape (version 3.7.2) on the basis of the interaction generated from the miRcode, miRTarBase, miRDB, and TargetScan databases. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed with DAVID 6.8 bioinformatics resources and plotted by using the ggplot2 package in R. The effect of genes on LUAD prognosis was assessed by applying the survival package in R in accordance with the Kaplan–Meier curve. Results In total, 1645 DElncRNAs, 117 DEmiRNAs, and 2729 DEmRNAs were identified in LUAD. The LUAD-specific ceRNA network was composed of 157 nodes and 378 edges (329 DElncRNA–DEmiRNA interactions and 49 DEmiRNA–DEmRNA interactions). GO and KEGG pathway annotations suggested that the LUAD-specific ceRNA network was related to tumor-related molecular functions and pathways. Seven lncRNAs (DISC1-IT1, SYNPR-AS1, H19, LINC00460, LINC00518, DSCR10, and STEAP2-AS1), one miRNA (hsa-mir-31), and 16 mRNAs (ATAD2, OSCAR, KIF23, E2F7, PFKP, MCM4, CEP55, CBX2, CCNE1, CLSPN, CCNB1, CDC25A, EZH2, CHEK1, SLC7A11, and PBK) were revealed to be significantly correlated with overall survival. Conclusion In this study, we described the potential regulatory mechanism of the progression of LUAD. We proposed a new lncRNA–miRNA–mRNA ceRNA network that could help further explore the molecular mechanisms of LUAD.


INTRODUCTION
Lung cancer is one of the most frequent cancers worldwide (1). It is classified into small cell lung carcinoma (accounting for 15% of cases) and non-small cell lung carcinoma (NSCLC, accounting for 85% of cases) (2). NSCLC is further divided into three types: lung adenocarcinoma (LUAD), squamous cell carcinoma, and large cell carcinoma; LUAD is the main histological subtype of lung cancer and accounts for 40% of all cases (3). Despite advances in diagnosis and treatment in recent years, the overall survival rate of patients with LUAD remains unsatisfactory. Statistics show that the average 5-year survival rate is less than 20% (4). The insufficient understanding of the biological mechanism of LUAD limits the improvement of therapeutic effects. Therefore, further elucidating tumor pathogenesis and finding new biomarkers to improve prognosis are urgent.
Non-coding RNAs (ncRNAs), which constitute more than 90% of the RNA transcripts of the human genome, lack protein coding potential but are considered the key regulators of normal cell function and many diseases, including cancer (5). ncRNAs include microRNAs (miRNAs, 21-24 base pairs) and long noncoding RNAs (lncRNAs, longer than 200 base pairs) (6)(7)(8)(9). As endogenous gene expression inhibitors, miRNAs can bind to the 3 untranslated region of target RNAs to promote mRNA degradation or inhibit protein translation (10,11). At present, miRNAs are widely accepted as oncogenes or tumor-suppressor genes that promote or inhibit the occurrence of tumors.
Long non-coding RNAs have been proven to be capable of regulating gene expression at multiple levels, namely, epigenetic, transcriptional, and post-transcriptional (12), and they are an indispensable participant in the development of cancer (13). In 2011, Salmena et al. (14) first proposed the competitive endogenous RNAs (ceRNA) hypothesis, which states that some RNAs, as ceRNAs, can regulate the expression of downstream mRNA by combining shared miRNAs. This hypothesis describes that ceRNAs can transform the function of target miRNAs by competing for the common binding sites of mRNAs on the target miRNAs. An increasing number of studies have shown that lncRNAs can act as ceRNAs to regulate target mRNA expression by competing for shared miRNAs. The schematic diagram of lncRNA-mediated ceRNA regulatory network is shown in Supplementary Figure 1. Several dysregulated expressed lncRNAs play an indispensable role in the development of tumors, showing their potential role in carcinogenesis and tumor inhibition (15,16). An increasing amount of research reported that lncRNAs, as ceRNAs, are involved in the development of a variety of tumors, such as lymphoma (17), colorectal cancer (18), and gastric cancer (19). However, to date, research on the role of lncRNA-related ceRNA regulator networks in LUAD is rare.
In our study, the genome-wide expression profiles of lncRNAs, miRNAs, and mRNAs in patients with LUAD were screened Abbreviations: ceRNA, competing endogenous RNA; DElncRNAs, differentially expressed lncRNAs; DEmiRNAs, differentially expressed miRNAs; DEmRNAs, differentially expressed mRNAs; DERNAs, differentially expressed RNAs; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma; MREs, miRNA response elements; ncRNAs, non-coding RNAs; TCGA, The Cancer Genome Atlas. through The Cancer Genome Atlas (TCGA) data portal. In addition, a LUAD-specific ceRNA network was established through comprehensive analysis. This network will help find new therapeutic targets and pathways for the treatment of patients and prolong the survival time of patients. Finally, we analyzed prognostic RNAs in the ceRNA network and found a biomarker that could be used to predict the survival of patients with LUAD.

Data Collection and Preprocessing
RNAseq data, miRNAseq data, and corresponding LUAD clinical data were downloaded from the TCGA data portal 1

DElncRNAs, DEmiRNAs, and DEmRNAs in LUAD
We identified mRNAs and lncRNAs by using the Ensembl database. Differentially expressed lncRNAs (DElncRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) between LUAD samples and corresponding paracancerous tissues were analyzed and normalized by using the edgeR package in R. FDR <1%, | logFC| > 2, and P < 0.05 were used as thresholds. We used gplots package in R to generate volcano plots and heatmaps.

Construction of a ceRNA Regulatory Network
The miRcode database 2 was used to match DElncRNAs and DEmiRNAs. MiRNA target genes were predicted on the basis of three databases: miRTarBase 3 , miRDB 4 , and TargetScan 5 . Subsequently, we integrated the interaction between DEmiRNAs and DElncRNAs or DEmRNAs to construct a ceRNA regulatory network. Cytoscape (version 3.7.2) was used to visualize the ceRNA network.

Functional Enrichment and
Protein-Protein Interaction Analysis DAVID 6.8 Bioinformatics Resources was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations. In addition, bubble charts were plotted by using the ggplot2 package of R software. P < 0.05 was considered a threshold. Protein-protein interaction (PPI) network was established with The STRING online search tool, with a combined score ≥0.4. Cytoscape v3.7.2 was used to visualized the PPI network.

Survival Analysis
Survival analysis was performed to assess the prognostic value of differentially expressed RNAs in LUAD by using the survival package in R. P < 0.05 was regarded as statistically significant.

Construction of the ceRNA Network in LUAD
A total of 105 DElncRNAs and 13 DEmiRNAs were paired into 329 DElncRNA-DEmiRNA interactions, whereas 13 DEmiRNAs and 39 DEmRNAs were matched to form 49 pairs of DEmiRNA-DEmRNA interactions. Figure 4A shows the number of DEmRNA in the ceRNA network. Finally, the LUAD-specific lncRNA-miRNA-mRNA ceRNA regulatory network, which contained 157 nodes and 378 edges (Figure 4B), was constructed.   Table 3.

Functional Enrichment Analysis and PPI Network Construction
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses revealed that the DEmRNAs involved in the ceRNA network were remarkably associated with 14 BPs, including DNA replication, cell cycle regulation, leukocyte migration, and G2 DNA damage checkpoint. The enrichment of MFs is mainly related to protein binding.    We found that the most enriched CC was the nucleus. ceRNA network-related genes were significantly enriched in three KEGG pathways, namely, the p53 signaling pathway, miRNAs in cancer, and cell cycle. The results of functional enrichment analysis are listed in Supplementary Table 3 and shown in Figure 5. Enrichment analysis suggested that the LUAD-specific ceRNA network might be involved in the tumor process by regulating these biological processes and pathways. In addition, to identify interactions of proteins translated from the mRNAs in the ceRNA network, The PPI network was constructed (Figure 6). We found that some genes with high combined score, including CHEK1, CDC25A, CCNE1, CCNB1, and MCM4 were mainly enriched in the "cell cycle" pathway.

DISCUSSION
Given the unsatisfactory survival rate and high mortality of LUAD, identifying specific biomarkers for the diagnosis and treatment of patients with LUAD is an urgent concern. An increasing number of studies have demonstrated that lncRNAs contribute to tumorigenesis and tumor progression through a variety of pathways. The ceRNA hypothesis proposes that lncRNAs, as ceRNAs, can regulate gene expression in LUAD by competing for shared miRNAs (20,21). For instance, lncRNA TTN-AS1 is overexpressed in LUAD compared with that in paracancerous tissues, and its overexpression can promote the malignant development of LUAD by regulating the miR-142-5p/cyclin-dependent kinase 5 signaling pathway (22). In addition, lncRNA LINC00355 and lncRNA LINC00466 can target CCNE1 and HOXA10 via sponging miR-195 and miR-144, respectively, thereby promoting the progress of LUAD (23,24). However, the genome-wide screening of the lncRNA-mediated ceRNA network in LUAD remains lacking.
In this study, we used LUAD expression data downloaded from the TCGA database to identify DElncRNAs, DEmiRNAs, and DEmRNAs between LUAD and corresponding paracancerous tissues. By integrating the interaction between DEmiRNAs and DEmRNAs or DElncRNAs, we constructed a LUAD-specific ceRNA network, which included 157 nodes and 378 edges. The results of GO enrichment analysis revealed that ceRNA-related RNAs were mainly involved in cancerrelated biological processes, such as DNA replication, cell cycle regulation, and leukocyte migration. In addition, the DEmRNAs involved in the ceRNA network were significantly enriched in three KEGG pathways, namely, the cell cycle, miRNA, and p53 signaling pathways. The enrichment results suggested that the lncRNA-miRNA-mRNA ceRNA network might regulate the biological processes and pathways of LUAD.
Although we constructed an LUAD-specific ceRNA regulatory network and screened for potential prognostic biomarkers, several limitations remain. First, this study was based merely on the gene expression data downloaded from the TCGA database. Prospective research involving different populations and centers with large patient sizes are needed to validate our results. Second, our research did not involve the clinical characteristics, such as TNM stage and gene mutation, of patients with LUAD. Third, further experimental research is needed to verify the potential biological mechanisms of these ceRNAs in LUAD in the future.

CONCLUSION
An LUAD-specific lncRNA-related ceRNA regulatory network was constructed by using bioinformatics methods. We also identified potential prognostic biomarkers, which provided novel insights into the tumorigenesis mechanism of LUAD. Further experimental studies are needed to validate this underlying biological regulatory mechanism in the future.

AUTHOR CONTRIBUTIONS
XW and ZS: design and initiation of the study, quality control of data, data analysis and interpretation, and manuscript preparation and editing. HZ and YW: data acquisition. ZY: study concept and design and initiation of the study. All authors: revision and approval of the final version of the manuscript.