Abstract
Background
Competing endogenous RNAs (ceRNAs) are a newly identified type of
regulatory RNA. Accumulating evidence suggests that ceRNAs play an
important role in the pathogenesis of diseases such as cancer. Thus,
ceRNA dysregulation may represent an important molecular mechanism
underlying cancer progression and poor prognosis. In this study, we
aimed to identify ceRNAs that may serve as potential biomarkers for
early diagnosis of lung adenocarcinoma (LUAD).
Methods
We performed differential gene expression analysis on TCGA-LUAD
datasets to identify differentially expressed (DE) mRNAs, lncRNAs, and
miRNAs at different tumor stages. Based on the ceRNA hypothesis and
considering the synergistic or feedback regulation of ceRNAs, a
lncRNA–miRNA–mRNA network was constructed. Functional analysis was
performed using gene ontology term and KEGG pathway enrichment analysis
and KOBAS 2.0 software. Transcription factor (TF) analysis was carried
out to identify direct targets of the TFs associated with LUAD
prognosis. Identified DE genes were validated using gene expression
omnibus (GEO) datasets.
Results
Based on analysis of TCGA-LUAD datasets, we obtained 2,610 DE mRNAs,
915 lncRNAs, and 125 miRNAs that were common to different tumor stages
(|log[2](Fold change)| ≥ 1, false discovery rate < 0.01), respectively.
Functional analysis showed that the aberrantly expressed mRNAs were
closely related to tumor development. Survival analyses of the
constructed ceRNA network modules demonstrated that five of them
exhibit prognostic significance. The five ceRNA interaction modules
contained one lncRNA (FENDRR), three mRNAs (EPAS1, FOXF1, and EDNRB),
and four miRNAs (hsa-miR-148a, hsa-miR-195, hsa-miR-196b, and
hsa-miR-301b). The aberrant expression of one lncRNA and three mRNAs
was verified in the LUAD GEO dataset. Transcription factor analysis
demonstrated that EPAS1 directly targeted 13 DE mRNAs.
Conclusion
Our observations indicate that lncRNA-related ceRNAs and TFs play an
important role in LUAD. The present study provides novel insights into
the molecular mechanisms underlying LUAD pathogenesis. Furthermore, our
study facilitates the identification of potential biomarkers for the
early diagnosis and prognosis of LUAD and therapeutic targets for its
treatment.
Keywords: ceRNA, Biomarkers, Molecular mechanisms, LUAD
Introduction
Lung cancer, which is one of the most common malignancies worldwide,
poses a serious threat to human health and life ([30]Jemal et al.,
2011; [31]Torre et al., 2015). In recent years, there has been an
increase in the morbidity of lung cancer owing to deterioration of the
environment ([32]Cui et al., 2017; [33]David et al., 2016; [34]Haugen
et al., 2000; [35]Jin et al., 2014; [36]Kooperstein, Schifrin & Leahy,
1965; [37]Weiss, 1983). Furthermore, lung cancer is difficult to
diagnose early and associated with poor prognosis and low survival rate
([38]Bernaudin, 2010; [39]Hardavella, George & Sethi, 2016;
[40]Onganer, Seckl & Djamgoz, 2005; [41]Prim et al., 2010; [42]Reinmuth
et al., 2014; [43]Tamura et al., 2015). The prognosis of patients with
lung cancer is not only related to its biological factors and immune
status, but is also closely related with the therapeutic factors
([44]Fan, Zhang & Liu, 2014; [45]Prim et al., 2010; [46]Ringbaek et
al., 1999; [47]Truong, Viswanathan & Erasmus, 2011); association
between multiple factors leads to complex changes in the disease and
makes prognosis difficult. Tumor staging is an important factor in
determining prognosis, including survival in cancer patients ([48]Gong
et al., 2017; [49]Li et al., 2017). Studies have shown that patients
with late stage lung cancer display extremely poor prognosis, with
invasion and metastasis, and additionally respond poorly to clinical
treatment ([50]Mattern, Koomagi & Volm, 2002; [51]Prim et al., 2010).
Therefore, early detection of cancer should enable administration of
effective treatment and considerably improve patient prognosis.
Although the prognosis of non-small cell lung cancer (NSCLC) has been
greatly improved following the introduction of multidisciplinary
treatment, the prognosis of lung cancer patients and treatment
efficiency remain poor ([52]Fan, Zhang & Liu, 2014; [53]Huang et al.,
2012; [54]Miura et al., 1996; [55]Pujol et al., 1994); therefore, the
identification of biological markers and therapeutic targets for the
early detection and treatment of recurrent lung tumor is critical.
Competing endogenous RNAs (ceRNAs) are a group of regulatory RNA
molecules that compete with other RNA molecules to bind specific
miRNAs, thereby regulating target gene expression ([56]Phelps et al.,
2016; [57]Shukla, Singh & Barik, 2011; [58]Tay et al., 2011). Many
studies have suggested that miRNA-mediated ceRNA regulatory mechanisms
play crucial roles in tumor occurrence and development ([59]Bartel,
2009; [60]Hansen et al., 2013; [61]Tay, Rinn & Pandolfi, 2014).
Therefore, the identification of ceRNA biomarkers related to lung
cancer should enable the development of effective strategies for the
diagnosis and treatment of lung cancer.
In this study, we aimed to identify differentially expressed genes
(DEGs) associated with lung adenocarcinoma (LUAD) development and
construct a ceRNA regulatory network using the TCGA-LUAD dataset.
Further, we examined the relationship of the identified genes and ceRNA
interaction modules with overall survival and prognosis. Through our
study, we sought to gain new insights into the molecular mechanisms
underlying LUAD and identify potential biomarkers for early diagnosis
of this disease.
Materials and Methods
Data source and preprocessing
Lung adenocarcinoma RNA-seq and lncRNA gene expression level three data
were obtained from the database Cancer Genome Atlas (TCGA,
[62]https://tcga-data.nci.nih.gov/docs/publications/tcga/). After an
initial filter (excluding samples with uncertain tumor stage), 569
samples were included for further study. Patient clinical data,
including outcome and staging information, were also downloaded. Of the
569 samples, 510 were LUAD samples (168 with T1 stage, 276 with T2
stage, 47 with T3 stage, and 19 with T4 stage) and 59 were normal
samples.
Screening of genetically altered genes
Gene raw read counts were used to perform differential expression
analysis with DESeq2 (v 1.18.1) ([63]Love, Huber & Anders, 2014), a R
package that uses a model based on the negative binomial distribution
and is widely used for RNA-seq data differential analysis.
Differentially expressed genes were defined by two criteria: false
discovery rate (FDR) < 0.01 and fold change (FC) ≥ 2. In order to
filter out the low-quality DEGs, we performed the following additional
filter strategy. For mRNA and miRNA, we only retained DEGs that showed
FPKM or RPKM larger than the threshold of one for at least 10% samples;
for lncRNA DEGs, the threshold was set to 0.1. We then compared the
normal tissue sample with different tumor stage samples (T1, T2, T3,
and T4) to identify common DEGs for further analysis.
Construction of lncRNA-associated ceRNA network (gene-set construction)
Based on the relationship between lncRNA, mRNA, miRNA, a ceRNA network
was constructed using the following steps: 1) Both up-regulated and
down-regulated LUAD-specific RNAs were chosen to construct the ceRNA
network; 2) miRanda software was used to predict the miRNA–lncRNA
interactions, and interactions in starBase ([64]Li et al., 2014) were
also included; 3) miRNA–mRNA interactions in miRTarbase and starBase
were used; 4) Correlation between LUAD-specific lncRNA and mRNA was
calculated with WGCNA ([65]Zhang & Horvath, 2005) corAndPvalue function
and correlation pairs with coefficient > 0.8 and p-value < 0.05 were
retained; 5) The hypergeometric method was used to test the enrichment
significance of miRNAs between mRNA and lncRNA. The enrichment
significance was calculated by using the following formula:
[MATH: P=∑i=c<
/mi>min(K
,n)(Ki)(N−K<
/mtd>n−i<
/mtr>)(Nn) :MATH]
In the above formula, N represents the total number of miRNAs, K is the
number of miRNAs targeting mRNA, n is the number of miRNAs targeting
lncRNA, and c is the number of miRNAs shared by the mRNA and lncRNA.
Only ceRNA modules with FDR < 0.05 were considered.
Gene functional enrichment analysis
In order to investigate key mRNAs at molecular and functional level,
gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)
pathway function enrichment analysis was performed. GOseq, which uses
Wallenius’ noncentral hypergeometric distribution model taking gene
length bias into account, was used to perform GO enrichment analysis.
Gene ontology terms with p-value less than 0.05 were considered
significantly enriched. KOBAS 2.0 software was used to test statistical
enrichment in KEGG pathways, and those with a Fisher’s exact test
p-value of less than 0.05 was considered to be significantly enriched
([66]Xie et al., 2011).
Transcription factor (TF) analysis
Co-expression networks of LUAD-mRNA expression data sets were generated
using Genie3 with EPAS1 and FOXF1 as Transcription factors (TFs) of
interest. A threshold of 0.005 was used to filter the co-expression
network. Additionally, Spearman correlations between the TFs and
LUAD-specific mRNAs were calculated, and TF-mRNA pairs with coefficient
> 0.03 were retained as this value represents positive correlation.
Using this procedure, we obtained a TF co-expression draft regulatory
network. As the analysis was only based on co-expression, the results
may include numerous indirect targets (such as downstream effects). To
identify genes that are most likely the direct targets, we used
RcisTarget ([67]Aibar et al., 2017; [68]Imrichova et al., 2015) to
perform cis-regulatory motif analysis of each TF regulon. RcisTarget is
an R package that identifies TF binding motifs that are
over-represented in a gene list. We used the database, scoring 500 bp
upstream of TSS, and only motifs with a Normalized Enrichment Score
(NES) > 3.0 were considered associated with TFs of interest.
Validation of gene expression with gene expression omnibus data
In order to validate the key differentially expressed (DE) molecules,
we examined four datasets from among the gene expression omnibus (GEO)
datasets. Of the selected datasets, the [69]GSE10072 study contained 58
LUAD tumor and 49 non-tumor tissue mRNA expression data, [70]GSE32863
contained 58 LUAD tumor and 58 non-tumor tissue mRNA expression data,
[71]GSE85716 contained six LUAD tumor and six non-tumor tissue lncRNA
expression data, and [72]GSE104854 contained three LUAD (pooled from
nine tissues) and three non-tumor (pooled from nine tissues) expression
data. [73]GSE10072, [74]GSE32863, and [75]GSE85716 datasets, which were
derived from microarray data, were analyzed with GEO2R online. The
[76]GSE104854 dataset, which was obtained by high-throughput
sequencing, was analyzed with the DESeq2 R package based on read
counts.
Statistical analysis
To investigate the impact of the expression level of RNAs on the
prognostic survival of patients, survival analysis was performed. In
the analysis, Kaplan–Meier survival analyses and log-rank test were
performed to study the relationship between RNA expression states
(cutoff point: median value) and survival time. Univariate cox
proportional hazards regression was applied to identify the RNAs
associated with overall survival. In order to evaluate the potential of
dysregulated ceRNA modules as biomarkers, the risk scoring classifier
was constructed. For each ceRNA module, we calculated the risk score
for multivariate survival analysis to determine the prognostic
influence of the RNAs as a whole. All statistical analyses were
performed using R version 3.4.3 software ([77]R Core Team, 2017).
Results
Differentially expressed RNAs in LUAD
Among the LUAD datasets, expression data for 19,669 mRNAs, 7,309
lncRNAs, and 1,882 miRNAs were extracted from TCGA and analyzed with R
package (DESeq2). In this analysis, we compared adjacent normal lung
tissue with lung cancer subtypes defined by pathology stage (T1-T4),
respectively. When we combined these four groups and analyzed for DE
RNAs, 2,160 mRNAs (1,527 up- and 1,083 down-regulated), 915 lncRNAs
(662 up- and 253 down-regulated), and 125 miRNAs (73 up- and 52
down-regulated) showed consistently differential expression (FC ≥ 2,
FDR < 0.01) ([78]Figs. 1A–[79]1C). [80]Figs. 2A–[81]2C shows the
heatmap of the common DEGs; the lung tumor tissue samples could be
easily distinguished from the adjacent non-tumor lung tissues from the
heat map. Based on these data, DE mRNAs, lncRNAs, and miRNAs were
selected for further analysis.
Figure 1. Significantly expressed genes in LUAD.
[82]Figure 1
[83]Open in a new tab
(A) The Venn diagram of mRNA DEGs, the circle represents mRNA DEGs
between normal tissue and tumor stage tissue; (B) the Venn diagram of
lncRNA DEGs; (C) the Venn diagram of miRNA DEGs.
Figure 2. Heatmap of significantly expressed common genes between different
pathology stages of LUAD.
[84]Figure 2
[85]Open in a new tab
(A) mRNA DEGs; (B) lncRNA DEGs; (C) miRNA DEGs.
Functional enrichment of mRNA DEGs
In order to understand the functions of aberrantly expressed genes in
this study, a total of 2,610 aberrantly expressed genes were analyzed.
We used GOseq for GO annotation and KOBAS 2.0, which includes the KEGG
database, to classify and analyze the potential gene functions in the
pathways. Our analysis revealed an enrichment of 2,879 GO terms and 30
KEGG pathways (p value < 0.05) ([86]Table S1). Gene ontology analyses
were classified into three functional groups: molecular function group,
biological process group, and cellular component group ([87]Fig. 3).
Enrichment of these DEGs represents a measure of significance of a
function. As shown in [88]Fig. 3, the most enriched GO biological
process terms were “single-organism cellular process” and
“single-multicellular organism process.” Other significant GO terms
included “cell proliferation,” “cell adhesion,” “system development,”
“cell migration,” and “cell adhesion.” The significant pathways
identified by KEGG pathway analysis included viral carcinogenesis,
systemic lupus erythematosus, cell adhesion molecules, and p53
signaling pathway ([89]Fig. 4). Our analysis indicated that 2,610
dysregulated mRNAs were involved in signaling pathways related to
environmental information processing and human diseases such as cancer.
Figure 3. Gene Ontology analysis and significant enriched GO terms of DEGs in
LUAD.
[90]Figure 3
[91]Open in a new tab
GO analysis classified the DEGs into three groups (i.e., molecular
function, biological process and cellular component).
Figure 4. Significantly enriched pathway terms in LUADs and 5ST.
[92]Figure 4
[93]Open in a new tab
The functional and signaling pathway enrichment were conducted using
KOBAS 2.0 software.
Competing endogenous RNA network analysis and integrated ceRNA network
construction
Based on the miRNAs and mRNAs described in [94]Fig. 2, starBase v2.0
database and miRanda were used to predict miRNA-targeted RNAs. Then,
combing miRNA–lncRNA interactions with the miRNA–mRNA interactions, an
integrated lncRNA–miRNA–mRNA ceRNA network was established. As shown in
[95]Fig. 5, the network contained 107 interactions, and included 5
lncRNAs, 20 mRNAs, and 34 miRNAs ([96]Table S2). After ranking the
molecules according to their degree, FENDRR, which is a lncRNA, showed
the highest degree (degree = 38). FOXF1, EDNRB, and EPAS1 are all mRNA
molecules, ranking the second, third, and fourth respectively. In other
studies, FENDRR was demonstrated to be down-regulated in breast and
gastric cancers and to participate in regulating tumor malignancy
([97]Li et al., 2018; [98]Xu et al., 2014). FOXF1 and EPAS1 are TFs
that have been demonstrated to participate in tumor progression
([99]Putra et al., 2015; [100]Tamura et al., 2014). EDNRB is a G
protein-coupled receptor that activates a phosphatidylinositol-calcium
second messenger system, and is also closely related with cancer
([101]Chen et al., 2013; [102]Wuttig et al., 2012).
Figure 5. The lncRNA–miRNA–mRNA ceRNA network constructed from DEGs.
[103]Figure 5
[104]Open in a new tab
Purple represents up-regulated genes and yellow represents
down-regulated genes. Circles represent mRNA, rectangles represent
lncRNA, and triangles represent miRNA.
Key ceRNAs and their association with clinical features
In order to further explore the key role of ceRNAs in the occurrence
and development of LUAD, we aimed to identify prognostic-specific
ceRNAs. We calculated expression of each ceRNA and performed prognosis
analyses (e.g. of overall survival time) to determine the
survival-significant lncRNA–miRNA–mRNA interactions based on Cox
analysis (p value < 0.05). As presented in [105]Fig. 6 and [106]Table
1, univariate Cox regression analysis showed that
EDNRB-FENDRR-hsa-mir-196b, EPAS1-FENDRR-hsa-mir-148a,
FOXF1-FENDRR-hsa-mir-148a, FOXF1-FENDRR-hsa-mir-195, and
FOXF1-FENDRR-hsa-mir-301b ceRNA interactions were associated with the
overall survival of patients with LUAD. Kaplan–Meier survival curves
indicated that EDNRB-FENDRR-hsa-mir-196b interaction was negatively
correlated with overall survival, whereas EPAS1-FENDRR-hsa-mir-148a,
FOXF1-FENDRR-hsa-mir-148a, FOXF1-FENDRR-hsa-mir-195, and
FOXF1-FENDRR-hsa-mir-301b were positively correlated with survival.
Figure 6. Kaplan–Meier survival curves for five ceRNA pairs associated with
overall survival.
[107]Figure 6
[108]Open in a new tab
(A) EDNRB-FENDRR-hsa-mir-196b. (B) EPAS1-FENDRR-hsa-mir-148a. (C)
FOXF1-FENDRR-hsa-mir-148a. (D) FOXF1-FENDRR-hsa-mir-195. (E)
FOXF1-FENDRR-hsa-mir-301b.
Table 1. Gene fold change in five ceRNA pairs.
T1_vs_Normal T2_vs_Normal T3_vs_Normal T3_vs_Normal
Gene log2FC FDR log2FC FDR log2FC FDR log2FC FDR
EPAS1 −2.64 5.65E-97 −2.76 1.1E-127 −2.49 1.3E-64 −2.61 9.45E-41
FOXF1 −2.44 6.68E-55 −2.96 2E-126 −2.58 1.2E-45 −2.63 8.73E-30
EDNRB −3.27 2.09E-78 −3.75 9.2E-123 −3.38 3.4E-60 −2.94 8.47E-38
FENDRR −3.48 2.85E-54 −4.37 2.7E-112 −3.78 1.62E-42 −3.89 1.94E-42
hsa-mir-148a 1.51 1.31E-24 1.46 1.51E-23 1.45 4.7E-11 1.51 1.07E-07
hsa-mir-195 −1.98 6.47E-37 −2.25 5.03E-68 −2.28 2.9E-30 −2.08 6.23E-13
hsa-mir-196b 3.00 1.59E-21 4.42 7.94E-35 2.98 5.6E-15 4.21 1.08E-20
hsa-mir-301b 3.48 2.62E-22 4.14 7.95E-28 2.94 9.5E-15 3.02 1.82E-05
[109]Open in a new tab
Transcription factor cis-target analysis
Transcription factors play important roles in the regulation of
biological processes by binding enhancer or promoter regions of DNA,
thereby activating or suppressing gene expression. This allows genes to
be expressed in the right cell at the appropriate time. In this study,
we found that EPAS1 and FOXF1 showed a high degree of representation in
the ceRNA network. This finding suggests that these two TFs play an
important role in LUAD and may be regulated by a ceRNA-based mechanism.
In order to verify the direct targets of the two TFs, we analyzed them
further. Co-expressed genes may be regulated by the same or similar
transcription factors; however, they may include numerous indirect
targets. Therefore, we performed binding motif enrichment analysis to
identify direct target genes of the two TFs. In this study, we use
GENIE3 to perform TF-genes co-expression analysis and RcisTarget to
analyze enriched motifs in the TSS upstream of the gene-sets. As shown
in [110]Fig. 7, the EPAS1 co-expressed gene-set was enriched for
the“cisbp__M6212” motif (NES = 3.21); likely direct gene targets of
EPAS1 include ADRB2, CALCRL, EMP2, FHL1, GRIA1, LGI3, LIMS2, NCKAP5,
RXFP1, SLC6A4, SMAD6, STX11, and TMEM100. However, the FOXF1
co-expressed gene-set did not show enrichment for any motifs.
Figure 7. Co-expressed genes analysis for EPAS1.
[111]Figure 7
[112]Open in a new tab
(A) EPAS1 co-expressed gene-set was enriched for motif “cisbp__M6212.”
(B) EPAS1-regulated genes in LUAD.
Validation of the expression of key RNA molecules with GEO data
Finally, four reported studies were screened out from the GEO to verify
the differential expression of key mRNAs and lncRNAs in LUAD.
[113]GSE10072 and [114]GSE32863 were used for EPAS1, FOXF1, and EDNRB
mRNA verification. [115]GSE85716 and [116]GSE104854 were used for
FENDRR differential expression. As shown in [117]Fig. 8 and [118]Table
2, all four molecules showed significantly lower levels of expression
in tumor datasets, which is consistent with the results of our study.
Figure 8. Validation of gene expression in GEO dataset.
[119]Figure 8
[120]Open in a new tab
(A) EPAS1, FOXF1, and EDNRB expression level in [121]GSE10072. (B)
EPAS1, FOXF1 and EDNRB expression level in [122]GSE32863. (C) FENDRR
expression level in [123]GSE85716 and [124]GSE104854.
Table 2. EPAS1, FOXF1, EDNRB and FENDRR expression in GEO dataset.
Study mRNA lncRNA
EPAS1 FOXF1 EDNRB FENDRR
log2FC FDR log2FC FDR log2FC FDR log2FC FDR
[125]GSE10072 −2.087 3.08E-33 −2.228 2.17E-35 −2.299 4.30E-35 – –
[126]GSE32863 −2.152 3.05E-33 −1.616 2.69E-28 −2.385 2.35E-44 – –
[127]GSE85716 – – – – – – −4.25009 0.000343
[128]GSE104854 – – – – – – −4.12757 3.96E-08
[129]Open in a new tab
Discussion
Non-small cell lung cancer, which mainly includes adenocarcinoma and
squamous cell carcinoma, accounts for about one-sixth of cancer deaths
in the global population ([130]Jemal et al., 2011; [131]Torre et al.,
2015). Despite numerous advances in the treatment of lung cancer, the
prognosis of this disease remains poor, and the 5-year overall survival
rate is only 11% ([132]Huang et al., 2012; [133]Miura et al., 1996).
Therefore, elucidation of the underlying mechanisms is indispensable
for diagnosing, preventing, and treating NSCLC.
With the development of high-throughput sequencing technology,
accumulating omics data reveal that the occurrence and development of
the disease is closely related to multiple factors, including genomic
mutations, expression profiles, copy number variants, and structural
variations ([134]Balbin et al., 2013; [135]Kikutake & Yahara, 2016;
[136]Kim et al., 2015; [137]Yan et al., 2017). In recent years, several
studies have found that multiple RNA molecules are closely related to
tumor progression ([138]Chen et al., 2018; [139]He et al., 2018;
[140]Men, Liu & Ren, 2018; [141]Shang et al., 2017; [142]Song et al.,
2018; [143]Tang et al., 2017; [144]Wang et al., 2018; [145]Yu et al.,
2017). For example, lncRNA-DANCR enhances the stemness feature of
cancer cells to increase the degree of malignancy of lung cancer
([146]Lu et al., 2018; [147]Zhen et al., 2018); miRNA-221 was reported
to regulate the process of tumor differentiation, proliferation,
migration, and invasion ([148]Lv et al., 2014; [149]Yin, Xu & Li,
2017). These RNA molecules represent biological markers for the
prognosis of lung cancer. In this study, we found that multiple novel
ceRNA molecules can be used as novel prognostic markers; further, we
examined the putative molecular mechanisms underlying the development
of tumors.
miRNAs modulate the expression of target genes by regulating the
transcription and stability of their mRNAs or lncRNAs. Investigation of
the lncRNA-mRNA co-expression network is important for analysis of the
function and regulatory mechanisms of lncRNAs ([150]Fang, Wang & Li,
2018). In this study, analysis showed that a number of mRNAs, lncRNAs,
and miRNAs that were DE were common at each stage of lung cancer. Gene
ontology and KEGG analysis showed that the DEGs were mainly involved in
“single-organism cellular process,” “single-multicellular organism
process,” “cell proliferation,” “cell adhesion,” “system development,”
“cell migration,” and “cell adhesion biological processes.” Further
analysis revealed that the up-regulated genes were involved in the
mitotic cell cycle and nuclear division process ([151]Table S3), while
down-regulated genes were involved in vasculature development, cell
motility, cell migration, and cell proliferation ([152]Table S4). Kyoto
encyclopedia of genes and genomes analysis found that the up-regulated
genes were mainly involved in cancer and other related signaling
pathways such as the p53 signaling pathway.
Survival curve analysis reveals the correlation between gene expression
and patient prognosis. Our study found that the five most important
ceRNA-related RNA molecules were significantly correlated with the
survival rate of patients with lung cancer (p < 0.01), and most of the
genes were highly correlated with cancer. The five ceRNA interactions
consisted of one lncRNA, three mRNAs, and four miRNAs. The lncRNA,
FENDRR ([153]Grote et al., 2013), which is produced from a spliced long
non-coding RNA transcribed bidirectionally with FOXF1 on the opposite
strand, is down-regulated in both breast cancer and gastric cancer. In
breast cancer, lower expression level of FENDRR was associated with
shorter overall survival and progression-free survival in breast cancer
patients. FENDRR knockdown promoted breast cancer cell proliferation
and migration, and suppressed cell apoptosis; in contrast, its
overexpression inhibited tumor growth in a xenograft model ([154]Li et
al., 2018). FENDRR was also reported to be down-regulated in gastric
cancer cell lines and tissues. Further, FENDRR expression was
negatively correlated with tumor metastasis and stage ([155]Xu et al.,
2014). In our study, FENDRR was also down-regulated in LUAD, and showed
the highest degree in ceRNA network, indicating its important role in
LUAD. The three mRNAs, FOXF1, EDNRB, and EPAS1, ranked 2nd, 3rd, and
4th, respectively in the ceRNA network. FOXF1 and EPAS1 are both TFs
that have been associated with tumor malignancy ([156]Putra et al.,
2015; [157]Tamura et al., 2014). Our TF analysis indicated that EPAS1
directly regulates 13 molecules associated with LUAD, namely ADRB2,
CALCRL, EMP2, FHL1, GRIA1, LGI3, LIMS2, NCKAP5, RXFP1, SLC6A4, SMAD6,
STX11, and TMEM100. EDNRB, a G protein-coupled receptor, is also
closely related with tumor ([158]Chen et al., 2013; [159]Wuttig et al.,
2012). Our findings suggest that these molecules represent potential
cancer-related biomarkers.
Conclusions
In summary, in the present study, we comprehensively analyzed the
expression of key RNA molecules in the progression of lung cancer and
identified key RNA molecules (including mRNA/lncRNA/miRNA) that were
abnormally expressed in different stages of this disease. Further,
using the above-mentioned key RNA molecules to construct a ceRNA–gene
interaction network, it was found that several FOXF1- and EDNRB-related
ceRNA molecules play an important role in the development and
progression of lung cancer, and can affect the prognosis of this
disease. This series of related ceRNA molecules is expected to provide
a novel basis for diagnosis and evaluation of prognosis, and represents
potential new targets for the treatment of lung cancer.
Supplemental Information
Supplemental Information 1. Table S1. GO and KEGG enrichment analysis
of 2,610 dysregulated mRNAs.
[160]Click here for additional data file.^ (924.8KB, xlsx)
DOI: 10.7717/peerj.6694/supp-1
Supplemental Information 2. Table S2. Interactions of the ceRNA network
inLUAD.
[161]Click here for additional data file.^ (12.7KB, xlsx)
DOI: 10.7717/peerj.6694/supp-2
Supplemental Information 3. Table S3. GO and KEGG enrichment analysis
of up-regulated mRNAs.
[162]Click here for additional data file.^ (388.6KB, xlsx)
DOI: 10.7717/peerj.6694/supp-3
Supplemental Information 4. Table S4. GO and KEGG enrichment analysis
of down-regulated mRNAs.
[163]Click here for additional data file.^ (485.3KB, xlsx)
DOI: 10.7717/peerj.6694/supp-4
Supplemental Information 5. Clinical characteristics of patients on
LUAD of TCGA.
[164]Click here for additional data file.^ (132KB, zip)
DOI: 10.7717/peerj.6694/supp-5
Supplemental Information 6. Raw expression data of miRNA on LUAD of
TCGA.
[165]Click here for additional data file.^ (2.3MB, zip)
DOI: 10.7717/peerj.6694/supp-6
Supplemental Information 7. Raw expression data of mRNA on LUAD of
TCGA.
[166]Click here for additional data file.^ (17.3MB, zip)
DOI: 10.7717/peerj.6694/supp-7
Supplemental Information 8. Raw expression data of lncRNA on LUAD of
TCGA.
[167]Click here for additional data file.^ (15.5MB, zip)
DOI: 10.7717/peerj.6694/supp-8
Acknowledgments