- Open Access
Identification of a novel prognosis-associated ceRNA network in lung adenocarcinoma via bioinformatics analysis
BioMedical Engineering OnLine volume 20, Article number: 117 (2021)
Lung adenocarcinoma (LUAD) is the most common subtype of nonsmall-cell lung cancer (NSCLC) and has a high incidence rate and mortality. The survival of LUAD patients has increased with the development of targeted therapeutics, but the prognosis of these patients is still poor. Long noncoding RNAs (lncRNAs) play an important role in the occurrence and development of LUAD. The purpose of this study was to identify novel abnormally regulated lncRNA–microRNA (miRNA)–messenger RNA (mRNA) competing endogenous RNA (ceRNA) networks that may suggest new therapeutic targets for LUAD or relate to LUAD prognosis.
We used the SBC human ceRNA array V1.0 to screen for differentially expressed (DE) lncRNAs and mRNAs in four paired LUAD samples. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to annotate the DE lncRNAs and mRNAs. R bioinformatics packages, The Cancer Genome Atlas (TCGA) LUAD database, and Kaplan–Meier (KM) survival analysis tools were used to validate the microarray data and construct the lncRNA–miRNA–mRNA ceRNA regulatory network. Then, quantitative real-time PCR (qRT-PCR) was used to validate the DE lncRNAs in 7 LUAD cell lines.
A total of 2819 DE lncRNAs and 2396 DE mRNAs (P < 0.05 and fold change ≥ 2 or ≤ 0.5) were identified in four paired LUAD tissue samples. In total, 255 of the DE lncRNAs were also identified in TCGA. The GO and KEGG analysis results suggested that the DE genes were most enriched in angiogenesis and cell proliferation, and were closely related to human cancers. Moreover, the differential expression of ENST00000609697, ENST00000602992, and NR_024321 was consistent with the microarray data, as determined by qRT-PCR validation in 7 LUAD cell lines; however, only ENST00000609697 was associated with the overall survival of LUAD patients (log-rank P = 0.029). Finally, through analysis of ENST00000609697 target genes, we identified the ENST00000609697–hsa-miR-6791-5p–RASL12 ceRNA network, which may play a tumor-suppressive role in LUAD.
ENST00000609697 was abnormally expressed in LUAD. Furthermore, downregulation of ENST00000609697 and its target gene RASL12 was associated with poor prognosis in LUAD. The ENST00000609697–hsa-miR-6791-5p–RASL12 axis may play a tumor-suppressive role. These results suggest new potential prognostic and therapeutic biomarkers for LUAD.
Lung cancer is the second most commonly diagnosed cancer worldwide, with 2.21 million new cases annually, and is the most common cause of cancer death (1.79 million deaths annually) . Approximately 85% of lung cancer cases are nonsmall-cell lung cancer (NSCLC), and lung adenocarcinoma (LUAD) is currently the most common histological subtype of NSCLC [2,3,4]. Patients with late-stage disease at diagnosis have a poor prognosis. However, the occurrence, development, and prognosis of tumors are not only related to pathological type and clinical stage, but also closely related to abnormal gene expression in tumor cells  Early prevention and continuous improvements in targeted drugs support the clinical translation of the lung cancer treatment model, prolong the progression-free survival (PFS) and overall survival (OS) of patients, and improve their prognosis . However, the prognosis of LUAD is still poor due to local recurrence or distant metastasis. Therefore, there is an urgent need to predict and explore more biomarkers for early diagnosis and therapeutic targets in LUAD.
Long noncoding RNAs (lncRNAs) are noncoding RNA (ncRNA) molecules of more than 200 nucleotides that are most commonly not translated into protein, but are crucial players in diverse cellular and physiological functions [7, 8]. Recently, with the development and application of high-throughput sequencing and gene chip technologies, researchers have found that lncRNAs play important roles in the occurrence and development of a variety of tumors . LncRNAs that are abnormally expressed in tumor tissues can not only be used as specific tumor biomarkers for early diagnosis and prognosis, but also directly interact with DNA, messenger RNA (mRNA), or protein to regulate chromatin modification or structure or to affect gene transcription, splicing, and translation . In general, lncRNAs can regulate a variety of physiological and pathological processes in tumor development, such as cell proliferation, differentiation, migration, and invasion; stem cell reprogramming; tumorigenesis; and drug resistance [11,12,13,14,15]. Competing endogenous RNAs (ceRNAs) are RNAs containing microRNA (miRNA) recognition elements (MREs)  that can regulate the expression of genes harboring the corresponding MRE or the expression of proteins by competitively binding to miRNAs . LncRNAs can also act as ceRNAs and function within lncRNA–miRNA–mRNA ceRNA networks [17, 18]. Perturbation of ceRNA networks may affect diseases and explain disease processes; thereby, presenting opportunities for new therapies . For instance, multiple studies have demonstrated that patients with various cancers with high HOTAIR expression exhibit higher lymphatic invasion and shorter survival times [19,20,21,22,23,24]. However, there is a need to discover new functional lncRNA–miRNA–mRNA ceRNA networks in LUAD. Therefore, research to identify more ceRNA networks related to LUAD diagnosis and prognosis is warranted.
In this study, we used gene chip technology to screen four paired LUAD tissue samples and analyzed and then predicted the multiple crucial functions of the identified DE lncRNAs and mRNAs. This analysis identified significantly more downregulated genes than upregulated genes. The data were also analyzed using R bioinformatics tools to construct the ceRNA regulatory network. Finally, we successfully screened the novel lncRNA ENST00000609697 and its target DE mRNA RASL12, which were negatively correlated with poor prognosis in LUAD. The results suggested that the ENST00000609697–hsa-miR-6791-5p–RASL12 ceRNA network may play a tumor-suppressive role and may contain new therapeutic targets and pathways related to LUAD survival.
Identification of DE lncRNAs and mRNAs in LUAD and adjacent tissues
In our study, microarray profiling was performed with four paired LUAD tissue samples (tumor and paracancerous tissues) to identify DE lncRNAs and mRNAs. An overview of this study is shown in Fig. 1. Overall, we identified 2819 lncRNAs and 2396 mRNAs with significant differential expression (P < 0.05 and fold change (FC) ≥ 2 or ≤ 0.5): 859 upregulated lncRNAs, 1960 downregulated lncRNAs, 757 upregulated mRNAs, and 1639 downregulated mRNAs. The hierarchical clustering heatmap showed the expression levels of the DE lncRNAs (Fig. 2A) and mRNAs (Fig. 2B) and the molecular signatures of these DE lncRNAs and mRNAs distinguish cancer from adjacent tissues. By volcano plot and scatter plot analyses evaluating the overall distribution of the two datasets, these DE RNAs were divided into up- and downregulated lncRNAs (Fig. 2C) and up- and downregulated mRNAs (Fig. 2D). The number of downregulated genes was significantly greater than that of upregulated genes. The top 20 significantly DE lncRNAs and mRNAs identified according to FC values are presented in Tables 1 and 2.
Validation of the DE lncRNAs via the Cancer Genome Atlas (TCGA) database
To verify the microarray data in a large cohort of clinical samples, we downloaded the TCGA LUAD database, which contains both gene expression and patient survival data for the screened cohort, and obtained 573 samples (including 514 LUAD tissue samples and 59 adjacent tissue samples). The clinical patient information is presented in Table 3. As shown in the hierarchical clustering heatmap (Fig. 2E) and the volcano plot (Fig. 2F), 1916 DE lncRNAs (|log2FC|> 1, P < 0.05); namely, 1271 upregulated and 645 downregulated lncRNAs, were identified. As shown in the Venn diagram (Fig. 2G), the intersection of the 2819 DE lncRNAs identified by microarray analysis with the 1916 DE lncRNAs identified by TCGA database analysis contained 255 overlapping DE lncRNAs (Additional file 1: 255 overlapping genes).
Annotation analyses of the DE lncRNAs and mRNAs
GO analysis was used to annotate gene functions and standardize the descriptions of the DE genes according to the biological process (BP), cellular component (CC), and molecular function (MF) categories. We analyzed the cis-regulated lncRNAs and found that most of the top 30 GO terms enriched by the upregulated and downregulated genes (i.e., DE lncRNAs and mRNAs) were in the BP and CC categories (Fig. 3A, B). The top three descriptive terms enriched by the DE lncRNAs were atomic septum development, structural molecule activity conferring elasticity, and embryonic digestive tract morphogenesis (Fig. 3C). However, condensed chromosome outer kinetochore, cell migration involved in heart development, and regulation of vasculogenesis were the top three descriptive terms enriched by the DE mRNAs (Fig. 3D). Moreover, all the DE lncRNAs and mRNAs were involved in angiogenesis and cell proliferation.
The KEGG database, which enables pathway analysis of DE genes to identify biological functions, is divided into the following six classifications: cellular processing, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. Comprehensive analysis of the KEGG classification results for the DE lncRNAs (Fig. 3G) and DE mRNAs (Fig. 3H) showed enrichment mainly in the terms signal transduction, immune system, and cancer: overview. Moreover, KEGG pathway enrichment analysis suggested that the DE genes were enriched mainly in vascular smooth muscle contraction, focal adhesion, and the TGF beta signaling pathway (Fig. 3E, F). Further analysis of the KEGG pathway term human diseases showed that these DE genes were closely related to small-cell lung cancer, NSCLC, melanoma, glioma, prostate cancer, thyroid cancer, colorectal cancer (CRC), and other tumors (Fig. 3G, H).
Analysis of lncRNA target genes
To further clarify the functional annotations of the DE genes, we determined the intersection of the target genes of the 1302 DE lncRNAs and the 2396 DE mRNAs, and the resulting 523 common DE genes were selected with Venn diagram software (Fig. 4A). GO analysis showed that these genes were also enriched in the CC and BP categories. Moreover, the main enriched terms were extracellular matrix, myosin complex, and cytoskeleton in the CC category (Fig. 4B); signal transduction in the BP category (Fig. 4C); and peptidase activity in the MF category (Fig. 4D). The KEGG analysis results showed that these genes were enriched mainly in the pathways focal adhesion, axon guidance, differentiated cardiomyopathy, and melanoma (Fig. 4E). These results suggested that the identified DE genes may play important roles in cell morphology, adhesion, intercellular connections, and signal transduction.
Candidate DE lncRNA validation in LUAD cell lines and OS analysis
We selected four candidate DE lncRNAs from the 255 overlapping genes: two downregulated genes (ENST00000609697 and ENST00000443224) and two upregulated genes (ENST00000602992 and NR_024321). To confirm the screening results, the expression of the 4 DE lncRNAs was validated in 7 LUAD cell lines and compared with that in the BEAS-2B cell line using qRT-PCR (Fig. 5A). The expression of ENST00000609697 and ENST00000443224 showed a significant decreasing trend in almost all the LUAD cell lines, consistent with the microarray data (P < 0.05), while ENST00000443224 was upregulated in H1993 cells (P < 0.05) (Fig. 5A). The significant increasing trend (P < 0.05) in ENST00000602992 and NR_024321 expression was also consistent with the microarray data, but the increasing trend in NR_024321 expression was not obvious in H2228 cells (Fig. 5B). Next, we downloaded gene expression data and patient follow‐up data from the TCGA dataset to elucidate whether these candidate genes are potential prognostic markers for LUAD. Through TCGA dataset analysis, we found that ENST00000609697 was downregulated (P < 0.001) (Fig. 5C) and was the only candidate gene related to the prognosis of LUAD (log-rank P = 0.029) (Fig. 5D). ENST00000602992 and NR_024321 were upregulated in the TCGA dataset (P < 0.001) (Additional file 2: Fig. S1A, B). However, ENST00000602992 was not associated with the prognosis of LUAD (P = 0.24) (Additional file 2: Fig. S1C), and NR_024321 upregulation was not positively correlated with good prognosis in LUAD (P = 0.018) (Additional file 2: Fig. S1D). However, the downregulation of ENST00000609697 was positively correlated with good prognosis in LUAD; and therefore, this gene was considered a candidate biomarker that may function as a tumor suppressor.
CeRNA regulatory network involving the DE lncRNAs
To further illustrate the potential interactions among the DE lncRNAs, miRNAs DE involved in LUAD, lncRNA–miRNA–mRNA ceRNA networks were constructed with a total of 188 DE lncRNAs, 444 DE miRNAs and 410 DE mRNAs (Additional file 2: Fig. S1E). Moreover, we found that most DE lncRNAs in the ceRNA regulatory network were downregulated. We selected the ceRNA network of the candidate gene ENST00000609697 for further analysis and found that it included 7 miRNAs (hsa-miR-3191-3p, hsa-miR-4731-5p, hsa-miR-598-5p, hsa-miR-6791-5p, hsa-miR-4292, hsa-miR-4446-3p, and hsa-miR-1827) and 20 DE mRNAs (COLGALT2, MYOCD, TNS1, RASL12, CNN1, etc.) (Fig. 6A). The enriched miRNA hsa-miR-4731-5p targeted most DE mRNAs in the ENST00000609697 ceRNA network, indicating that it may play a critical role in LUAD.
Functional and survival analyses considering the target DE mRNAs in the ENST00000609697 ceRNA network
We next conducted GO enrichment analysis of the 20 targeted DE mRNAs in three ontologies: BP, CC, and MF. The 30 GO terms most enriched by the 20 targeted DE mRNAs are shown in Fig. 6B. The most enriched GO terms in the BP, CC, and MF categories were smooth muscle cell differentiation, focal adhesion, and actin binding, respectively (Fig. 6B). Most DE mRNAs mapped to the BP category; thus, we generated a BP cnetplot that showed the DE mRNAs associated with the top 10 BP terms; this analysis identified, 4 DE mRNAs (CNN1, FLNC, FOXF1, and MYOCD) (Fig. 6C). FOXF1 and MYOCD are related to multiple biological processes, suggesting that they may be critical genes in LUAD. The BP emapplot showed the overlapping relationship between each pair of terms (Fig. 6D) and suggested that smooth muscle cell differentiation was a very important biological process. To screen ceRNA networks related to LUAD prognosis, we downloaded expression and survival data related to the 20 target DE mRNAs in the ENST00000609697 ceRNA network from the UCSC Xena database. RASL12 expression was downregulated in LUAD (P < 0.0001) (Fig. 6E), and this downregulation was positively correlated with good prognosis (P = 0.034) (Fig. 6F). These results suggested that the ENST00000609697–hsa-miR-6791-5p–RASL12 axis may play a tumor-suppressive role in LUAD.
Here, we identified 2819 DE lncRNAs and 2396 DE mRNAs, including 859 upregulated lncRNAs, 1960 downregulated lncRNAs, 757 upregulated mRNAs, and 1639 downregulated mRNAs. More genes were downregulated than upregulated, indicating that downregulated genes may play important roles in the biology of LUAD. To explore the potential mechanisms of the DE genes, we performed GO and KEGG analyses of the aberrantly expressed lncRNAs and mRNAs. GO analysis showed that the DE lncRNAs were enriched mainly in atomic septum development, structural molecule activity conferring elasticity, and embryonic digestive tract morphogenesis and that the DE mRNAs were enriched mainly in condensed chromosome outer kinetochore, cell migration involved in heart development, and regulation of vasculogenesis. However, all the DE lncRNAs DE mRNAs were involved in angiogenesis and cell proliferation. Abnormalities in these two processes are closely related to the occurrence and development of cancers [25, 26]. The KEGG classification results for the DE lncRNAs and mRNAs showed that they were enriched mainly in signal transduction, the immune system, and cancers. Moreover, KEGG pathway enrichment analysis suggested that these DE genes were enriched mainly in vascular smooth muscle contraction, focal adhesion, and the TGF beta signaling pathway, and were also closely related to small cell lung cancer, NSCLC, melanoma, glioma, prostate cancer, thyroid cancer, CRC, and other cancers. According to the previous study, the focal adhesion and the TGF beta signaling pathways play essential roles in cell proliferation, and dysregulation of these two pathways is closely associated with oncogenesis [27, 28]. In addition, to further verify whether the functional annotations of the DE lncRNAs and mRNAs were basically consistent, we used Venn diagram software to intersect the DE lncRNA target genes and DE mRNAs and found 523 overlapping genes, which were involved mainly in the extracellular matrix, myosin complex, cytoskeleton, and signal transduction, among other processes. Moreover, the KEGG analysis results showed that these genes were enriched mainly in the focal adhesion and melanoma pathways. The main functional annotations of the DE lncRNAs and mRNAs suggested that these genes may play important roles in cell morphology, adhesion, intercellular connections, and signal transduction and are highly related to cancer. Thus, they are worthy of further analysis and verification.
An increasing number of aberrantly expressed lncRNAs are being identified as novel key regulators of the development of multiple human cancers [29, 30]. Aberrantly expressed lncRNAs may serve as biomarkers or function as oncogenes or tumor suppressors ; however, most studies have focused on lncRNAs as oncogenes. For instance, Song et al. reported that the protein claudin-4 encoded by the CLDN4 gene was upregulated in gastric cancer and related to poor prognosis . The expression levels of the lncRNAs CCAT1 and CCAT2 were found to be significantly increased in CRC, and both lncRNAs were significantly correlated with poor relapse-free survival (RFS) and OS; these lncRNAs could thus be used independently or jointly as important prognostic biomarkers in CRC . In addition, lncTCF7 was found to be significantly overexpressed in liver tumor tissues and liver cancer stem cells (CSCs); this lncRNA can recruit the SWI/SNF complex to the TCF7 gene promoter to regulate its expression, thus activating the Wnt signaling pathway . In our study, the intersection of the 2819 DE lncRNAs identified by microarray analysis with the 1916 DE lncRNAs identified by TCGA database analysis revealed 255 overlapping DE lncRNAs: 161 downregulated and 94 upregulated. Then, we selected 4 candidate DE lncRNAs—2 downregulated (ENST00000609697 and ENST00000443224) and 2 upregulated (ENST00000602992 and NR_024321)—for validation in 7 LUAD cell lines. The expression of ENST00000609697, ENST00000602992, and NR_024321 was consistent with the microarray data. However, analysis of the relative expression levels of the candidate genes and the associations of these genes with patient survival in the TCGA dataset revealed that ENST00000609697 was downregulated and was the only candidate gene positively correlated with good prognosis in LUAD. Therefore, we considered ENST00000609697 a candidate gene that may be a novel tumor suppressor.
Recent studies have revealed that lncRNAs can act as ceRNAsby competitively binding to miRNAs to form lncRNA–miRNA–mRNA ceRNA networks and in turn play a critical role in the diagnosis, prognosis, and treatment of cancer [18, 34]. For example, lncRNA-KRTAP5-AS1 and lncRNA-TUBB2A can competitively bind miR-596 and miR-3620-3p as ceRNAs to promote CLDN4 expression, enhance cell proliferation and invasion, and promote epithelial–mesenchymal transition (EMT) . To identify the potential interactions among the DE lncRNAs, DE miRNAs, and DE mRNAs, we constructed lncRNA–miRNA–mRNA ceRNA networks involving a total of 188 DE lncRNAs, 444 DE miRNAs, and 410 DE mRNAs. Interestingly, most of the DE lncRNAs in the ceRNA regulatory network were downregulated. We then screened the ENST00000609697 ceRNA network, which was downregulated and positively correlated with good prognosis in LUAD. This network contained seven miRNAs (hsa-miR-3191-3p, hsa-miR-4731-5p, hsa-miR-598-5p, hsa-miR-6791-5p, hsa-miR-4292, hsa-miR-4446-3p, and hsa-miR-1827) and 20 DE mRNAs (COLGALT2, MYOCD, TNS1, RASL12, CNN1, etc.). We performed an in-depth analysis of the functions related to the ENST00000609697 ceRNA network and found that the most enriched GO terms in the BP, CC, and MF categories were smooth muscle cell differentiation, focal adhesion, and actin binding, respectively. Smooth muscle cell differentiation is very important for the stability and repair of the vascular system, and abnormalities in this biological process can directly or indirectly affect the growth, proliferation, and migration of tumor cells and the tumor immune microenvironment [35,36,37,38]. Focal adhesions are the center of cellular mechanical sensation and serve as bridges between integrin, the extracellular matrix and the cytoskeleton, which is correlated with the tumor microenvironment. Changes in signal transmission through focal adhesions of malignant cells are very important for tumor cell metastasis [38,39,40]. Actin binding-related proteins participate in the formation of the cytoskeleton and regulate cell adhesion and migration . The proliferation, migration, and invasion of tumor cells are dependent on proteins related to angiogenesis, focal adhesions, and actin binding. Therefore, the ENST00000609697 ceRNA network may play an important role in the tumor microenvironment of LUAD, and its functions are worth further exploration.
Subsequently, we downloaded expression and survival data for the 20 target DE mRNAs in the ENST00000609697 ceRNA network from the UCSC Xena database and found that RASL12 expression was downregulated in LUAD (P < 0.0001) and was positively correlated with good prognosis (P = 0.034). RASL12, a member of the RAS-like GTPase family, is localized in the cytoplasm . However, evidence that RASL12 functions as a small GTP-binding protein is lacking; In fact, studies have reported that RASL12 could be homologous to the RAS-like GTPases RERG, RASL11A, RASL11B, RASL10A and RASL10B, which play tumor-suppressive roles in human cancers [43,44,45]. In addition, a recent study reported that the tumor suppressor RASSF1 can form a complex with RASL12 and recruit RASL12 to microtubules . When combining these findings with our results, we inferred that RASL12 may be a tumor suppressor and that the ENST00000609697–hsa-miR-6791-5p–RASL12 axis may play a tumor-suppressive role in LUAD. More experiments should be performed to verify the role and regulatory mechanism of this axis.
Our study identified DE lncRNAs and mRNAs in LUAD tissue samples via microarray profiling and bioinformatics analysis approaches. Our results showed that downregulation of ENST00000609697 and its target gene RASL12 was associated with poor prognosis in LUAD. We identified a novel ceRNA network (ENST00000609697–hsa-miR-6791-5p–RASL12) that might play a tumor-suppressive role. These results might indicate potential molecular therapeutic targets and biomarkers for LUAD.
Materials and methods
Patient selection and tumor tissue collection
None of the patients with newly diagnosed LUAD received radiotherapy or chemotherapy before surgery. LUAD tissues—both tumor and paracancerous (> 5 cm from the tumor) tissues—were obtained from patients during thoracic surgery at the Affiliated Hospital of Hebei University. Thirty-four pairs of tissue samples were collected and pathologically confirmed. We selected four pairs of tissue samples for microarray screening: four tumor tissues (1-C, 5-C, 9-C, and 12-C) and four matched adjacent normal tissues (1-N, 5-N, 9-N, and 12-N). Patient characteristics are presented in Table 4, and the quality results of tissue RNA are present in Additional file 3.
All tissue samples were kept in liquid nitrogen prior to RNA extraction. We selected the four samples based on the time of sample collection, RNA quality, and sex ratio. The four samples were stored in liquid nitrogen for no more than 60 days, and the RNA quality was A1. The male to female ratio was 1:1. The study was approved by the Clinical Research Ethics Committee of the Affiliated Hospital of Hebei University, and all tissue samples were collected with written informed consent from the patients (Additional file 3).
The human bronchial epithelial cell line BEAS-2B and NSCLC cell lines (H1299, A549, H1975, HCC78, HCC827, H2228 and H1993) were purchased from the Typical Culture Preservation Commission Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and had no mycoplasma contamination. The BEAS-2B cell line was cultured in BEBM supplemented with bronchial epithelial cell growth factor (BEGM Kit, LONZA Corporation, USA). The NSCLC cell lines were cultured in RPMI 1640 medium (H1299, H1975, HCC78, HCC827, H2228, and H1993) or F12K medium (A549) supplemented with 10% fetal bovine serum and 1% penicillin–streptomycin at 37 °C in a humidified atmosphere containing 5% CO2.
Total RNA extraction and quantitative real-time PCR (qRT-PCR)
The total RNA was extracted from cells and tissues using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions. First-strand cDNA was synthesized using a Revert Aid First Strand cDNA Synthesis Kit (Roche, USA). qRT-PCR analysis of lncRNAs was performed in an Applied Biosystems 7500 Fast Real-Time PCR system (Applied Biosystems, USA) using One-Step SYBR PrimeScript (Roche, USA) according to the manufacturer’s instructions. The primer sequences are shown in Table 5.
The expression levels of lncRNAs were normalized to those of β-actin, and the relative lncRNA expression levels were calculated by the 2−ΔΔCt method.
LncRNA microarray analysis
The total RNA was amplified and labeled with a Low Input Quick Amp Labeling Kit, One-Color (Agilent Technologies, USA, Cat. #5190-2305), following the manufacturer’s instructions. Labeled cRNA was purified with a RNeasy Mini Kit (QIAGEN, GmBH, Germany, Cat. #74106). Each slide was hybridized in a hybridization oven with 1.65 μg of Cy3-labeled cRNA using a Gene Expression Hybridization Kit (Agilent Technologies, USA, Cat. #5188-5242) according to the manufacturer’s instructions. RNA samples from each group were analyzed using the SBC human ceRNA array V1.0 (4 × 180 K, Shanghai Biotechnology Co. Ltd., Shanghai, China), which contains approximately 68,423 lncRNA and 18,853 mRNA probes. LncRNA and mRNA probes were designed based on the latest genome version (human/grch38 (hg38)) with a probe length of 60 nt covering the GENCODE v21, Ensembl, UCSC, NONCODE, LNCipedia, lncRNAdb, and deepBase databases. Please refer to product information for details. The GeoPlatform number was GPL26192. Slides were scanned with an Agilent Microarray Scanner (Agilent Technologies, USA, Cat. #G2565CA) with default settings: dye channel, green; scan resolution, 3 μm; PMT, 100%; and range, 20 bits. The data were extracted with Feature Extraction software 10.7 (Agilent Technologies, Santa Clara, CA, USA).
Raw data were normalized with the quantile algorithm in the Limma package in R. Significantly DE transcripts were screened by FC ≥ 2 or ≤ − 2 and P ≤ 0.05. The screening principles to identify the four candidate genes were as follows: (1) DE lncRNAs were sorted by FC; (2) the difference in signal intensity between one pair of tumor tissue and normal control tissue was greater than 7 (downregulated DE lncRNAs, normal > 7; upregulated DE lncRNAs, cancer > 7) to facilitate later detection and verification; (3) the DE lncRNAs were located on chromosomes other than chrX and chrY; (4) the DE lncRNAs were shorter than 2 kb; and (5) exonic sense lncRNAs were excluded.
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses
We performed GO enrichment (http://www.geneontology.org/) and KEGG pathway analysis (https://www.kegg.jp/). GO and KEGG enrichment analyses were performed with the functions enrichGO and enrichKEGG in the R package clusterProfiler  with the following parameters: (1) the gene of interest was assigned the gene vector using the Entrez Gene ID for database annotation mapping; (2) the QrgDb object (annotation db) was set with org.Hs.eg.db, the default object for Homo sapiens; (3) the P value was calculated by a hypergeometric distribution:
where N is the number of genes in the background, M is the number of genes annotated in the background to the biological term, n indicates the length of the gene vector provided in parameter one, and k is the number of genes annotated to the corresponding biological term in set n; and 4) the P value was adjusted for multiple comparisons with the BH (Benjamini and Hochberg) method. The terms identified by this analysis were arranged in descending order according to the enrichment factor value, and the top 30 terms were considered.
Construction of the regulatory network
We selected the top correlated mRNA/lncRNA pairs in the normal and cancer datasets based on the correlation threshold of the 99th percentile of the corresponding overall correlation distribution in both cases. Then, we built up two regression models:
Seed match analysis
A perfect match at positions 2 to 7 of the 5ʹ-end of the mature miRNA sequence (6-mer miRNA seed) is the minimal pairing requirement considered predictive for miRNA target recognition. We used seed match analysis to restrict the above selected lncRNA/miRNA/mRNA triplets to those in which both the lncRNA and mRNA have at least one perfect 6-mer seed match for the shared miRNA.
By integrating the statistical analysis and seed match analysis results, we built the miRNA-mediated interaction (MMI) networks for both normal and cancer tissues. The nodes in these networks represent mRNAs and lncRNAs with highly correlated expression profiles, while the edges represent the miRNAs that mediate the interactions. The regulatory network was visualized using the Cytoscape tool  (3.9.0). The configurations (node color, position, and size) were chosen manually.
All data are presented as the mean ± standard error. Paired samples compared using a paired two-tailed Student’s t test. Multiple group comparisons were performed with an one-way ANOVA followed by Dunnett’s multiple comparison test. Statistical analyses were performed with GraphPad Prism 5 (GraphPad Software, USA). P values less than 0.05 were considered to indicate statistical significance.
Availability of data and materials
The datasets analyzed in this study are available from the corresponding authors on request.
Long noncoding RNAs
Competing endogenous RNA
Kyoto Encyclopedia of Genes and Genomes
Nonsmall-cell lung cancer
MicroRNA recognition elements
Acute myeloid leukemia
Pearson correlation coefficient
Ferlay J, Colombet M, Soerjomataram I, Parkin DM, Piñeros M, Znaor A, et al. Cancer statistics for the year 2020: an overview. Int J Cancer. 2020;149(4):778–89. https://doi.org/10.1002/ijc.33588.
Brahmer JR, Govindan R, Anders RA, Antonia SJ, Sagorsky S, Davies MJ, et al. The Society for Immunotherapy of Cancer consensus statement on immunotherapy for the treatment of non-small cell lung cancer (NSCLC). J Immunother Cancer. 2018;6(1):75. https://doi.org/10.1186/s40425-018-0382-2.
Liu Z, Sun D, Zhu Q, Liu X. The screening of immune-related biomarkers for prognosis of lung adenocarcinoma. Bioengineered. 2021;12(1):1273–85. https://doi.org/10.1080/21655979.2021.1911211.
Miyazawa T, Marushima H, Saji H, Kojima K, Hoshikawa M, Takagi M, et al. PD-L1 expression in non-small-cell lung cancer including various adenocarcinoma subtypes. Ann Thorac Cardiovasc Surg. 2019;25(1):1–9. https://doi.org/10.5761/atcs.oa.18-00163.
Zheng Q, Min S, Zhou Q. Identification of potential diagnostic and prognostic biomarkers for LUAD based on TCGA and GEO databases. Biosci Rep. 2021;41(6):BSR20204370. https://doi.org/10.1042/BSR20204370.
Zsákai L, Sipos A, Dobos J, Erős D, Szántai-Kis C, Bánhegyi P, et al. Targeted drug combination therapy design based on driver genes. Oncotarget. 2019;10(51):5255–66. https://doi.org/10.18632/oncotarget.26985.
Khalila AM, Guttmana M, Huartea M, Garbera M, Rajd A, Moralesa DR, et al. Many human large intergenic noncoding RNAs associate with chromatin-modifying complexes and affect gene expression. PNAS. 2009;106(28):11667–72. https://doi.org/10.1073/pnas.0904715106.
Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, et al. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009;458(7235):223–7. https://doi.org/10.1038/nature07672.
Chi Y, Wang D, Wang J, Yu W, Yang J. Long non-coding RNA in the pathogenesis of cancers. Cells. 2019;8(9):1015. https://doi.org/10.3390/cells8091015.
Bhan A, Mandal SS. LncRNA HOTAIR: a master regulator of chromatin dynamics and cancer. Biochim Biophys Acta. 2015;1856(1):151–64. https://doi.org/10.1016/j.bbcan.2015.07.001.
Bao Z, Yang Z, Huang Z, Zhou Y, Cui Q, Dong D. LncRNADisease 2.0: an updated database of long non-coding RNA-associated diseases. Nucleic Acids Res. 2019;47(D1):D1034–7. https://doi.org/10.1093/nar/gky905.
Geisler S, Coller J. RNA in unexpected places: long non-coding RNA functions in diverse cellular contexts. Nat Rev Mol Cell Biol. 2013;14(11):699–712. https://doi.org/10.1038/nrm3679.
Ponting CP, Oliver PL, Reik W. Evolution and functions of long noncoding RNAs. Cell. 2009;136(4):629–41. https://doi.org/10.1016/j.cell.2009.02.006.
Ye JR, Liu L, Zheng F. Long noncoding RNA bladder cancer associated transcript 1 promotes the proliferation, migration, and invasion of nonsmall cell lung cancer through sponging miR-144. DNA Cell Biol. 2017;36(10):845–52. https://doi.org/10.1089/dna.2017.3854.
Bhan A, Soleimani M, Mandal SS. Long noncoding RNA and cancer: a new paradigm. Cancer Res. 2017;77(15):3965–81. https://doi.org/10.1158/0008-5472.CAN-16-2634.
Fabian MR, Sonenberg N. The mechanics of miRNA-mediated gene silencing: a look under the hood of miRISC. Nat Struct Mol Biol. 2012;19(6):586–93. https://doi.org/10.1038/nsmb.2296.
Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011;146(3):353–8. https://doi.org/10.1016/j.cell.2011.07.014.
Zhao M, Feng J, Tang L. Competing endogenous RNAs in lung cancer. Cancer Biol Med. 2021;18(1):1–20. https://doi.org/10.20892/j.issn.2095-3941.2020.0203.
Ono H, Motoi N, Nagano H, Miyauchi E, Ushijima M, Matsuura M, et al. Long noncoding RNA HOTAIR is relevant to cellular proliferation, invasiveness, and clinical relapse in small-cell lung cancer. Cancer Med. 2014;3(3):632–42. https://doi.org/10.1002/cam4.220.
Liu X, Liu Z, Sun M, Liu J, Wang Z, De W. The long non-coding RNA HOTAIR indicates a poor prognosis and promotes metastasis in non-small cell lung cancer. Cancer Med. 2014;3(3):632–42. https://doi.org/10.1002/cam4.220.
Ding C, Cheng S, Yang Z, Lv Z, Xiao H, Du C, et al. Long non-coding RNA HOTAIR promotes cell migration and invasion via down-regulation of RNA binding motif protein 38 in hepatocellular carcinoma cells. Int J Mol Sci. 2014;15(3):4060–76. https://doi.org/10.3390/ijms15034060.
Heubach J, Monsior J, Deenen R, Niegisch G, Szarvas T, Niedworok C, et al. The long noncoding RNA HOTAIR has tissue and cell type-dependent effects on HOX gene expression and phenotype of urothelial cancer cells. Mol Cancer. 2015;14:108. https://doi.org/10.1186/s12943-015-0371-8.
Kim K, Jutooru I, Chadalapaka G, Johnson G, Frank J, Burghardt R, et al. HOTAIR is a negative prognostic factor and exhibits pro-oncogenic activity in pancreatic cancer. Oncogene. 2013;32(13):1616–25. https://doi.org/10.1038/onc.2012.193.
Kogo R, Shimamura T, Mimori K, Kawahara K, Imoto S, Sudo T, et al. Long noncoding RNA HOTAIR regulates polycomb-dependent chromatin modification and is associated with poor prognosis in colorectal cancers. Cancer Res. 2011;71(20):6320–6. https://doi.org/10.1158/0008-5472.CAN-11-1021.
Otto T, Sicinski P. Cell cycle proteins as promising targets in cancer therapy. Nat Rev Cancer. 2017;17(2):93–115. https://doi.org/10.1038/nrc.2016.138.
Zhao Z, Sun W, Guo Z, Zhang J, Yu H, Liu B. Mechanisms of lncRNA/microRNA interactions in angiogenesis. Life Sci. 2020;254:116900. https://doi.org/10.1016/j.lfs.2019.116900.
Gough NR, Xiang X, Mishra L. TGF-beta signaling in liver, pancreas, and gastrointestinal diseases and cancer. Gastroenterology. 2021. https://doi.org/10.1053/j.gastro.2021.04.064.
Mishra YG, Manavathi B. Focal adhesion dynamics in cellular function and disease. Cell Signal. 2021;85:110046. https://doi.org/10.1016/j.cellsig.2021.110046.
Zhang H, Chen Z, Wang X, Huang Z, He Z, Chen AY. Long non-coding RNA: a new player in cancer. J Hematol Oncol. 2013;6(37):1–7. https://doi.org/10.1186/1756-8722-6-37.
Kai-Xin L, Cheng C, Rui L, Zheng-Wei S, Wen-Wen T, Peng X. Roles of lncRNA MAGI2-AS3 in human cancers. Biomed Pharmacother. 2021;141:111812. https://doi.org/10.1016/j.biopha.2021.111812.
Song YX, Sun JX, Zhao JH, Yang YC, Shi JX, Wu ZH, et al. Non-coding RNAs participate in the regulatory network of CLDN4 via ceRNA mediated miRNA evasion. Nat Commun. 2021;12(1):3149. https://doi.org/10.1038/s41467-021-23211-y.
Ozawa T, Matsuyama T, Toiyama Y, Takahashi N, Ishikawa T, Uetake H, et al. CCAT1 and CCAT2 long noncoding RNAs, located within the 8q.24.21 “gene desert”, serve as important prognostic biomarkers in colorectal cancer. Ann Oncol. 2017;28(8):1882–8. https://doi.org/10.1093/annonc/mdx248.
Wang Y, He L, Du Y, Zhu P, Huang G, Luo J, et al. The long noncoding RNA lncTCF7 promotes self-renewal of human liver cancer stem cells through activation of Wnt signaling. Cell Stem Cell. 2015;16(4):413–25. https://doi.org/10.1016/j.stem.2015.03.003.
Abdollahzadeh R, Daraei A, Mansoori Y, Sepahvand M, Amoli MM, Tavakkoly-Bazzaz J. Competing endogenous RNA (ceRNA) cross talk and language in ceRNA regulatory networks: a new look at hallmarks of breast cancer. J Cell Physiol. 2019;234(7):10080–100. https://doi.org/10.1002/jcp.27941.
Thompson AM, Martin KA, Rzucidlo EM. Resveratrol induces vascular smooth muscle cell differentiation through stimulation of SirT1 and AMPK. PLoS ONE. 2014;9(1):e85495. https://doi.org/10.1371/journal.pone.0085495.
Sun HM, Chen XL, Chen XJ, Liu J, Ma L, Wu HY, et al. PALLD regulates phagocytosis by enabling timely actin polymerization and depolymerization. J Immunol. 2017;199(5):1817–26. https://doi.org/10.4049/jimmunol.1602018.
Yoshio T, Morita T, Kimura Y, Tsujii M, Hayashi N, Sobue K. Caldesmon suppresses cancer cell invasion by regulating podosome/invadopodium formation. FEBS Lett. 2007;581(20):3777–82. https://doi.org/10.1016/j.febslet.2007.06.073.
Liu Z, Liu X, Cai R, Liu M, Wang R. Identification of a tumor microenvironment-associated prognostic gene signature in bladder cancer by integrated bioinformatic analysis. Int J Clin Exp Pathol. 2021;14(5):551–66.
Branis J, Pataki C, Sporrer M, Gerum RC, Mainka A, Cermak V, et al. The role of focal adhesion anchoring domains of CAS in mechanotransduction. Sci Rep. 2017;7:46233. https://doi.org/10.1038/srep46233.
Paluch EK, Aspalter IM, Sixt M. Focal adhesion-independent cell migration. Annu Rev Cell Dev Biol. 2016;32:469–90. https://doi.org/10.1146/annurev-cellbio-111315-125341.
Zhou J, Kang X, An H, Lv Y, Liu X. The function and pathogenic mechanism of filamin A. Gene. 2021;784:145575. https://doi.org/10.1016/j.gene.2021.145575.
Ueda K, Fujiki K, Shirahige K, Gomez-Sanchez CE, Fujita T, Nangaku M, et al. Genome-wide analysis of murine renal distal convoluted tubular cells for the target genes of mineralocorticoid receptor. Biochem Biophys Res Commun. 2014;445(1):132–7. https://doi.org/10.1016/j.bbrc.2014.01.125.
Finlin BS, Gau CL, Murphy GA, Shao H, Kimel T, Seitz RS, et al. RERG is a novel ras-related, estrogen-regulated and growth-inhibitory gene in breast cancer. J Biol Chem. 2001;276(45):42259–67. https://doi.org/10.1074/jbc.M105888200.
Zhao W, Ma N, Wang S, Mo Y, Zhang Z, Huang G, et al. RERG suppresses cell proliferation, migration and angiogenesis through ERK/NF-kappaB signaling pathway in nasopharyngeal carcinoma. J Exp Clin Cancer Res. 2017;36(1):88. https://doi.org/10.1186/s13046-017-0554-9.
Zou H, Hu L, Li J, Zhan S, Cao K. Cloning and characterization of a novel small monomeric GTPase, RasL10B, with tumor suppressor potential. Biotechnol Lett. 2006;28(23):1901–8. https://doi.org/10.1007/s10529-006-9176-6.
Dhanaraman T, Singh S, Killoran RC, Singh A, Xingjian Xu, Shifman JM, et al. RASSF effectors couple diverse RAS subfamily GTPases to the Hippo pathway. Cell Biol. 2020;13(635):1–15. https://doi.org/10.1126/scisignal.abb4778.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.
We thank Shanghai Biotechnology Corporation and Dr. Mingjie Chen (NewCore Biotech Co., Ltd., Shanghai) for microarray profiling and bioinformatics analysis related to this study.
This study was supported by The Youth Scientific Research Fund Project of the Affiliated Hospital of Hebei University (Grant number 2017Q002), The Key Scientific Research Fund Project of the Affiliated Hospital of Hebei University (Grant number 2015Z005), The Hebei Province Government Foundation for Clinical Medical Talent Training and Basic Research Projects (Grant number 361007), and The Outstanding Young Scientific Research and Innovation Team of Hebei University (Grant number 605020521007).
Ethics approval and consent to participate
The study was approved by the Clinical Research Ethics Committee of the Affiliated Hospital of Hebei University, and all tissue samples were collected with written informed consent from the patients.
Consent for publication
All authors declare that there are no conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
255 overlapping genes.
(A) The relative expression of the candidate DE lncRNA ENST00000602992 in the TCGA dataset. (B) The relative expression of the candidate DE lncRNA NR_024321 in the TCGA dataset. (C) Kaplan–Meier (KM) survival analysis based on the candidate DE lncRNA ENST00000602992. (D) Kaplan–Meier (KM) survival analysis based on the candidate DE lncRNA NR_024321. x axis: overall survival (years); y axis: survival rate. Green and red represent the low and high DE lncRNA expression groups, respectively. (E) CeRNA and coexpression regulatory networks of DE lncRNAs. CeRNA network: green and red represent down- and upregulated lncRNAs, respectively; □ represents lncRNA; ○ represents mRNA; the line represents miRNA; and the size of the circle or square represents the ability of the gene to interact with other genes.
Tissue RNA quality.
About this article
Cite this article
Li, Y., Yu, X., Zhang, Y. et al. Identification of a novel prognosis-associated ceRNA network in lung adenocarcinoma via bioinformatics analysis. BioMed Eng OnLine 20, 117 (2021). https://doi.org/10.1186/s12938-021-00952-x
- Lung adenocarcinoma