Comprehensive Analysis of Competing Endogenous RNA (ceRNA) Network Based on RNAs Differentially Expressed in Lung Adenocarcinoma Using The Cancer Genome Atlas (TCGA) Database

Background

Lung adenocarcinoma is the leading cause of cancer deaths in the world, and its mortality rate is still rising [1]. The incidence and mortality of male lung cancer patients rank first in all malignant tumors, and female lung cancer patients account for the second place [2]. Non-small cell lung cancer accounts for about 85% of lung cancer [3] and adenocarcinoma is more common than squamous cell carcinoma [4].

The sequencing of the human genome shows that only a small fraction of the DNA could encode the protein in the entire genome, most of the DNA could not. Most of the transcription products from DNA of these noncoding proteins are noncoding RNAs (ncRNAs) and long noncoding RNAs (lncRNAs) [5]. LncRNA is a type of ncRNA with a transcript longer than 200 nucleotides and lncRNA, although it lacks protein coding ability, plays an important role in cell differentiation and cell metabolism and other related life activities [6]. In 2011, Salmena et al. first proposed the hypothesis of competing endogenous RNA (ceRNA), believed that messenger RNA (mRNA), lncRNA, and other transcription products can be use as ceRNAs to compete and combinate with microRNAs (miRNAs) through the same miRNA response element (MRE) to regulate the expression level of genes and achieve mutual regulation, thereby affecting the function of cells [7]. In addition, ceRNA can competitively bind miRNAs to regulate the number of free miRNAs to influence the translation and silencing effect of miRNAs on their target mRNAs, thereby regulating the expression of target genes [8]. Currently, many experiments have proven that lncRNA can serve as a competitive platform for miRNAs and mRNAs, which allow it to play an important role in the regulation of cell cycle and apoptosis in various malignant tumors such as breast, gastric, liver, lung, and colorectal cancer, thereby affecting tumor invasion and migration. Thus, lncRNA plays an important role in tumorigenesis [9–11], however, its function in the tumor mechanism is still unclear.

The Cancer Genome Atlas (TCGA) is a publicly funded project with the main goal to provide publicly available datasets to help improve diagnostic methods, treatment standards, and ultimately prevent cancer [12]. In this study, we obtained RNA expression profiles and clinical information of lung adenocarcinoma patients from TCGA database, and further constructed the lncRNA-miRNA-mRNA ceRNA network of lung adenocarcinoma. Important biological functions were identified by Gene Ontology (GO) enrichment analysis of differentially expressed mRNAs in ceRNA. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was used to analyze signal transduction and metabolic pathways. The Kaplan-Meier survival analysis was used to identify prognostic characteristics of mRNA, lncRNA, and miRNA involved in the ceRNA network, and to select potential prognostic factors for patients with lung adenocarcinoma.

Material and Methods

TCGA DATA COLLECTION:

Based on GDC (Genomic Data Commons; https://portal.gdc.cancer.gov/), 594 patients with lung adenocarcinoma mRNA, lncRNA, and miRNA expression data (including 535 cases with lung adenocarcinoma tissues and 59 cases of normal tissues) were downloaded and extracted from TCGA database (up to May 1, 2019) by using “TCGA biolinks” package of R software; the clinical data of patients were extracted. Gene expression levels were used for differential gene analysis and ceRNA network construction, and patient clinical information was used for survival analysis. This study was performed using TCGA publication guidelines. Since all patient data were from TCGA, no ethics approval was required.

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED RNAS:

MRNAs and lncRNAs were encoded according to Gencode Release v30 (GRCh38.p12) (https://www.gencodegenes.org/human/), which was converted Ensembl ID to gene symbol by using Perl language. To filter out all low expression RNAs, we set the file path to read mRNAs, lncRNAs, and miRNAs expression data to remove all average reading ≤2 using “edgeR” package of R software. Fold change (FC) ≥2 and P-value adjusted to false discovery rate (FDR) <0.01 were considered as 2 screening criteria for the analysis of differentially expressed genes of lung adenocarcinoma. All data were subjected to a “2-tailed test”. The volcano plots were visualized by using “gplots” package of R software.

CERNA NETWORK OF LNCRNA-MIRNA-MRNA:

Based on the screened differentially expressed mRNAs, lncRNAs, and miRNAs, the miRcode (http://www.mircode.org/) database was first used to predict the interaction between the differential lncRNA and the differential miRNA. Second, 4 databases: miRDB (http://www.mirdb.org/), mirTarBase (http://www.mirbase.org/), miRanda (http://www.microrna.org/microrna/home.do), and TargetScan (http://www.targetscan.org/) were used to retrieve and predict target mRNAs for differential miRNAs. Finally, the Venny plots were visualized by using the “VennDiagram” package of R software, which was then used to intersected with the extracted differentially expressed mRNA. The ceRNA network of lncRNA-miRNA-mRNA for lung adenocarcinoma was constructed by Cytoscape v3.7.1.

GO ENRICHMENT ANALYSIS AND KEGG PATHWAYS ANALYSIS:

To explore the potential biological mechanisms of mRNA in the network, a database for annotation, visualization and integrated discovery (DAVID) v6.8 (https://david.ncifcrf.gov/) were used to perform GO enrichment analysis and KEGG pathway analysis of target mRNAs involved in ceRNA networks. The relationship between GO terms and target genes were visualized by using “ggplot2” package of R software.

SURVIVAL AND PROGNOSIS ANALYSIS:

The effective survival status and survival time were extracted from the clinical data of lung adenocarcinoma patients in TCGA database. By established Lasso regression model and Cox regression model, the differential mRNAs, lncRNAs, and miRNAs related to survival in the network were screened. Survival analysis was performed using the “survival” software package of the R software (set P<0.05).

Results

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED GENES IN LUNG ADENOCARCINOMA:

The mRNAs, lncRNAs, and miRNAs expression data of lung adenocarcinoma were downloaded and extracted from TCGA database. Using the “edgeR” package of R software, we set criteria as fold change (FC) ≥2, FDR (adjusted P-value) <0.01 and average expression values greater than 2. After analysis, 2551 differentially expressed mRNAs were screened, including 2033 upregulated and 518 downregulated, 1359 differentially expressed lncRNAs, including 1185 upregulated and 174 downregulated, 99 differential expression miRNAs, including 78 upregulated and 21 downregulated. The volcano plots were visualized by using “gplots” package of R software (Figure 1A–1C).

CONSTRUCTION OF CERNA NETWORK:

To construct the ceRNA network of lncRNA-miRNA-mRNA, we search 100 differentially expressed lncRNAs through the miRcode database, 43 differential miRNAs corresponding to the 100 differential lncRNAs were found. The predictive target mRNAs for 43 differential miRNAs were retrieved and extracted from the 4 databases of TargetScan, miRDB, miRanda, mirTarBase, respectively. The Venny plot was visualized by using “VennDiagram” package of R software (Figure 2A). Finally, 294 target mRNAs were obtained though intersection with the extracted differentially expressed mRNAs (Figure 2B). Then, we further study 100 lncRNAs, 43 miRNAs, and 294 mRNAs which were involved in the ceRNA network, and 2 sub-networks were extracted. The ceRNA and 2 sub-networks were visualized by Cytoscape v3.7.1 (Figure 3A–3C).

ENRICHMENT ANALYSIS OF DIFFERENTIALLY EXPRESSED MRNA IN CERNA NETWORK:

In order to analyze the potential biology of the target mRNAs functions in the ceRNA network, we input the differential mRNAs and remove the background according to the settings to extract the GO enrichment analysis data (set FDR <0.05) by using the DAVID database. 5 enriched GO terms in the molecular function (MF) and biological process (BP) were identified (Figure 4, Table 1). Ten significant pathways identified by the KEGG pathway analysis (set P<0.05) (Figure 5, Table 2).

From the aforementioned, we found that the target mRNAs were significant enrichment in molecular function.

UNIVARIATE COX PROPORTIONAL HAZARD REGRESSION ANALYSIS:

To clarify the prognostic relevance of differential RNAs in lung adenocarcinoma, we performed a univariate Cox proportional hazard regression analysis of 100 lncRNAs, 49 miRNAs, and 294 mRNAs differentially expressed in the ceRNA network to initially screened for differential RNAs expression in lung adenocarcinoma. P<0.05 was set a significant correlation between the patient’s overall survival (OS). A total of 144 mRNAs, 29 lncRNAs, and 5 miRNAs were screened for the possible prognosis of lung adenocarcinoma.

BUILDING A LASSO LOGISTIC REGRESSION MODEL: For differential RNAs screened by univariate Cox proportional hazard regression analysis, the data was normalized using a Lasso logistic regression model and overfitted differential RNA models were removed. Before construct the model, it was necessary to further matrix the independent variables (expression data of each sample), set the dependent variables (survival state and survival time) to 2 categorical variables, and applied the Lasso logistic regression model to the independent variables. Variables and dependent variables (Figure 6A–6C).

From the Lasso logistic regression model, λ was a very important parameter in the Lasso regression. The number of independent variables included in the model could be adjusted according to the choice of different λ values. Therefore, determining the optimal λ value of the model was very important for the study. This study used a 10-fold cross-validation method to determine an optimal λ value (Figure 7A–7C).

This study was chosen the values of lambda.min which were 0.06632688, 0.01893901 and 0.00230649, respectively, corresponding to the dotted line on the left side of each figure in Figure 7; and 7 mRNAs, 18 lncRNAs, and 5 miRNAs were screened from the model finally. The variables entering the model at this time have a survival state and a survival time.

MULTI-FACTOR COX RISK REGRESSION ANALYSIS: There were 7 mRNAs, 18 lncRNAs, and 5 miRNAs selected by Lasso regression analysis. Multivariate Cox proportional hazard regression analysis was performed (P<0.05). The forest maps were visualized by “survminer” of R software (Figure 8A–8C).

There were 3 mRNAs (ANLN, IGFBP1, and TFAP2A), 7 lncRNA (AC015923.1, AC017048.3, FGF12-AS2, HOXA11-AS, LINC00211, LINC00319, and MED4-AS1) and 3 miRNAs (hsa-mir-142, hsa-mir-31, and hsa-mir-490) associated with the prognostic of lung adenocarcinoma. Significantly, has-mir-31 shows the greatest prognostic value related with the adverse effect of the prognostic of lung adenocarcinoma (P<0.001).

KAPLAN-MEIER SURVIVAL ANALYSIS: We performed a Kaplan-Meier survival analysis for the 3 mRNAs, 7 lncRNAs, and 3 miRNAs that showed significant prognosis in the forest map and plotted the survival curve by using “survival” package of R software. There were 3 mRNAs (ANLN, IGFBP1, and TFAP2A), 4 lncRNAs (AC015923.1, FGF12-AS2, LINC00211, and MED4-AS1) and 2 miRNAs (hsa-mir-31, and hsa-mir-490) found to be associated with the prognosis of lung adenocarcinoma (Figure 9A–9I).

Then, we further evaluated the risk score of each patient, and divided the patients into a high-risk score group and a low-risk score group with the corresponding median risk score as the boundary, and plotted the risk score curve (Figure 10A–10C). The results showed that the prognosis of the lung adenocarcinoma low-risk score group was significantly better than that of the high-risk score group (P<0.05).

Discussion

Lung adenocarcinoma is the most common histological type of lung cancer, and its morbidity and mortality continue to increase [13,14]. Lung adenocarcinoma treatment includes surgery, chemotherapy, and radiotherapy. Although considerable progress has been made in the treatment of lung adenocarcinoma, the overall treatment effect is still unsatisfying, and the prognosis of patients with lung adenocarcinoma remains poor. The overall 5-year survival rate is still below 15% [15]. Therefore, the exploration of biomarkers that can provide potentially effective target therapy has become a hot spot in current lung adenocarcinoma research [16]. In recent years, the continuous improvement of high-throughput transcriptome sequencing technology and bioinformatics analysis methods has opened up new ideas and provided new methods for studying lung adenocarcinoma [17].

The miRNA is a small molecule ncRNAs of about 22 nucleotides, and about one-third of all human genes are regulated by these RNAs [18]. The gene expression is often not regulated by a single miRNA, but the result of co-regulation by multiple miRNAs [19]. MiRNA specifically binds to target gene mRNA at specific binding sites and participates in the regulation of many important life activities (such as cell differentiation, proliferation, apoptosis, and hormone secretion), thereby inhibiting protein synthesis at the post-transcriptional level [20].

At present, many studies have shown that ceRNA is related to many important biological functions of cancer incidence and development [21]. In addition, a novel regulatory mechanism, the ceRNA hypothesis, which could be mediated through lncRNA, has been proposed [22]. LncRNA XIST has been shown to promote epithelial-to-mesenchymal transition (EMT) of retinoblastoma via sponge miR-101 [23]. Silent lncRNA HOXA-AS2 has been found to inhibit breast cancer cell proliferation and invasion via miR-520c-3p sponge [24]. With the further development and analysis of ceRNA networks, miRNAs, lncRNAs, or mRNAs cannot be considered as single elements in tumors; the expression of mRNAs changes with miRNAs, eventually leading to the occurrence of tumors [25]. Because of the complex mechanisms of cancer development, the role of a single gene or pathway is likely very limited. Therefore, systematic construction and analysis of ceRNA networks can provide more effective research directions for us to study potential cancer mechanisms.

In this study, we used R software to screen differentially expressed genes in lung adenocarcinoma from TCGA database, and further constructed the ceRNA network of lncRNA-miRNA-mRNA in lung adenocarcinoma. We also found that LINC00461 and has-mir-139 are 2 important network nodes in lung adenocarcinoma ceRNA network. Moreover, we predict the enrichment function and pathway of differential mRNAs involved in ceRNA networks by GO enrichment analysis and KEGG pathway analysis. The correlation between differential genes in the network and overall survival was further explored. We found 3 mRNAs (ANLN, IGFBP1, and TFAP2A) in differentially expressed genes, 4 lncRNAs (AC015923.1, FGF12-AS2, LINC00211, and MED4-AS1) and 2 miRNAs (has-mir-31 and has-mir-490) associated with the prognostic of lung adenocarcinoma. Significantly, has-mir-31 showed the greatest prognostic value related with the adverse effect of the prognostic of lung adenocarcinoma (P<0.001). Finally, we applied risk scores to predict the association between lung adenocarcinoma patients of different risk groups and OS; we found that lung adenocarcinoma patients in the low-risk score group had significantly better prognosis than those in the high-risk score group.

Studies have shown that hsa-miR-31 has been identified as a new prognostic biomarker for SCC in Chinese patients [26]. Li et al. found that miRNA-490 regulates lung cancer metastasis by targeting poly (c) binding protein 1, which confirmed that hsa-mir-490-3p may be a potential prognostic marker for lung cancer [27]. Xu et al. showed that the expression of ANLN in lung adenocarcinoma was related to the metastasis of cancer cells. It has also been reported that the expression of ANLN in lung adenocarcinoma is associated with metastasis of cancer cells [28]. Long et al. used immunohistochemical staining in lung adenocarcinoma patients to show high expression of ANLN in adenocarcinoma tissues, and found that the expression level of ANLN was associated with TNM stage and OS rate of lung adenocarcinoma patients, suggesting that ANLN is a potential biomarker for diagnosis and prognosis of lung adenocarcinoma [29]. TFAP2A is an AP-2 family transcription factor that binds the consensus sequence 5′-GCCNNNGGC-3′ and is overexpressed in lung adenocarcinoma. TFAP2A regulates the expression of multiple downstream target genes involved in human cancer development and progression. Liu et al. designed a migration rescue experiment and found that TFAP2A gene knockout could affect EMT and reduce cancer cell migration in lung adenocarcinoma [30]. Shi et al. in a study of TFAP2A found it to be highly expressed in various nasopharyngeal carcinoma cells and tumors, and that it was related to the expression of hypoxia-inducible factor-1α (HIF-1α). In addition, TFAP2A overexpression was positively correlated with tumor staging, local invasion, clinical progress, and poor prognosis in patients with nasopharyngeal carcinoma, which suggests that TFAP2A may be a potential prognostic marker for nasopharyngeal carcinoma [31].

Conclusions

In summary, by analyzing the expression data of lung adenocarcinoma in TCGA database, we found that 3 mRNAs, 4 lncRNAs, and 2 miRNAs were screened as potential prognostic biomarkers for lung adenocarcinoma. In addition, by constructing the lncRNA-miRNA-mRNA ceRNA network we provided new ideas for the competitive interaction between different RNAs in lung adenocarcinoma patients, and we identified LINC00461 and has-mir-139 as 2 important regulatory network nodes in lung adenocarcinoma ceRNA. These findings provided more insights into the regulatory mechanisms underlying lung adenocarcinoma.