# Construction and Analysis of Survival-Associated Competing Endogenous RNA Network in Lung Adenocarcinoma

Authors: Lixian Chen; Zhonglu Ren; Yongming Cai
Journal: BioMed Research International
 (2021)

---

## Abstract

Increasing evidence has shown that noncoding RNAs play significant roles in the initiation, progression, and metastasis of tumours via participating in competing endogenous RNA (ceRNA) networks. However, the survival-associated ceRNA in lung adenocarcinoma (LUAD) remains poorly understood. In this study, we aimed to investigate the regulatory mechanisms underlying ceRNA in LUAD to identify novel prognostic factors. mRNA, lncRNA, and miRNA sequencing data obtained from the GDC data portal were utilized to identify differentially expressed (DE) RNAs. Survival-related RNAs were recognized using univariate Kaplan-Meier survival analysis. We performed functional enrichment analysis of survival-related mRNAs using the clusterProfiler package of R and STRING. lncRNA-miRNA and miRNA-mRNA interactions were predicted based on miRcode, Starbase, and miRanda. Subsequently, the survival-associated ceRNA network was constructed for LUAD. Multivariate Cox regression analysis was used to identify prognostic factors. Finally, we acquired 15 DE miRNAs, 49 DE lncRNAs, and 843 DE mRNAs associated with significant overall survival. Functional enrichment analysis indicated that survival-related DE mRNAs were enriched in cell cycle. The survival-associated lncRNA-miRNA-mRNA ceRNA network was constructed using five miRNAs, 49 mRNAs, and 21 lncRNAs. Furthermore, seven hub RNAs (LINC01936, miR-20a-5p, miR-31-5p,TNS1, TGFBR2, SMAD7, and NEDD4L) were identified based on the ceRNA network. LINC01936 and miR-31-5p were found to be significant using the multifactorial Cox regression model. In conclusion, we successfully constructed a survival-related lncRNA-miRNA-mRNA ceRNA regulatory network in LUAD and identified seven hub RNAs, which provide novel insights into the regulatory molecular mechanisms associated with survival of LUAD, and identified two independent prognostic predictors for LUAD.


---

## Body

## 1. Introduction

Lung cancer is the most commonly diagnosed and lethal malignancy worldwide [1]. Lung adenocarcinoma (LUAD) is a common subtype of lung cancer [2]. Despite the recent advances in targeted therapeutic strategies, the outcomes of the available treatment strategies for LUAD remain unsatisfactory owing to the drug resistance and relapse, and the five-year overall survival is less than 20% [3]. Therefore, there is an urgent need to understand the molecular mechanisms underlying the pathogenesis of LUAD and identify novel potential prognostic biomarkers to improve prognosis of the disease.Genetic mutations and dysregulation that can contribute to the pathogenesis of cancer are served as biomarkers. Mutations in epidermal growth factor receptor (EGFR) occur in approximately 20% cases of lung cancer [4], and epidermal growth factor receptor-tyrosine kinase inhibitors (EGFR-TKIs) are indispensable in the treatment of EGFR-mutant advanced LUAD. Next-generation sequencing technology has been used to study the role of various RNAs in greater depth. Long noncoding RNAs (lncRNAs) have been considered as potential biomarkers and therapeutic targets due to their unique expression in various cells [5]. Several studies have indicated that dysregulation of lncRNAs, such as MIR31HG [6] and LINC01512 [7], promotes the progression and proliferation of tumour cells in LUAD. microRNAs (miRNAs) play a crucial role in the regulation of protein expression and therefore are considered potential biomarkers in cancer diagnosis. Wang et al. [8] identified a four-miRNA signature comprising miR-142-5p, miR-409-3p, miR-223-3p, and miR-146a-5p, for an early detection of LUAD. Xu et al. [9] reported that miRNA-21, miRNA-125b, and miRNA-224 are associated with chemotherapy sensitivity in patients with LUAD. However, the RNA biomarkers require a critical review before their application for clinical decision-making.The competing endogenous RNA (ceRNA) network hypothesis, which states that noncoding RNAs (ncRNAs), miRNAs, and mRNAs communicate with each other through microRNA response elements (MREs), has been implicated in posttranscriptional regulation [10, 11]. miRNAs repress the translation of target mRNAs by partial or complete complementary binding to MREs on their target RNA transcripts [12, 13]. ncRNAs can act as endogenous miRNA sponges to competitively bind miRNAs through shared MREs to regulate the expression levels of mRNAs, thereby forming specific ceRNA regulatory network comprising ncRNA-miRNA-mRNA interactions [14]. In the past decade, the study of ceRNAs has gained increased attention and several studies have reported their involvement in tumorigenesis [15], migration [16], and prognosis [17]. For example, lncRNA MAFG-AS1 regulates the expression of MAFG to facilitate proliferation of LUAD cells via miR-744-5p [18]. The ceRNA hypothesis provides novel insights into tumorigenesis [19] and biomarker identification [11] at the system biology level.Bioinformatics techniques are used to integrate and analyse large-scale genomic data, such as RNA-Seq and microarray, to discover potential molecular mechanisms and identify biomarkers, and to guide further experiments. Kumar et al. identified hub genes as potential biomarkers from a large number of differentially expressed (DE) genes by protein-protein interaction (PPI) network [20] and enrichment analysis [21, 22], providing valuable ideas for further study. Wan et al. [23] identified a prognosis-associated ceRNA axes in prostate cancer based on RNA sequencing data using bioinformatics approaches, and validated their regulatory mechanisms by cell proliferation and dual luciferase reporter assay.In this study, we aimed to construct a ceRNA network associated with survival in DE genes to reveal the molecular mechanisms underlying LUAD and initially identify prognostic factors, thereby to provide new ideas for further biological experiments. In addition to identifying the relationship among various RNAs based on the RNA interaction database, we also performed three statistical tests on lncRNA-mRNA pairs to screen for significant ceRNA interactions based on the ceRNA network hypothesis. The concise LUAD ceRNA network proposed by us would provide accurate and reliable results for subsequent studies.


## 2. Materials and Methods

### 2.1. Data Source and Preprocessing

The GDC data portal (https://portal.gdc.cancer.gov/) [24] is an accessible high-quality cancer genome data-sharing platform that provides primary processed genomic data (level 3 data). We acquired level 3 RNA-Seq (including mRNA and lncRNA) and miRNA-Seq RNA expression data (HTSeq-counts), and clinical information of LUAD patients, who were part of TCGA project from GDC on October 15, 2019. After excluding duplicate samples and other tissue samples, mRNA and lncRNA dataset included 524 cancer samples and 59 adjacent nontumour tissue samples and isoform quantification data from miRNA-Seq included 516 cancer samples and 46 adjacent nontumour tissue samples.The RNA-Seq and miRNA-Seq data were processed using the R package GDCRNATools [25]. The raw RNA counts were normalised using the trimmed mean of M value (TMM) method [26] and transformed via the voom method [27], wherein the RNAs with lower expression, where the log CPM was found to be lower than 1 in more than half of the samples, were filtered out. The procedure followed in this study is demonstrated in Figure 1.Figure 1
Flow chart of ceRNA network construction.


### 2.2. Differential Expression Analysis

The limma [28] method was used to identify DE RNAs. Fold change (FC) refers to the differences in RNA expression within samples, and the ∣FC∣>2 as the threshold value was set based on previous studies on the ceRNA network [29, 30]. ∣FC∣>2 and false discovery rate FDR<0.01 were considered statistically significant. Compared to adjacent nontumour tissue samples, RNAs with a higher expression level (FC>2) in tumour tissue samples were considered upregulated DE RNAs, whereas RNAs with lower expression level (FC<−2) were considered downregulated DE RNAs.


### 2.3. Survival Analysis of DE RNAs

Univariate Kaplan-Meier survival analysis was performed to determine the correlation between the expression level of each DE RNA and the survival time of patients with LUAD. LUAD patients were categorised into high- and low-expression groups based on the median expression of certain DE RNAs. The hazard ratio (HR) of the two groups was evaluated using the Kaplan-Meier plot, and their difference was assessed by performing the log-rank test using the survival package of R [31]. Results with p<0.05 were considered statistically significant.


### 2.4. Functional Enrichment Analysis

The clusterProfiler package [32] of R software is a widely used method for functional enrichment analysis. This package performs overrepresentation and hypergeometric tests to identify DE mRNAs enriched in biological functions or processes. Several enrichment methods ignore the numerical information of DE mRNAs. However, the STRING database (https://string-db.org/) [33] provides another platform to analyse the numerical data via the two-sided Kolmogorov-Smirnov test and aggregate fold change test that perform well in various settings [33, 34]. We used aforementioned two tools to perform functional enrichment analysis on survival-related DE mRNAs and their FC value (for STRING). The cut-off value was set as p adjusted < 0.01 for clusterProfiler and FDR<0.01 for STRING.


### 2.5. Construction of Survival-Associated ceRNA Networks and Identification of Prognostic Predictors

We used miRcode (http://www.mircode.org/) [35] to predict the potential interactions between survival-related miRNAs and lncRNAs. Starbase (http://starbase.sysu.edu.cn/) [36] and miRanda (http://www.microrna.org/microrna/home.do) [37] were used to predict target genes of survival-related miRNAs. Starbase uses multiple algorithms and Ago-binding sites to predict miRNA target sites and their target genes. miRanda predicts the miRNA-mRNA interactions with optimal sequence complementarity using a weighted dynamic programming algorithm and thermodynamic analysis. Therefore, we chose these two databases to improve the reliability of the prediction outcomes.The ceRNA hypothesis proposed that lncRNAs and their target mRNAs had a positive correlation and they shared miRNAs. Therefore, the competing endogenous interactions between lncRNA and mRNA were evaluated by performing three different statistical tests using GDCRNATools package to select ceRNA pairs matching the ceRNA hypothesis. First, a hypergeometric test was performed to test whether lncRNA and mRNA significantly share a number of miRNAs. Second, Pearson correlation analysis was performed to test the positive correlation between lncRNA and mRNA expressions. Third, regulation pattern analysis [38] was used to measure the regulatory role of miRNAs on lncRNAs and mRNAs. The test criteria were p<0.05, and regulation similarity was not equal to 0.The lncRNA-miRNA-mRNA ceRNA regulatory network associated with the survival of LUAD was constructed using Cytoscape 3.7.1 [39]. Hub nodes of the ceRNA network were identified using Cytoscape plugin cytoHubba [40]. We imported the mRNAs from the ceRNA network into the STRING database [33] and selected “Homo sapiens” in organism and medium confidence in the minimum required interaction score to obtain a PPI network. The hub RNAs were subjected to multivariate Cox regression analysis to identify independent prognostic predictors.


## 2.1. Data Source and Preprocessing

The GDC data portal (https://portal.gdc.cancer.gov/) [24] is an accessible high-quality cancer genome data-sharing platform that provides primary processed genomic data (level 3 data). We acquired level 3 RNA-Seq (including mRNA and lncRNA) and miRNA-Seq RNA expression data (HTSeq-counts), and clinical information of LUAD patients, who were part of TCGA project from GDC on October 15, 2019. After excluding duplicate samples and other tissue samples, mRNA and lncRNA dataset included 524 cancer samples and 59 adjacent nontumour tissue samples and isoform quantification data from miRNA-Seq included 516 cancer samples and 46 adjacent nontumour tissue samples.The RNA-Seq and miRNA-Seq data were processed using the R package GDCRNATools [25]. The raw RNA counts were normalised using the trimmed mean of M value (TMM) method [26] and transformed via the voom method [27], wherein the RNAs with lower expression, where the log CPM was found to be lower than 1 in more than half of the samples, were filtered out. The procedure followed in this study is demonstrated in Figure 1.Figure 1
Flow chart of ceRNA network construction.


## 2.2. Differential Expression Analysis

The limma [28] method was used to identify DE RNAs. Fold change (FC) refers to the differences in RNA expression within samples, and the ∣FC∣>2 as the threshold value was set based on previous studies on the ceRNA network [29, 30]. ∣FC∣>2 and false discovery rate FDR<0.01 were considered statistically significant. Compared to adjacent nontumour tissue samples, RNAs with a higher expression level (FC>2) in tumour tissue samples were considered upregulated DE RNAs, whereas RNAs with lower expression level (FC<−2) were considered downregulated DE RNAs.


## 2.3. Survival Analysis of DE RNAs

Univariate Kaplan-Meier survival analysis was performed to determine the correlation between the expression level of each DE RNA and the survival time of patients with LUAD. LUAD patients were categorised into high- and low-expression groups based on the median expression of certain DE RNAs. The hazard ratio (HR) of the two groups was evaluated using the Kaplan-Meier plot, and their difference was assessed by performing the log-rank test using the survival package of R [31]. Results with p<0.05 were considered statistically significant.


## 2.4. Functional Enrichment Analysis

The clusterProfiler package [32] of R software is a widely used method for functional enrichment analysis. This package performs overrepresentation and hypergeometric tests to identify DE mRNAs enriched in biological functions or processes. Several enrichment methods ignore the numerical information of DE mRNAs. However, the STRING database (https://string-db.org/) [33] provides another platform to analyse the numerical data via the two-sided Kolmogorov-Smirnov test and aggregate fold change test that perform well in various settings [33, 34]. We used aforementioned two tools to perform functional enrichment analysis on survival-related DE mRNAs and their FC value (for STRING). The cut-off value was set as p adjusted < 0.01 for clusterProfiler and FDR<0.01 for STRING.


## 2.5. Construction of Survival-Associated ceRNA Networks and Identification of Prognostic Predictors

We used miRcode (http://www.mircode.org/) [35] to predict the potential interactions between survival-related miRNAs and lncRNAs. Starbase (http://starbase.sysu.edu.cn/) [36] and miRanda (http://www.microrna.org/microrna/home.do) [37] were used to predict target genes of survival-related miRNAs. Starbase uses multiple algorithms and Ago-binding sites to predict miRNA target sites and their target genes. miRanda predicts the miRNA-mRNA interactions with optimal sequence complementarity using a weighted dynamic programming algorithm and thermodynamic analysis. Therefore, we chose these two databases to improve the reliability of the prediction outcomes.The ceRNA hypothesis proposed that lncRNAs and their target mRNAs had a positive correlation and they shared miRNAs. Therefore, the competing endogenous interactions between lncRNA and mRNA were evaluated by performing three different statistical tests using GDCRNATools package to select ceRNA pairs matching the ceRNA hypothesis. First, a hypergeometric test was performed to test whether lncRNA and mRNA significantly share a number of miRNAs. Second, Pearson correlation analysis was performed to test the positive correlation between lncRNA and mRNA expressions. Third, regulation pattern analysis [38] was used to measure the regulatory role of miRNAs on lncRNAs and mRNAs. The test criteria were p<0.05, and regulation similarity was not equal to 0.The lncRNA-miRNA-mRNA ceRNA regulatory network associated with the survival of LUAD was constructed using Cytoscape 3.7.1 [39]. Hub nodes of the ceRNA network were identified using Cytoscape plugin cytoHubba [40]. We imported the mRNAs from the ceRNA network into the STRING database [33] and selected “Homo sapiens” in organism and medium confidence in the minimum required interaction score to obtain a PPI network. The hub RNAs were subjected to multivariate Cox regression analysis to identify independent prognostic predictors.


## 3. Results and Discussion

### 3.1. Survival-Related DE RNAs

Differential expression analysis identified 1097 (37.15%) upregulated and 1856 (62.85%) downregulated DE mRNAs (Figure2(a)), 104 (55.91%) upregulated and 82 (44.09%) downregulated DE lncRNAs (Figure 2(b)), and 93 (62.00%) upregulated and 57 (38.00%) downregulated DE miRNAs (Figure 2(c)) between tissue samples and adjacent nontumour tissue samples. Heat maps of the three DE RNAs are shown in Figures S1–S3 in the Supplementary Materials. We further analysed the associations between these DE RNAs and survival time using univariate Kaplan-Meier survival analysis. In total, we identified 15 DE miRNAs, 49 DE lncRNAs, and 843 DE mRNAs with significant overall survival in 511 patients with LUAD for subsequent analysis.Figure 2
Volcano plots of the DE RNAs: (a) DE mRNAs, (b) DE lncRNAs, and (c) DE miRNAs. The green dots indicated downregulated DE RNAs and the red indicated upregulated DE RNAs in tumour samples. The black dots indicated excluded RNAs.
(a)(b)(c)


### 3.2. Functional Enrichment Analysis of Survival-Related mRNAs

The 843 mRNAs with significant overall survival were analysed using clusterProfiler and STRING for Gene Ontology (GO) enrichment analysis and identified the first 10 terms withp values among three different categories. GO included three different aspects: biological process (BP), cellular component (CC), and molecular function (MF). Figure 3 shows that 36 upregulated mRNAs with high log2FC values were enriched in six terms associated with BP and four cellular components as identified via clusterProfiler. The four terms associated with BP (sister chromatid segregation, nuclear chromosome segregation, mitotic nuclear division, and mitotic sister chromatid segregation) were involved in cell cycle. The four terms associated with CC were involved in chromosome. Figure 4 shows 60 mRNAs enriched in eight terms associated with BP and the two terms associated with CC as identified by STRING. Similar to the clusterProfiler results, three terms associated with BP involved in cell cycle (mitotic cell cycle, cell cycle process, and mitotic cell cycle process) and two terms associated with CC involved in chromosome. For the MF ontology, the catalytic activity acting on DNA was identified using two different methods.Figure 3
The first 10 GO terms of survival-related DE mRNAs in LUAD from clusterProfiler.Figure 4
The first 10 terms of survival-related DE mRNAs from STRING.The clusterProfiler Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway results (Figure5) indicated that mRNAs were mainly enriched in cell cycle, human T cell leukaemia virus 1 infection, and cell adhesion molecules (CAMs). STRING identified two different pathways, namely, cell cycle and oocyte meiosis.Figure 5
The KEGG pathways of survival-related DE mRNAs in LUAD. Thex-axis indicates the number of mRNAs participated in the pathway.The two methods revealed that DE mRNAs were enriched in the cell cycle process and downstream terms. Cell cycle involves progression of cell and nuclear replication and dysregulation of cell replication and division contributes to tumorigenesis. Several studies have reported that overexpression of cell division cycle-associated genes is associated with tumour cell proliferation indicating poor survival in lung cancer patients [41–43]. Rac3 induces apoptosis of LUAD cells via cell cycle pathway and is associated with longer survival [44]. Other KEGG pathways have rarely been mentioned in earlier studies on LUAD. A previous study indicated that human T cell leukaemia virus type I infection induces gene expression of CAMs in lung epithelial cells [45]. However, the association between these pathways and LUAD has not been reported previously.


### 3.3. Survival-Related lncRNA-miRNA-mRNA ceRNA Network in LUAD

To investigate the regulatory interaction in survival-related RNAs, we acquired information on lncRNA-miRNA interactions from miRcode and miRNA-mRNA pairs from Starbase and miRanda. Three statistical tests were performed on lncRNA-mRNA pairs to confirm significant ceRNA pairs. The aforementioned results intersected with survival-significant RNAs. Finally, we considered five miRNAs, 49 mRNAs, and 21 lncRNAs (Table1) to construct a survival-related ceRNA network (Figure 6) comprising 37 pairs of miRNA-lncRNA interaction and 61 pairs of miRNA-mRNA interaction. This network suggests a potential regulatory relationship between lncRNA-miRNA-mRNA in LUAD prognosis. Several RNAs in the ceRNA network have been verified for their regulatory role in lung cancer or other cancers. LINC00857 regulates cell growth, glycolysis, and apoptosis in LUAD [46]. lncRNA MNX1-AS1 regulates the progression of oesophageal squamous cell carcinoma by targeting the miR-34a/SIRT1 axis [47]. However, most of these interactions have not been previously reported to be associated with LUAD.Table 1
The characteristic of RNAs in the survival-related ceRNA network.
SymbolDifferential expression analysisKaplan-Meier survival analysisLog2FCFDRHRLower 95Upper 95p valuemiRNAmiR-31-5p3.03909.65E−121.41501.05811.89230.0200miR-21-5p2.56642.29E−471.43471.07131.92140.0133miR-148a-3p1.42225.86E−190.62380.46590.83520.0014miR-20a-5p1.12512.26E−111.47031.09651.97150.0084miR-133a-3p-2.55161.71E−280.73330.54840.98060.0380lncRNAMNX1-AS14.35902.83E−401.33671.00051.78580.0483VPS9D1-AS12.49684.97E−231.36031.01851.81690.0364LINC008572.03584.28E−251.41441.05941.88830.0196AP004608.11.70269.91E−040.70470.52780.94100.0176AC012236.11.41834.53E−100.64750.48470.86500.0032MELTF-AS11.27534.49E−131.34541.00731.79680.0438SNHG121.04741.98E−100.68470.51280.91430.0106CYP1B1-AS1-1.03878.08E−100.74640.55860.99750.0464HAGLR-1.11451.09E−050.73300.54900.97870.0358AC105942.1-1.30811.45E−280.70630.52880.94330.0183LINC01852-1.36488.65E−370.71090.53190.95000.0202AC021016.2-1.72205.30E−590.65660.49150.87720.0043MIR99AHG-1.79722.33E−260.57270.42830.76580.0002COLCA1-1.80734.92E−150.58090.43510.77570.0003AC090559.1-1.89114.83E−360.69590.52130.92920.0152AC093278.2-1.98081.96E−640.73360.54950.97940.0366AC125807.2-2.27456.44E−701.45211.08671.94030.0111LINC00261-2.60411.66E−180.69570.52080.92930.0138C8orf34-AS1-2.82818.13E−240.70720.52930.94500.0184LINC01936-2.83181.24E−610.69480.52040.92770.0141LHFPL3-AS2-4.08094.86E−530.65550.49080.87540.0044mRNAZFPM2-AS14.34128.62E−261.33741.00181.78550.0497COL1A12.71771.04E−391.45561.08931.94500.0104CCNA22.55892.64E−221.74941.30912.33790.0002E2F72.21241.04E−181.69521.26652.26900.0003RALGPS21.52461.15E−201.53041.14412.04700.0036PTGFRN1.03632.55E−181.41891.06231.89510.0173PSRC11.02557.85E−111.48731.11381.98600.0075PDE4B-1.01522.67E−140.65640.49160.87640.0046TMEM64-1.02003.00E−141.39771.04601.86780.0225SH2B3-1.07342.01E−270.73330.54930.97910.0371AKAP13-1.11361.52E−190.72330.54180.96570.0294SATB1-1.14136.49E−180.68960.51620.92120.0117MAP3K8-1.14186.74E−260.72730.54470.97110.0313ZC3H12C-1.14415.24E−221.50901.12942.01630.0051SMAD7-1.14656.02E−340.71370.53450.95290.0227PIK3R1-1.16019.04E−280.67370.50440.89960.0076ZNF704-1.16696.09E−130.62190.46450.83260.0011ACSL4-1.17653.70E−181.39091.04091.85860.0246HLA-DQB1-1.19233.72E−080.71610.53620.95640.0237DHCR24-1.21715.25E−220.73680.55110.98510.0370NEDD4L-1.22271.94E−200.67130.50260.89640.0070PRDM6-1.23773.06E−160.73960.55390.98760.0413LATS2-1.24164.69E−461.37261.02681.83480.0307VAMP2-1.26089.82E−360.74390.55710.99330.0445ELMO1-1.26973.39E−240.60720.45440.81150.0007CRY2-1.28111.89E−340.63540.47590.84850.0023KAT2B-1.30005.49E−380.63580.47600.84910.0022PREX1-1.33281.98E−270.74110.55500.98960.0446DNAJB4-1.40161.31E−321.59981.19812.13620.0016ST6GALNAC6-1.44231.71E−460.65090.48730.86930.0037SYNE1-1.44353.67E−170.67580.50620.90230.0088BTG2-1.48921.34E−180.64890.48580.86660.0034SH3BP5-1.58201.04E−400.73270.54870.97840.0351APBB1-1.58528.53E−390.70270.52620.93850.0169DOCK4-1.73937.68E−500.72970.54660.97430.0338DPYSL2-1.79512.05E−350.62230.46590.83110.0014FAM184A-1.80981.53E−250.65440.48980.87420.0040TGFBR2-1.85278.59E−440.71180.53310.95050.0214EGR2-1.90115.23E−380.74710.55960.99740.0490CADM1-1.97841.95E−210.69750.52230.93160.0145ADGRD1-2.11153.18E−210.58090.43480.77620.0002ERG-2.26313.57E−830.68200.51080.91060.0099TNS1-2.37442.53E−470.67080.50240.89570.0071WWC2-2.48012.28E−841.41891.06091.89760.0168KLF4-2.60711.81E−451.38001.03341.84280.0290SEMA6D-2.69469.92E−610.63910.47840.85360.0025ANKRD29-3.03746.21E−500.67550.50520.90310.0074LAMP3-3.44111.33E−420.74870.56070.99970.0496NCKAP5-4.47801.21E−1160.68210.51090.91060.0101FC: fold change; HR: hazard ratio; FDR: false discovery rate.Figure 6
Survival-related lncRNA-miRNA-mRNA ceRNA network of LUAD. The diamonds represent lncRNAs, ellipse mRNAs, and square miRNAs. The red represents upregulated genes, and the blue represents downregulated genes. The shade of the colour represents the magnitude of the value.The number of mRNAs in the ceRNA network was insignificant to perform functional enrichment analysis via clusterProfiler or STRING. Therefore, we used the PANTHER classification system [48] available on the GO website (http://geneontology.org/). The cellular components and KEGG pathways showed no significant results. All 48 genes, except ZFPM2-AS1, were identified and enriched in 38 terms associated with BP and three terms associated with MF. The results were sorted based on hierarchical relation of terms via PANTHER, with the parent term indented below the subclass (Table 2). Major genes were enriched in the regulation of cellular process (GO:0051244), response to stimulus (GO:0050896), signalling (GO:0023052), and their subclass. The regulation of cellular process involves the regulation of the rate, frequency, and extent of cellular processes. The signalling is a process that transmits information in biological systems. Moreover, the end of signal transduction (GO:0007165) regulates the initiation of transcription [49, 50]. In general, genes in the ceRNA network regulate the activity of various enzymes, participate in signal transduction, and indirectly regulate the initiation of transcription.Table 2
GO enrichment analysis of mRNAs in survival-related ceRNA network.
GO termsCountFDRBiological processRegulation of cyclin-dependent protein serine/threonine kinase activity44.42E−02Regulation of cyclin-dependent protein kinase activity44.57E−02Regulation of cellular process403.16E−02Biological regulation423.99E−02Regulation of protein kinase activity94.20E−02Regulation of kinase activity116.06E−03Regulation of transferase activity123.64E−03Regulation of catalytic activity193.22E−03Regulation of molecular function224.67E−03Regulation of protein metabolic process174.94E−02Regulation of developmental growth71.11E−02Regulation of developmental process173.08E−02Regulation of growth103.67E−03Response to peptide hormone84.51E−03Response to peptide92.89E−03Response to chemical264.58E−03Response to stimulus363.70E−03Response to organonitrogen compound111.12E−02Response to organic substance212.89E−03Response to nitrogen compound112.17E−02Response to hormone113.42E−03Response to endogenous stimulus143.83E−03Regulation of cell growth73.11E−02Negative regulation of catalytic activity94.37E−02Negative regulation of transcription by RNA polymerase II103.35E−02Negative regulation of nitrogen compound metabolic process164.36E−02Negative regulation of protein metabolic process112.51E−02Intracellular signal transduction162.70E−03Signal transduction282.62E−03Signalling292.70E−03Cell communication292.93E−03Cellular response to stimulus302.19E−02Negative regulation of multicellular organismal process114.31E−02Regulation of multicellular organismal process201.06E−02Regulation of cell differentiation151.14E−02Cell differentiation229.65E−03Cellular developmental process221.17E−02Positive regulation of biological process291.61E−02Molecular functionGuanyl-nucleotide exchange factor activity62.16E−02Enzyme binding214.54E−05Kinase binding102.18E−02FDR: false discovery rate.


### 3.4. Hub RNAs of ceRNA Network and Prognostic Predictor

A subnetwork with 15 hub nodes (Figure7, Table 3) was identified using maximal clique centrality (MCC) in cytoHubba plug-ins. There were a total of six ceRNA pairs in the subnetwork, among which LINC01936-TNS1 exhibited the highest correlation coefficient (Figure 8) indicating they might have the same expression patterns. LINC01936 and TNS1 were the highest scoring nodes in their respective categories. The ceRNA network suggested LINC01936 and TNS1 interacted with miR-20a-5p and miR-31-5p. Our study demonstrated that lower expression of LINC01936 was associated with longer overall survival (Figure 9(a)). However, the role of LINC01936 in LUAD remains unclear. miR-20a-5p exhibited highest topological parameters, indicating that it plays a crucial role in the ceRNA network. Overexpression of miR-20a-5p promotes the migration and invasion of tumour cells [51] and correlates with a shorter survival [52], which is consistent with the findings of our study (Figure 9(c)). The low expression group of miR-31-5p showed better survival (Figure 9(b)). Wei et al. reported that miR31-5p was upregulated in LUAD patients with lymph node metastasis, and low expression of miR-31 was associated with good prognosis in patients with T2N0 stage [53]. TNS1 participates in fibrillar adhesion formation and cell migration [54] and is involved in signal transduction [55]. However, the role of TNS1 in tumours remains controversial. TNS1 negatively regulates tumour migration and invasion, and its high expression is associated with longer metastasis-free survival in breast cancer [56]. Moreover, higher expression of TNS1 is associated with worse prognosis in colon adenocarcinoma [57]. In this study, we observed that TNS1 was downregulated in LUAD tissues and that higher expression was associated with better prognosis (Figure 9(d)). Based on our analysis, we predicted that LINC01936 regulates TNS1 via miR-20a-5p and miR31-5p. However, the role of TNS1 and its regulatory interaction in LUAD remains unclear and warrants further in vitro and in vivo studies.Figure 7
The subnetwork of survival-related ceRNA network. The shade of the colour represents the magnitude of MCC score.Table 3
10 hub RNAs in the ceRNA network ranked by the MCC method.
RankHub RNAsMCC scoreBetweenness centralityCloseness centralityDegree1miR-20a-5p340.57050.4933342miR-148a-3p250.40530.4405253miR-21-5p180.24920.4066184miR-31-5p140.18580.3895145miR-133a-3p70.07270.362776LINC0193630.07240.448536HAGLR30.05070.391536CYP1B1-AS130.04820.371936TNS130.05060.4134310PREX120.01910.3507210MIR99AHG20.0160.3231210COLCA120.01040.3348210AC021016.220.04140.4134210SNHG1220.00860.3203210AC125807.220.02060.36822MCC: maximal clique centrality.Figure 8
Correlation analysis between lncRNAs and mRNAs of subnetwork: (a) LINC01936 and TNS1, (b) LINC01936 and PREX1, (c) MIR99HG and TNS1, and (d) MIR99HG and PREX1.
(a)(b)(c)(d)Figure 9
Kaplan-Meier curves of 7 hub RNAs: (a) LINC01936, (b) miR-31-5p, (c) miR-20a-5p, (d) TNS1, (e) TGFBR2, (f) SMAD7, and (g) NEDD4L.
(a)(b)(c)(d)(e)(f)(g)The scale of the ceRNA network constructed in this study was small. Some information may have been lost during the identification of hub genes using network topological parameters alone. Therefore, we included genes from the ceRNA network into the STRING to analyse protein interactions to identify hub genes. The PPI network (Figure10) demonstrated the protein interactions of the ceRNA network. TGFBR2, SMAD7, and NEDD4L with the highest degree were identified as hub genes in the PPI network. TGFBR2 is the crucial receptor for transforming growth factor-β 1 (TGF-β1). The TGF-β1 ligand binding to TGFBR2 depends on the serine and threonine residues of the receptor, which in turn binds to the TGF-β receptor I to initiate downstream signalling such as Smad and non-Smad signalling pathways to regulate cell proliferation, migration, and apoptosis [58, 59]. TGFBR2 is downregulated in various cancers [60]. Borczuk et al. reported that low expression of TGFBR2 associated with lymph node metastasis in patients with LUAD and increased risk of death [61]. Smad complexes translocate to the nucleus to initiate gene transcription. SMAD7 is an inhibitory Smad molecule that inhibits the formation of Smad complex [62]. Inhibition of miR-21 leads to SMAD7 upregulation, which inhibits cell invasion via TGF-β receptor signalling in non-small-cell lung cancer [63]. A previous study demonstrated that NEDD4L can limit TGF-β signalling by activating SMAD2/3 [64]. Downregulated NEDD4L enhances tumour metastasis and results in poor prognosis [65].Figure 10
Protein-protein interaction network of survival-related mRNAs in ceRNA network. Size of node represents degree value of mRNA. Line thickness represents the strength of data support (combine score).The three hub genes were downregulated in LUAD, and their high expression levels indicated a longer survival (Figures9(e)–9(g)). Moreover, the ceRNA network showed that LINC01936 is a ceRNA of the three hub genes and mediated through its interaction with miR-20a-5p. In conclusion, several interactions regulate the three hub genes by competitive binding of LINC01936 to miR-20a-5p, which in turn regulate TGF-β signalling and downstream signalling pathways. This affects LUAD progression and patient prognosis. To the best of our knowledge, these interactions have not been reported earlier. Thus, our study outcomes lay a strong foundation for future research studies in this field. Moreover, further studies are required to confirm whether silencing or overexpression of LINC01936 affects the expression of hub genes.Multifactorial Cox regression analysis was performed to identify independent prognostic factors from the abovementioned seven hub RNAs. Our results (Table4), based on the multifactorial Cox regression model, indicate that LINC01936 and miR-31-5p are independent prognostic predictors of LUAD. LINC01936 was identified as a protective predictor for LUAD, while miR-31-5p was identified as a risk factor. The potential of miR-31-5p as a biomarker has been reported in oral carcinoma [66], colorectal cancer [67], and lung cancer [68]. This study is the first to demonstrate the prognostic potential of LINC01936 in LUAD.Table 4
Multifactorial Cox regression analysis in hub RNAs.
RNAβHR95% CI of HRp valueTNS1-0.024590.97570.8082-1.17790.79806LINC01936-0.204870.81480.7171-0.92570.00165TGFBR20.155141.16780.9465-1.44090.14790SMAD70.109761.1160.8812-1.41340.36253NEDD4L-0.097420.90720.7643-1.07680.26538miR-20a-5p0.127561.13610.9683-1.33290.11761miR-31-5p0.075021.07791.0136-1.14620.01675HR: hazard ratio; CI: confidence interval.Our study has several limitations. The three RNA expression data and clinical data for this study were based on TCGA database, and our findings lack biological validation. Computational prediction is only a preliminary step in ceRNA research. Therefore, these results need to be verified by studies involving large-scale clinical samples and laboratory methods such as qRT-PCR, luciferase reporter assay, and western blotting. The regulatory mechanism of the ceRNA network needs to be validated by furtherin vivo and in vitro research.


## 3.1. Survival-Related DE RNAs

Differential expression analysis identified 1097 (37.15%) upregulated and 1856 (62.85%) downregulated DE mRNAs (Figure2(a)), 104 (55.91%) upregulated and 82 (44.09%) downregulated DE lncRNAs (Figure 2(b)), and 93 (62.00%) upregulated and 57 (38.00%) downregulated DE miRNAs (Figure 2(c)) between tissue samples and adjacent nontumour tissue samples. Heat maps of the three DE RNAs are shown in Figures S1–S3 in the Supplementary Materials. We further analysed the associations between these DE RNAs and survival time using univariate Kaplan-Meier survival analysis. In total, we identified 15 DE miRNAs, 49 DE lncRNAs, and 843 DE mRNAs with significant overall survival in 511 patients with LUAD for subsequent analysis.Figure 2
Volcano plots of the DE RNAs: (a) DE mRNAs, (b) DE lncRNAs, and (c) DE miRNAs. The green dots indicated downregulated DE RNAs and the red indicated upregulated DE RNAs in tumour samples. The black dots indicated excluded RNAs.
(a)(b)(c)


## 3.2. Functional Enrichment Analysis of Survival-Related mRNAs

The 843 mRNAs with significant overall survival were analysed using clusterProfiler and STRING for Gene Ontology (GO) enrichment analysis and identified the first 10 terms withp values among three different categories. GO included three different aspects: biological process (BP), cellular component (CC), and molecular function (MF). Figure 3 shows that 36 upregulated mRNAs with high log2FC values were enriched in six terms associated with BP and four cellular components as identified via clusterProfiler. The four terms associated with BP (sister chromatid segregation, nuclear chromosome segregation, mitotic nuclear division, and mitotic sister chromatid segregation) were involved in cell cycle. The four terms associated with CC were involved in chromosome. Figure 4 shows 60 mRNAs enriched in eight terms associated with BP and the two terms associated with CC as identified by STRING. Similar to the clusterProfiler results, three terms associated with BP involved in cell cycle (mitotic cell cycle, cell cycle process, and mitotic cell cycle process) and two terms associated with CC involved in chromosome. For the MF ontology, the catalytic activity acting on DNA was identified using two different methods.Figure 3
The first 10 GO terms of survival-related DE mRNAs in LUAD from clusterProfiler.Figure 4
The first 10 terms of survival-related DE mRNAs from STRING.The clusterProfiler Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway results (Figure5) indicated that mRNAs were mainly enriched in cell cycle, human T cell leukaemia virus 1 infection, and cell adhesion molecules (CAMs). STRING identified two different pathways, namely, cell cycle and oocyte meiosis.Figure 5
The KEGG pathways of survival-related DE mRNAs in LUAD. Thex-axis indicates the number of mRNAs participated in the pathway.The two methods revealed that DE mRNAs were enriched in the cell cycle process and downstream terms. Cell cycle involves progression of cell and nuclear replication and dysregulation of cell replication and division contributes to tumorigenesis. Several studies have reported that overexpression of cell division cycle-associated genes is associated with tumour cell proliferation indicating poor survival in lung cancer patients [41–43]. Rac3 induces apoptosis of LUAD cells via cell cycle pathway and is associated with longer survival [44]. Other KEGG pathways have rarely been mentioned in earlier studies on LUAD. A previous study indicated that human T cell leukaemia virus type I infection induces gene expression of CAMs in lung epithelial cells [45]. However, the association between these pathways and LUAD has not been reported previously.


## 3.3. Survival-Related lncRNA-miRNA-mRNA ceRNA Network in LUAD

To investigate the regulatory interaction in survival-related RNAs, we acquired information on lncRNA-miRNA interactions from miRcode and miRNA-mRNA pairs from Starbase and miRanda. Three statistical tests were performed on lncRNA-mRNA pairs to confirm significant ceRNA pairs. The aforementioned results intersected with survival-significant RNAs. Finally, we considered five miRNAs, 49 mRNAs, and 21 lncRNAs (Table1) to construct a survival-related ceRNA network (Figure 6) comprising 37 pairs of miRNA-lncRNA interaction and 61 pairs of miRNA-mRNA interaction. This network suggests a potential regulatory relationship between lncRNA-miRNA-mRNA in LUAD prognosis. Several RNAs in the ceRNA network have been verified for their regulatory role in lung cancer or other cancers. LINC00857 regulates cell growth, glycolysis, and apoptosis in LUAD [46]. lncRNA MNX1-AS1 regulates the progression of oesophageal squamous cell carcinoma by targeting the miR-34a/SIRT1 axis [47]. However, most of these interactions have not been previously reported to be associated with LUAD.Table 1
The characteristic of RNAs in the survival-related ceRNA network.
SymbolDifferential expression analysisKaplan-Meier survival analysisLog2FCFDRHRLower 95Upper 95p valuemiRNAmiR-31-5p3.03909.65E−121.41501.05811.89230.0200miR-21-5p2.56642.29E−471.43471.07131.92140.0133miR-148a-3p1.42225.86E−190.62380.46590.83520.0014miR-20a-5p1.12512.26E−111.47031.09651.97150.0084miR-133a-3p-2.55161.71E−280.73330.54840.98060.0380lncRNAMNX1-AS14.35902.83E−401.33671.00051.78580.0483VPS9D1-AS12.49684.97E−231.36031.01851.81690.0364LINC008572.03584.28E−251.41441.05941.88830.0196AP004608.11.70269.91E−040.70470.52780.94100.0176AC012236.11.41834.53E−100.64750.48470.86500.0032MELTF-AS11.27534.49E−131.34541.00731.79680.0438SNHG121.04741.98E−100.68470.51280.91430.0106CYP1B1-AS1-1.03878.08E−100.74640.55860.99750.0464HAGLR-1.11451.09E−050.73300.54900.97870.0358AC105942.1-1.30811.45E−280.70630.52880.94330.0183LINC01852-1.36488.65E−370.71090.53190.95000.0202AC021016.2-1.72205.30E−590.65660.49150.87720.0043MIR99AHG-1.79722.33E−260.57270.42830.76580.0002COLCA1-1.80734.92E−150.58090.43510.77570.0003AC090559.1-1.89114.83E−360.69590.52130.92920.0152AC093278.2-1.98081.96E−640.73360.54950.97940.0366AC125807.2-2.27456.44E−701.45211.08671.94030.0111LINC00261-2.60411.66E−180.69570.52080.92930.0138C8orf34-AS1-2.82818.13E−240.70720.52930.94500.0184LINC01936-2.83181.24E−610.69480.52040.92770.0141LHFPL3-AS2-4.08094.86E−530.65550.49080.87540.0044mRNAZFPM2-AS14.34128.62E−261.33741.00181.78550.0497COL1A12.71771.04E−391.45561.08931.94500.0104CCNA22.55892.64E−221.74941.30912.33790.0002E2F72.21241.04E−181.69521.26652.26900.0003RALGPS21.52461.15E−201.53041.14412.04700.0036PTGFRN1.03632.55E−181.41891.06231.89510.0173PSRC11.02557.85E−111.48731.11381.98600.0075PDE4B-1.01522.67E−140.65640.49160.87640.0046TMEM64-1.02003.00E−141.39771.04601.86780.0225SH2B3-1.07342.01E−270.73330.54930.97910.0371AKAP13-1.11361.52E−190.72330.54180.96570.0294SATB1-1.14136.49E−180.68960.51620.92120.0117MAP3K8-1.14186.74E−260.72730.54470.97110.0313ZC3H12C-1.14415.24E−221.50901.12942.01630.0051SMAD7-1.14656.02E−340.71370.53450.95290.0227PIK3R1-1.16019.04E−280.67370.50440.89960.0076ZNF704-1.16696.09E−130.62190.46450.83260.0011ACSL4-1.17653.70E−181.39091.04091.85860.0246HLA-DQB1-1.19233.72E−080.71610.53620.95640.0237DHCR24-1.21715.25E−220.73680.55110.98510.0370NEDD4L-1.22271.94E−200.67130.50260.89640.0070PRDM6-1.23773.06E−160.73960.55390.98760.0413LATS2-1.24164.69E−461.37261.02681.83480.0307VAMP2-1.26089.82E−360.74390.55710.99330.0445ELMO1-1.26973.39E−240.60720.45440.81150.0007CRY2-1.28111.89E−340.63540.47590.84850.0023KAT2B-1.30005.49E−380.63580.47600.84910.0022PREX1-1.33281.98E−270.74110.55500.98960.0446DNAJB4-1.40161.31E−321.59981.19812.13620.0016ST6GALNAC6-1.44231.71E−460.65090.48730.86930.0037SYNE1-1.44353.67E−170.67580.50620.90230.0088BTG2-1.48921.34E−180.64890.48580.86660.0034SH3BP5-1.58201.04E−400.73270.54870.97840.0351APBB1-1.58528.53E−390.70270.52620.93850.0169DOCK4-1.73937.68E−500.72970.54660.97430.0338DPYSL2-1.79512.05E−350.62230.46590.83110.0014FAM184A-1.80981.53E−250.65440.48980.87420.0040TGFBR2-1.85278.59E−440.71180.53310.95050.0214EGR2-1.90115.23E−380.74710.55960.99740.0490CADM1-1.97841.95E−210.69750.52230.93160.0145ADGRD1-2.11153.18E−210.58090.43480.77620.0002ERG-2.26313.57E−830.68200.51080.91060.0099TNS1-2.37442.53E−470.67080.50240.89570.0071WWC2-2.48012.28E−841.41891.06091.89760.0168KLF4-2.60711.81E−451.38001.03341.84280.0290SEMA6D-2.69469.92E−610.63910.47840.85360.0025ANKRD29-3.03746.21E−500.67550.50520.90310.0074LAMP3-3.44111.33E−420.74870.56070.99970.0496NCKAP5-4.47801.21E−1160.68210.51090.91060.0101FC: fold change; HR: hazard ratio; FDR: false discovery rate.Figure 6
Survival-related lncRNA-miRNA-mRNA ceRNA network of LUAD. The diamonds represent lncRNAs, ellipse mRNAs, and square miRNAs. The red represents upregulated genes, and the blue represents downregulated genes. The shade of the colour represents the magnitude of the value.The number of mRNAs in the ceRNA network was insignificant to perform functional enrichment analysis via clusterProfiler or STRING. Therefore, we used the PANTHER classification system [48] available on the GO website (http://geneontology.org/). The cellular components and KEGG pathways showed no significant results. All 48 genes, except ZFPM2-AS1, were identified and enriched in 38 terms associated with BP and three terms associated with MF. The results were sorted based on hierarchical relation of terms via PANTHER, with the parent term indented below the subclass (Table 2). Major genes were enriched in the regulation of cellular process (GO:0051244), response to stimulus (GO:0050896), signalling (GO:0023052), and their subclass. The regulation of cellular process involves the regulation of the rate, frequency, and extent of cellular processes. The signalling is a process that transmits information in biological systems. Moreover, the end of signal transduction (GO:0007165) regulates the initiation of transcription [49, 50]. In general, genes in the ceRNA network regulate the activity of various enzymes, participate in signal transduction, and indirectly regulate the initiation of transcription.Table 2
GO enrichment analysis of mRNAs in survival-related ceRNA network.
GO termsCountFDRBiological processRegulation of cyclin-dependent protein serine/threonine kinase activity44.42E−02Regulation of cyclin-dependent protein kinase activity44.57E−02Regulation of cellular process403.16E−02Biological regulation423.99E−02Regulation of protein kinase activity94.20E−02Regulation of kinase activity116.06E−03Regulation of transferase activity123.64E−03Regulation of catalytic activity193.22E−03Regulation of molecular function224.67E−03Regulation of protein metabolic process174.94E−02Regulation of developmental growth71.11E−02Regulation of developmental process173.08E−02Regulation of growth103.67E−03Response to peptide hormone84.51E−03Response to peptide92.89E−03Response to chemical264.58E−03Response to stimulus363.70E−03Response to organonitrogen compound111.12E−02Response to organic substance212.89E−03Response to nitrogen compound112.17E−02Response to hormone113.42E−03Response to endogenous stimulus143.83E−03Regulation of cell growth73.11E−02Negative regulation of catalytic activity94.37E−02Negative regulation of transcription by RNA polymerase II103.35E−02Negative regulation of nitrogen compound metabolic process164.36E−02Negative regulation of protein metabolic process112.51E−02Intracellular signal transduction162.70E−03Signal transduction282.62E−03Signalling292.70E−03Cell communication292.93E−03Cellular response to stimulus302.19E−02Negative regulation of multicellular organismal process114.31E−02Regulation of multicellular organismal process201.06E−02Regulation of cell differentiation151.14E−02Cell differentiation229.65E−03Cellular developmental process221.17E−02Positive regulation of biological process291.61E−02Molecular functionGuanyl-nucleotide exchange factor activity62.16E−02Enzyme binding214.54E−05Kinase binding102.18E−02FDR: false discovery rate.


## 3.4. Hub RNAs of ceRNA Network and Prognostic Predictor

A subnetwork with 15 hub nodes (Figure7, Table 3) was identified using maximal clique centrality (MCC) in cytoHubba plug-ins. There were a total of six ceRNA pairs in the subnetwork, among which LINC01936-TNS1 exhibited the highest correlation coefficient (Figure 8) indicating they might have the same expression patterns. LINC01936 and TNS1 were the highest scoring nodes in their respective categories. The ceRNA network suggested LINC01936 and TNS1 interacted with miR-20a-5p and miR-31-5p. Our study demonstrated that lower expression of LINC01936 was associated with longer overall survival (Figure 9(a)). However, the role of LINC01936 in LUAD remains unclear. miR-20a-5p exhibited highest topological parameters, indicating that it plays a crucial role in the ceRNA network. Overexpression of miR-20a-5p promotes the migration and invasion of tumour cells [51] and correlates with a shorter survival [52], which is consistent with the findings of our study (Figure 9(c)). The low expression group of miR-31-5p showed better survival (Figure 9(b)). Wei et al. reported that miR31-5p was upregulated in LUAD patients with lymph node metastasis, and low expression of miR-31 was associated with good prognosis in patients with T2N0 stage [53]. TNS1 participates in fibrillar adhesion formation and cell migration [54] and is involved in signal transduction [55]. However, the role of TNS1 in tumours remains controversial. TNS1 negatively regulates tumour migration and invasion, and its high expression is associated with longer metastasis-free survival in breast cancer [56]. Moreover, higher expression of TNS1 is associated with worse prognosis in colon adenocarcinoma [57]. In this study, we observed that TNS1 was downregulated in LUAD tissues and that higher expression was associated with better prognosis (Figure 9(d)). Based on our analysis, we predicted that LINC01936 regulates TNS1 via miR-20a-5p and miR31-5p. However, the role of TNS1 and its regulatory interaction in LUAD remains unclear and warrants further in vitro and in vivo studies.Figure 7
The subnetwork of survival-related ceRNA network. The shade of the colour represents the magnitude of MCC score.Table 3
10 hub RNAs in the ceRNA network ranked by the MCC method.
RankHub RNAsMCC scoreBetweenness centralityCloseness centralityDegree1miR-20a-5p340.57050.4933342miR-148a-3p250.40530.4405253miR-21-5p180.24920.4066184miR-31-5p140.18580.3895145miR-133a-3p70.07270.362776LINC0193630.07240.448536HAGLR30.05070.391536CYP1B1-AS130.04820.371936TNS130.05060.4134310PREX120.01910.3507210MIR99AHG20.0160.3231210COLCA120.01040.3348210AC021016.220.04140.4134210SNHG1220.00860.3203210AC125807.220.02060.36822MCC: maximal clique centrality.Figure 8
Correlation analysis between lncRNAs and mRNAs of subnetwork: (a) LINC01936 and TNS1, (b) LINC01936 and PREX1, (c) MIR99HG and TNS1, and (d) MIR99HG and PREX1.
(a)(b)(c)(d)Figure 9
Kaplan-Meier curves of 7 hub RNAs: (a) LINC01936, (b) miR-31-5p, (c) miR-20a-5p, (d) TNS1, (e) TGFBR2, (f) SMAD7, and (g) NEDD4L.
(a)(b)(c)(d)(e)(f)(g)The scale of the ceRNA network constructed in this study was small. Some information may have been lost during the identification of hub genes using network topological parameters alone. Therefore, we included genes from the ceRNA network into the STRING to analyse protein interactions to identify hub genes. The PPI network (Figure10) demonstrated the protein interactions of the ceRNA network. TGFBR2, SMAD7, and NEDD4L with the highest degree were identified as hub genes in the PPI network. TGFBR2 is the crucial receptor for transforming growth factor-β 1 (TGF-β1). The TGF-β1 ligand binding to TGFBR2 depends on the serine and threonine residues of the receptor, which in turn binds to the TGF-β receptor I to initiate downstream signalling such as Smad and non-Smad signalling pathways to regulate cell proliferation, migration, and apoptosis [58, 59]. TGFBR2 is downregulated in various cancers [60]. Borczuk et al. reported that low expression of TGFBR2 associated with lymph node metastasis in patients with LUAD and increased risk of death [61]. Smad complexes translocate to the nucleus to initiate gene transcription. SMAD7 is an inhibitory Smad molecule that inhibits the formation of Smad complex [62]. Inhibition of miR-21 leads to SMAD7 upregulation, which inhibits cell invasion via TGF-β receptor signalling in non-small-cell lung cancer [63]. A previous study demonstrated that NEDD4L can limit TGF-β signalling by activating SMAD2/3 [64]. Downregulated NEDD4L enhances tumour metastasis and results in poor prognosis [65].Figure 10
Protein-protein interaction network of survival-related mRNAs in ceRNA network. Size of node represents degree value of mRNA. Line thickness represents the strength of data support (combine score).The three hub genes were downregulated in LUAD, and their high expression levels indicated a longer survival (Figures9(e)–9(g)). Moreover, the ceRNA network showed that LINC01936 is a ceRNA of the three hub genes and mediated through its interaction with miR-20a-5p. In conclusion, several interactions regulate the three hub genes by competitive binding of LINC01936 to miR-20a-5p, which in turn regulate TGF-β signalling and downstream signalling pathways. This affects LUAD progression and patient prognosis. To the best of our knowledge, these interactions have not been reported earlier. Thus, our study outcomes lay a strong foundation for future research studies in this field. Moreover, further studies are required to confirm whether silencing or overexpression of LINC01936 affects the expression of hub genes.Multifactorial Cox regression analysis was performed to identify independent prognostic factors from the abovementioned seven hub RNAs. Our results (Table4), based on the multifactorial Cox regression model, indicate that LINC01936 and miR-31-5p are independent prognostic predictors of LUAD. LINC01936 was identified as a protective predictor for LUAD, while miR-31-5p was identified as a risk factor. The potential of miR-31-5p as a biomarker has been reported in oral carcinoma [66], colorectal cancer [67], and lung cancer [68]. This study is the first to demonstrate the prognostic potential of LINC01936 in LUAD.Table 4
Multifactorial Cox regression analysis in hub RNAs.
RNAβHR95% CI of HRp valueTNS1-0.024590.97570.8082-1.17790.79806LINC01936-0.204870.81480.7171-0.92570.00165TGFBR20.155141.16780.9465-1.44090.14790SMAD70.109761.1160.8812-1.41340.36253NEDD4L-0.097420.90720.7643-1.07680.26538miR-20a-5p0.127561.13610.9683-1.33290.11761miR-31-5p0.075021.07791.0136-1.14620.01675HR: hazard ratio; CI: confidence interval.Our study has several limitations. The three RNA expression data and clinical data for this study were based on TCGA database, and our findings lack biological validation. Computational prediction is only a preliminary step in ceRNA research. Therefore, these results need to be verified by studies involving large-scale clinical samples and laboratory methods such as qRT-PCR, luciferase reporter assay, and western blotting. The regulatory mechanism of the ceRNA network needs to be validated by furtherin vivo and in vitro research.


## 4. Conclusions

In summary, our study constructed a survival-associated lncRNA-miRNA-mRNA ceRNA network in LUAD using bioinformatics approaches and identified seven hub RNAs (LINC01936, miR-20a-5p, miR-31-5p,TNS1, TGFBR2, SMAD7, and NEDD4L). LINC01936 and miR-31-5p were identified as independent prognostic predictors of LUAD. The ceRNA network identified in this study provides novel insights into the molecular regulatory mechanisms associated with LUAD progression. Further studies are required to explore the biological mechanisms of ceRNAs in LUAD and validate the prognostic value of LINC01936 and miR-31-5p in other cohorts.



---

*Source: 4093426-2021-02-12.xml*