Characterization of mRNA Expression and Endogenous RNA Profiles in Bladder Cancer Based on The Cancer Genome Atlas (TCGA) Database.

BACKGROUND Bladder cancer is a multifactorial disease with increasing incidence and mortality. Genetic alterations and altered expressions of mRNAs, long non-coding RNAs (lncRNAs), and miRNAs have been shown to play important roles in the tumorigenesis of bladder cancer. However, the functions of key RNAs and their regulatory network in bladder cancer are still to be elucidated. MATERIAL AND METHODS RNA profiles were downloaded from The Cancer Genome Atlas (TCGA) database. The differentially expressed mRNAs, lncRNAs, and miRNAs in bladder cancer were acquired through analyses of data from 414 bladder cancer tissues and 19 normal bladder tissues. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis was performed by using "DAVID6.8" and the R package "ClusterProfile". Protein-protein interaction and competing endogenous RNA (ceRNA) networks were constructed by using "STRING" database and Cytoscape 3.6.2. Based on the clinical data and Cox regression, a prognosis model was established, and survival analysis was performed. RESULTS A total of 1819 mRNAs, 659 lncRNAs, and 160 miRNAs were identified as significantly differentially expressed in bladder cancer of which 52 mRNAs, 58 lncRNAs, and 22 miRNAs were incorporated in the ceRNA network. CFL2 and TPM2 were found to be downregulated and showed significant correlation to each other in bladder cancer. HOXB5 and 6 lncRNAs (ADAMTS9-AS1, AC112721.1, LINC00460, AC110491.1, LINC00163, and HCG22) were strongly associated with high-grade, disease stages, and overall survival. CONCLUSIONS In this study, we have identified differentially expressed mRNAs, lncRNAs, and miRNAs in bladder cancer which were strongly associated with oncogenesis and prognosis. Further experimental studies are necessary to validate these results.

Background

Bladder cancer is one of the most malignant tumors with high incidence and mortality during the past decades [1]. In 2018, the incidence of bladder cancer was estimated to be approximately 81 000 new cases and the mortality to be 17 000 deaths in the United States (ranks sixth and eighth, respectively) [2]. Histologically, approximately 95% of bladder cancer occurs in the urothelial cells, and around three-quarters of bladder cancer patients are non-muscle invasive bladder cancer while the rest are muscle invasive bladder cancer [3-5]. The grading of bladder cancer is crucial, and the World Health Organization (WHO) 2004/2016 grading system is preferred by most pathologists. Based on the cytologic and architectural abnormalities, atypia, low-grade, and high-grade determine the degree of malignancy and prognosis of bladder cancer [6,7]. Although new strategies for cancer diagnosis, detection, and neo-treatments have been developed for patients with bladder cancer in the past 30 years, the 5-year relative survival rate is still unsatisfactory [8]. Further research is thus essential to better understand the underlying pathogenetic mechanisms to improve the outcome of bladder cancer.

Recently, genetic characteristics and acquired genetic alterations have been found to play an important role in the bladder cancer. And these distinct mechanisms have been found to construct a multifocal network. For example, inactivating mutations in the tumor suppressor TP53 and activating mutations in FGFR3 are found in both papillary and non-papillary bladder cancer [9,10]. Moreover, mutations in genes encoding transcription factors and chromatin-modifying enzymes and mutations in TERT promoter are also strongly implicated in some cases [11-13].

With the development of high-throughput sequencing methods, thousands of pseudogenes have been characterized [14]. Long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) have been revealed to play a critical role in different kinds of cancer including bladder cancer, by modulating and modifying their ancestral gene. Our research group previously reported that lncRNA ROR was significantly increased in bladder cancer and positively associated with its potential targeting gene ZBE1, thereby contributing to the progression of bladder cancer as well as promoting epithelial-to-mesenchymal transition (EMT) [15]. Liu et al. demonstrated that lncRNA SRY4-IT1 sponges mir-101-3p and upregulates EZH2 leading to aggressive phenotypes in bladder cancer [16]. Therefore, characterization of an integrated whole differential gene expressions network with related endogenous RNA profiles is important, as it might play a key role in the pathogenesis of bladder cancer.

In this study, we identified the differentially expressed mRNA, lncRNA, and miRNA expression profiles in bladder cancer from the TCGA database. In addition, we performed functional analyses of the differentially expressed RNAs and investigated their clinical significance in relation to prognosis in bladder cancer.

Material and Methods

Patients and pathological data

Transcriptome profiling data of 414 bladder cancer tissues and 19 normal bladder tissues were acquired from The Cancer Genome Atlas (TCGA) in October 2018. The RNA-seq data were generated from the Illumina HiSeqRNASeq and Illumina HiSeqmiRNASeq platforms. Using the GDC Data Transfer Tool (), all the gene expression profiles and clinical data of bladder cancer were downloaded. Ethical consent was not required as all the data in this study were obtained from TCGA database.

Identification of differentially expressed RNA

The “DESeq” package in R software [17] was utilized to identify the differentially expressed RNAs in bladder cancer when compared to normal bladder tissues. The significance level of the adjusted P-value was set at P<0.01 and the thresholds was set as |log2FoldChange| >2. Moreover, the differentially expressed mRNAs, lncRNAs and miRNAs were annotated using ENSEMBL (htps://).

CeRNA network construction

According to the ceRNA theory, lncRNAs can act as endogenous RNA and thereby regulating target gene transcripts by competing with shared miRNAs[18]. Based on the differentially expressed mRNAs, lncRNAs and miRNAs in bladder cancer, the target mRNAs of the miRNAs were predicted by Targetscan (), miRDB () and miRTarBase () [19]. The miRanda database () was used for the lncRNAs and miRNAs target predictions [20]. In addition, by combining the discriminatory expression profiles data, the interaction between lncRNAs and mRNAs was identified. The ceRNA network of mRNA-lncRNA-miRNA was visualized with Cytoscape v3.6.2.

Functional enrichment analysis

Gene Ontology (GO) enrichment analysis was performed by (DAVID 6.8) database () to group the differentially expressed RNAs into 3 categories including molecular function, biological process, and cellular component. The setting in the GO analysis was a false discovery rate (FDR) <0.01. The “ClusterProfiler” package in R software was utilized to analyze the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The barplots for GO and KEGG were generated by the “Goplot” package.

Prognosis risk scoring in differentially expressed RNAs

Univariate Cox regression analysis was conducted where P<0.0005 was considered statistically significant. The significant mRNAs in univariate Cox regression were subsequently analyzed in a multivariate Cox regression proportional hazards model. Furthermore, based on the median risk score, bladder cancer patients were divided into “high-risk” and “low-risk” groups. The risk scoring system was constructed by the formula as follows for predicting overall survival (OS): risk score=βg 1×exprg 1+βg 2×exprg 2+::: +βg n×exprg n [21].

Kaplan-Meier survival analysis and receiver operating characteristic (ROC) analysis were then performed to assess the risk scoring system with high and low scores. The survival and ROC analyses were accomplished by using the R package “survival” and “survivalROC”.

Protein–protein interaction (PPI) and correlation network construction

The differentially expressed mRNAs in the ceRNA network were analyzed through the protein–protein interaction (PPI) network. The PPI network was constructed with the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING). The minimum interaction value was set at low confidence (0.150). Based on the top 10 combined prediction score, we selected these gene expression profiles to perform correlation analyses for further validation.

Genes associated with grading and staging in ceRNA network and survival analysis

Differentially expressed mRNA, lncRNA, and miRNA in the ceRNA network were analyzed in relation to low-grade and high-grade bladder cancer. These aberrantly expressed RNAs were also evaluated in relation to tumor stage. The significance level of adjusted P-value was <0.01 and the thresholds was set as |log2FoldChange| >2. Kaplan-Meier plots and log-rank test was performed to evaluate the genes in relation to overall survival, where P<0.05 was considered as statistically significant.

Results

Characteristics of clinical features in TCGA database

Clinical information of bladder cancer patients from the TCGA database was available in 412 of 433 cases., The clinical characteristics including age, gender, tumor grade, tumor stage, and lymph node metastasis of the bladder cancer patients are presented in Table 1. Interestingly, except for gender, all other clinical features were significantly associated with survival. Patients with high-grade bladder cancer were associated with worse outcomes compared to patients with low-grade bladder cancer. The median age of the cohort was 69 years (range: 34–90 years).

Table 1

The clinicopathological characteristics of bladder cancer patients.

Characteristics	Subtype	Patinets n (%)	Dead	χ2	P-value
Age	>65	250 (60.7)	115	14.722	<0.01
	≤65	162 (39.3)	44
Gender	Male	304 (73.8%)	115	0.285	0.593
	Female	108 (26.2%)	44
Tumor grade	High grade	388 (94.2%)	159	16.016	<0.01
	Low grade	24 (5.8%)	0
Tumor Stage	Stage I and II	133 (32.3%)	27	27.728	<0.01
	Stage III and IV	279 (67.7%)	132
Lymph node	Negative	239 (58.0%)	63	37.810	<0.01
	Positive	131 (31.8%)	77
	Unknow	42 (10.2%)	19

N=412.

Differentially expressed RNAs and ceRNA network in bladder cancer

A total of 1819 mRNAs, 659 lncRNAs, and 160 miRNAs were found to be differentially expressed in bladder cancer, of which 1030 mRNAs (56.7%), 415 lncRNAs (63.0%) and 139 miRNAs (86.9%) were upregulated while others were downregulated. The differentially expressed RNAs are visualized in the volcano plot (Figure 1). We further validated the relationships and functions of these differentially expressed RNAs according to the ceRNA hypothesis. The miRNA-mRNA and lncRNA-miRNA interactions were first predicted to find 52 mRNAs targeted by 17 key miRNAs involved in ceRNAs network (Table 2) and 58 lncRNAs were predicted to be interacted with 22 miRNAs (Table 3). As a result, 52 mRNAs, 58 lncRNAs, and 22 miRNAs of the differentially expressed RNAs were incorporated in the established ceRNA network (Tables 4–6). In this ceRNA network, 16 mRNAs, 35 lncRNAs, and 17 miRNAs were downregulated; while 36 mRNAs, 23 lncRNAs, and 5 miRNAs were found to be upregulated (Figure 2).

Figure 1

Differentially expressed RNAs in bladder cancer are visualized by volcano plots. The red dots show upregulated genes while the green dots show downregulated genes. (A) Differentially expressed mRNAs in bladder cancer. (B) Differentially expressed lncRNAs in bladder cancer. (C) Differentially expressed miRNAs in bladder cancer.

Table 2

Differentially expressed miRNAs targeting mRNAs in the ceRNA network.

miRNA	mRNA
miR-141	EPHA7; ZEB1; ELAVL2; HOXB5
miR-145	MEST
miR-182	FGF9; THBS1; TCEAL7; ULBP2; PRKAA2
miR-183	AKAP12; ZEB1; CCNB1
miR-195	BTG2; TPM2; RUNX1T1; FGF2; TGFBR3; CBX2; PRICKLE2; MYB; WNT7A; MKX; ALOX12; ITPR1; CCNE1; RAB23; TMEM100; E2F7
miR-200a	HOXB5; EPHA7; ZEB1; ELAVL2; CCNE2
miR-205	ZEB1; SHISA6; LRRK2
miR-210	SERTM1; AIFM3; NR4A2
miR-217	NR4A2; MAP1B
miR-31	SELE; HOXC13
miR-372	DUSP2; CADM2; TMEM100; FBXL7; ELAVL2
miR-373	CADM2; TMEM100; FBXL7; ELAVL2; CFL2; DUSP2
miR-383	DIO1
miR-429	JUN; ZEB1; ZFPM2
miR-503	GREM2
miR-519d	NACC2; HMGB3; CYBRD1; SALL3; CFL2; FAM129A; DUSP2; POLQ; ELAVL2
miR-96	SLC25A25; ZEB1

Table 3

Differentially expressed miRNAs targeting lncRNAs in ceRNA network.

LncRNA	miRNA
IGF2-AS	miR-519d; miR-503
LINC00525	miR-301b; miR-96; miR-141; miR-200a; miR-182; miR-31; miR-383
PART1	miR-301b; miR-141; miR-200a; miR-143; miR-145; miR-195; miR-429; miR-205; miR-31
AC009065.1	miR-372; miR-373
C20orf166-AS1	miR-301b; miR-372; miR-373; miR-519d; miR-183; miR-429; miR-205; miR-489
GRIK1-AS1	miR-145; miR-205; miR-383
C20orf197	miR-372; miR-373; miR-143; miR-519d; miR-383
AP002478.1	miR-503; miR-372; miR-373; miR-195; miR-519d; miR-182; miR-192; miR-215; miR-205; miR-489
LINC00518	miR-141; miR-200a; miR-143; miR-145
LINC00482	miR-143
MIR22HG	miR-383; miR-489
C9orf163	miR-143; miR-195; miR-205; miR-489
LINC00336	miR-96; miR-143; miR-145; miR-217
AC008676.1	miR-301b; miR-372; miR-373; miR-141; miR-200a; miR-143; miR-519d; miR-31
LINC00487	miR-372; miR-373; miR-143; miR-183; miR-205; miR-31
AP000525.1	miR-503; miR-31
AC104472.1	miR-143; miR-183; miR-429; miR-31; miR-489
AL513123.1	miR-141; miR-200a; miR-183
LINC00473	miR-145; miR-195; miR-210
LINC00337	miR-372; miR-373; miR-145; miR-519d; miR-182; miR-217; miR-383
AC127496.3	miR-301b; miR-372; miR-373; miR-145; miR-183; miR-192; miR-215; miR-429
LINC00161	miR-145; miR-205
HCG22	miR-96; miR-145; miR-195; miR-182; miR-31; miR-383; miR-489
SACS-AS1	miR-503; miR-372; miR-143; miR-205
MIR137HG	miR-182; miR-192; miR-215;
RASA3-IT1	miR-205; miR-383
NALCN-AS1	miR-372; miR-373; miR-195; miR-182; miR-205; miR-31; miR-383
ERVH48-1	miR-301b; miR-96; miR-141; miR-200a; miR-145; miR-182
LINC00472	miR-503; miR-372; miR-373; miR-141; miR-200a; miR-143; miR-145; miR-195; miR-383; miR-489
AC110491.1	miR-141; miR-200a; miR-143; miR-182; miR-192; miR-215; miR-429; miR-205; miR-489
TLR8-AS1	miR-182; miR-31
LINC00460	miR-503; miR-143; miR-429; miR-489
JAZF1-AS1	miR-372; miR-373; miR-143; miR-519d; miR-205
MAGI2-AS3	miR-503; miR-372; miR-373; miR-141; miR-200a; miR-143; miR-145; miR-195; miR-429; miR-210; miR-217; miR-31; miR-489
LINC00163	miR-143; miR-183; miR-210
LINC00330	miR-503; miR-301b; miR-372; miR-373; miR-145; miR-195; miR-519d; miR-192; miR-215; miR-205; miR-383
LINC00402	miR-141; miR-200a; miR-143; miR-519d; miR-182; miR-429; miR-217; miR-383
TM4SF19-AS1	miR-141; miR-200a; miR-205
MYO16-AS1	miR-489
DLEU7-AS1	miR-96; miR-195; miR-182; miR-192; miR-215
AC009121.1	miR-141; miR200a
AC112721.1	miR-503; miR-195
AP004609.1	miR-383
ADAMTS9-AS1	miR-301b; miR-96; miR-145; miR-182; miR-31
HNF1A-AS1	miR-372; miR-373; miR-141; miR-200a; miR-195; miR-519d; miR-183; miR-217
ADAMTS9-AS2	miR-301b; miR-372; miR-373; miR-96; miR-141; miR-200a; miR-143; miR-145; miR-182; miR-183; miR-205; miR-31
AC078778.1	miR-301b
AC012640.1	miR-182; miR-205; miR-210; miR-383
AC073352.1	miR-96; miR-182
AC011453.1	miR-143; miR-205
AP000553.1	miR-192; miR-215; miR-217
AC128709.1	miR-503; miR-195; miR-183
LINC00534	miR-372; miR-373; miR-96; miR-192; miR-215; miR-205; miR-217; miR-489
NAV2-AS2	miR-96; miR-182; miR-31
ARAP1-AS1	miR-145
FRMD6-AS2	miR-143; miR-182
LINC00520	miR-503; miR-372; miR-373; miR-145; miR-195; miR-519d; miR-205; miR-217; miR-31

Table 4

Differentially expressed mRNAs in ceRNA network.

DEmRNAs	Regulation	Log fold change	FDR
CFL2	Down-regulation	−3.18787341	6.81E-64
SLC25A25	Down-regulation	−2.656571564	1.81E-63
FAM129A	Down-regulation	−3.684416254	2.52E-60
NACC2	Down-regulation	−2.417159742	5.74E-48
RAB23	Down-regulation	−2.902696027	1.86E-43
ITPR1	Down-regulation	−2.677318488	4.37E-38
ZEB1	Down-regulation	−2.830314677	2.64E-37
MAP1B	Down-regulation	−3.224125414	1.18E-34
NR4A2	Down-regulation	−2.806554541	9.69E-33
TPM2	Down-regulation	−3.056564907	1.34E-30
FBXL7	Down-regulation	−2.52708691	5.43E-28
THBS1	Down-regulation	−2.686696518	7.25E-25
PRICKLE2	Down-regulation	−2.241943352	9.47E-25
JUN	Down-regulation	−2.081270648	1.94E-23
ZFPM2	Down-regulation	−2.664653755	3.11E-23
AKAP12	Down-regulation	−2.592811618	8.39E-22
TMEM100	Down-regulation	−3.095618096	4.25E-20
TCEAL7	Down-regulation	−2.49946	5.49E-20
EPHA7	Down-regulation	−3.166	9.53E-20
FGF2	Down-regulation	−2.69593	9.99E-19
RUNX1T1	Down-regulation	−2.57174	1.52E-18
LRRK2	Down-regulation	−2.36765	2.55E-18
BTG2	Down-regulation	−2.13366	6.79E-18
CYBRD1	Down-regulation	−2.17256	8.81E-18
MKX	Down-regulation	−3.12218	1.49E-17
HMGB3	Up-regulation	2.057653	1.33E-16
PRKAA2	Down-regulation	−2.63399	1.31E-15
POLQ	Up-regulation	2.442571	1.38E-15
TGFBR3	Down-regulation	−2.0693	2.90E-15
DUSP2	Down-regulation	−2.49447	8.28E-15
GREM2	Down-regulation	−3.03523	2.37E-14
CCNB1	Up-regulation	2.066358	1.53E-12
SERTM1	Down-regulation	−3.83855	2.49E-12
SHISA6	Down-regulation	−2.77251	3.62E-11
CCNE1	Up-regulation	2.305447	2.37E-10
FGF9	Up-regulation	2.519249	2.45E-10
SELE	Down-regulation	−2.78825	7.82E-10
E2F7	Down-regulation	−2.47734	1.80E-09
CCNE2	Up-regulation	2.254535	6.75E-09
ALOX12	Up-regulation	2.116593	6.91E-09
ULBP2	Down-regulation	−2.00117	1.01E-08
MEST	Up-regulation	2.464683	7.47E-08
CBX2	Up-regulation	2.595458	1.57E-06
CADM2	Up-regulation	2.183717	3.23E-06
HOXC13	Down-regulation	−2.25318	2.10E-05
MYB	Up-regulation	2.796581	4.20E-05
SALL3	Up-regulation	2.034968	4.76E-05
AIFM3	Down-regulation	−2.55675	0.000173
HOXB5	Up-regulation	2.147751	0.000211
WNT7A	Up-regulation	2.024795	0.000334
ELAVL2	Up-regulation	6.532413	0.000532
DIO1	Up-regulation	3.143323	0.000581

Table 5

Differentially expressed lncRNAs in ceRNA network.

DELncRNAs	Regulation	Log fold change	FDR
HCG22	Down-regulation	−7.393236357	8.70E-86
ADAMTS9-AS1	Down-regulation	−5.173460141	1.28E-59
ADAMTS9-AS2	Down-regulation	−4.152710881	3.17E-50
LINC00330	Down-regulation	−5.357713971	8.39E-35
C20orf166-AS1	Down-regulation	−4.418702998	7.39E-34
MIR22HG	Down-regulation	−2.186866609	2.83E-29
JAZF1-AS1	Down-regulation	−3.248756128	2.20E-26
AC008676.1	Down-regulation	−2.677691052	1.61E-22
RASA3-IT1	Down-regulation	−4.319901412	8.29E-20
MAGI2-AS3	Down-regulation	−2.271922028	1.15E-18
FRMD6-AS2	Down-regulation	−3.536986684	6.70E-18
PART1	Down-regulation	−2.883633316	5.25E-16
AC110491.1	Down-regulation	−3.809886233	5.66E-16
AC078778.1	Up-regulation	2.171874229	2.62E-14
SACS-AS1	Down-regulation	−3.387765888	1.98E-13
AP004609.1	Down-regulation	−2.579281355	2.67E-11
C9orf163	Up-regulation	2.339600002	6.13E-11
LINC00472	Down-regulation	−2.086567223	5.03E-10
AC127496.3	Down-regulation	−2.916229805	7.40E-10
GRIK1-AS1	Down-regulation	−2.185078089	8.15E-10
AC073352.1	Up-regulation	2.540191221	4.51E-09
NALCN-AS1	Down-regulation	−2.49519035	8.70E-09
AP000553.1	Up-regulation	3.216908948	1.44E-08
LINC00163	Down-regulation	−2.829231994	1.25E-07
AC009121.1	Up-regulation	2.162583137	1.08E-06
LINC00460	Up-regulation	7.163558998	1.22E-06
DLEU7-AS1	Up-regulation	2.139118702	3.54E-06
LINC00337	Up-regulation	2.018174432	4.36E-06
AP000525.1	Up-regulation	2.750867634	1.26E-05
AL513123.1	Up-regulation	3.339752428	1.64E-05
LINC00161	Down-regulation	−2.004179345	2.13E-05
LINC00402	Down-regulation	−2.228396776	2.28E-05
C2orf48	Up-regulation	2.234385723	2.41E-05
AC112721.1	Up-regulation	5.679017823	2.64E-05
ERVH48-1	Up-regulation	6.171777638	3.07E-05
TM4SF19-AS1	Up-regulation	2.241213677	5.71E-05
LINC00525	Up-regulation	2.388170307	7.58E-05
AP002478.1	Up-regulation	2.499724256	7.88E-05
AC011453.1	Up-regulation	4.106228298	0.000145979
C20orf197	Up-regulation	3.222909721	0.000173251
LINC00487	Up-regulation	2.684913886	0.00019263
AC012640.1	Up-regulation	3.267430057	0.000198949
AC009065.1	Up-regulation	2.029597356	0.000245796
LINC00518	Up-regulation	4.436236799	0.000596273
AC104472.1	Up-regulation	2.743663788	0.001205612
MYO16-AS1	Up-regulation	4.783639374	0.001480169
HNF1A-AS1	Up-regulation	3.838602554	0.001551329
LINC00534	Up-regulation	2.943676684	0.002062953
LINC00520	Up-regulation	3.918845469	0.002229151
LINC00482	Up-regulation	2.211556791	0.003470315
LINC00473	Down-regulation	−2.010751587	0.003738568
LINC00336	Up-regulation	2.001908125	0.003803458
TLR8-AS1	Up-regulation	4.005367912	0.003995111
MIR137HG	Up-regulation	5.171531584	0.004865838
NAV2-AS2	Up-regulation	3.18541719	0.006297113
ARAP1-AS1	Up-regulation	3.214443621	0.007268128
AC128709.1	Up-regulation	2.114272765	0.008604833
IGF2-AS	Up-regulation	2.399553974	0.008980377

Table 6

Differentially expressed lncRNAs in ceRNA network.

DEmiRNAs	Regulation	Log fold change	FDR
hsa-mir-143	Down-regulation	−3.525679542	1.28E-41
hsa-mir-195	Down-regulation	−2.410426303	1.25E-32
hsa-mir-96	Up-regulation	3.767149943	8.97E-18
hsa-mir-210	Up-regulation	4.973834298	1.33E-17
hsa-mir-3199-2	Down-regulation	−2.187131639	1.55E-15
hsa-mir-145	Down-regulation	−2.133372137	3.72E-14
hsa-mir-183	Up-regulation	3.03501324	8.81E-14
hsa-mir-301b	Up-regulation	3.971849748	1.19E-12
hsa-mir-141	Up-regulation	2.631026555	3.23E-11
hsa-mir-503	Up-regulation	2.473934298	3.22E-10
hsa-mir-182	Up-regulation	2.348870325	3.51E-10
hsa-mir-429	Up-regulation	2.758180811	6.74E-09
hsa-mir-383	Down-regulation	−2.801221137	3.16E-08
hsa-mir-200a	Up-regulation	2.219506844	2.94E-07
hsa-mir-192	Up-regulation	2.40248191	6.49E-07
hsa-mir-205	Up-regulation	2.239011007	8.69E-06
hsa-mir-215	Up-regulation	4.242318245	2.68E-05
hsa-mir-31	Up-regulation	2.786222534	3.19E-05
hsa-mir-3136	Up-regulation	2.306880974	4.63E-05
hsa-mir-519d	Up-regulation	7.14700051	0.000164482
hsa-mir-372	Up-regulation	5.986379084	0.000186072
hsa-mir-489	Up-regulation	2.962732727	0.000303842

Figure 2

The ceRNA network of differentially expressed RNA profiles in bladder cancer, the upregulated genes are presented using red color while downregulated genes are shown using green color. Circles represent mRNAs, quadrangles represent lncRNAs and rhombus represent miRNAs.

GO and KEGG functional analysis

To further elucidate the functions and biological processes of the differentially expressed mRNAs, we conducted the GO and KEGG enrichment analysis. The barplot of the GO analysis showed significant enrichment in 26 molecular functions, biological processes, and cellular components (P-value <0.001) (Figure 3A). The results revealed that 313 genes were aggregated in extracellular region with the lowest P-value at 1.15E-36. KEGG pathway analysis showed that 20 pathways were associated with tumorigenesis of bladder cancer. As shown in Figure 3B, the neuroactive ligand-receptor interaction reflects the most significant cancer associated pathway which contained 61 genes, followed by alcoholism with 59 genes enriched.

Figure 3

GO and KEGG analysis of differentially expressed mRNAs in bladder cancer. (A) Twenty-six significant molecular functions, biological processes, and cellular components of differentially expressed mRNAs. (B) Twenty enrichment of KEGG pathways for differentially expressed mRNAs. GO – Gene Ontology; KEGG – Kyoto Encyclopedia of Genes and Genomes.

Prognostic risk score based on differentially expressed mRNAs

Of the 1820 differentially expressed mRNAs, we found 36 mRNAs to be significantly associated with survival using univariate Cox regression (Table 7). Subsequently, multivariate Cox regression was applied, and 12 mRNA expression profiles were identified as coupled to overall survival and included in the overall survival prediction model (Figure 4A, 4B). On the basis of the overall survival prediction model, a prognostic risk model was constructed where the high-risk group was significantly associated with lower overall survival (P<0.05) (Figure 4B). The discriminative evaluation of the risk scoring system was performed with ROC-curve analysis where we found an AUC=0.735 (Figure 4C).

Table 7

Significant genes demonstrated in univariate Cox regression model (P<0.005).

Gene	HR	z-score	P-value
EMP1	1.328762282	5.254914697	1.48E-07
MXRA7	1.432626177	4.705964739	2.53E-06
SERPINB12	1.162296991	4.580143916	4.65E-06
MAP1B	1.229296929	4.37622915	1.21E-05
LAMA2	1.267134732	4.327829003	1.51E-05
MAP1A	1.233776119	4.250957967	2.13E-05
KANK4	1.113571638	4.182201609	2.89E-05
ABCC9	1.231009469	4.178191876	2.94E-05
PEX5L	1.162984335	4.138450231	3.50E-05
PTGER3	1.144526818	4.026382546	5.66E-05
PTPRR	0.898457289	−3.961993355	7.43E-05
POU5F1	0.857411853	−3.875204799	0.000106535
PCSK9	1.102595243	3.853707011	0.000116343
TCHH	1.135841746	3.841592708	0.000122239
SAPCD1	0.791832106	−3.813548024	0.000136986
GABRG1	1.312711858	3.80120357	0.000143995
CNTN1	1.105824108	3.792658288	0.000149043
ADCYAP1R1	1.159131884	3.76366481	0.000167441
CCDC80	1.145437466	3.757189589	0.000171832
DTNA	1.138932514	3.75592982	0.000172699
PID1	1.179886946	3.740557523	0.000183613
RBMS3	1.203510811	3.693371944	0.0002213
XAGE2	1.105072802	3.689802104	0.000224429
PLEKHG4B	1.096473367	3.681482901	0.000231881
NTNG1	1.118858872	3.676085951	0.00023684
FBN2	1.108316893	3.656143873	0.000256038
FLNC	1.112773551	3.62254181	0.000291722
GHR	1.161660284	3.622453812	0.000291822
PRKG1	1.195566092	3.621918942	0.000292426
TNFAIP8L3	1.171369232	3.611019416	0.000304996
ADRA1D	1.17070966	3.553028545	0.000380823
AIFM3	0.853346363	−3.547409751	0.000389039
ABCA4	1.103837888	3.521258054	0.000429504
SRPX	1.146802308	3.515246382	0.000439346
CACNA2D1	1.156053229	3.499144698	0.000466753

Figure 4

Survival analysis and Cox regression. (A) A heatmap of 12 significant mRNA expression profile for prediction of overall survival by multivariate cox regression. (B) Twelve mRNA expression profiles for prediction of overall survival in bladder cancer by multivariate Cox regression (C) Kaplan-Meier plot shows statistical significance between high-risk and low-risk groups by risk scoring model (P<0.05). (D) The ROC curve for the risk scoring model (AUC=0.735). ROC – receiver operating characteristic; AUC – area under the curve.

PPI network

In order to identify the gene interactions of the differentially expressed mRNAs in the ceRNA network, a PPI network was constructed with these 52 aberrantly expressed genes. In total, 44 nodes and 144 edges were constituted in the single PPI network by STRING (Figure 5A). JUN was identified to have the largest number of edges interacting with 28 other genes. Nodes which connect with more than 4 genes are listed in the barplot (Figure 5B). Based on the predicted combined score for each pair of nodes from STRING, the top 10 correlating mRNAs are presented in Table 8. We then proceeded by validating the 10 correlations with Pearson’s correlation and as shown in Figure 6, we found a strong correlation between CFL2 and TPM2 (r=0.911) and between CCNE2 and E2F7 (r=0.734), respectively (Figure 6). The results indicated that there was an interaction between CFL2 and TPM2 and between CCNE2 and E2F7 respectively, in bladder cancer.

Figure 5

The protein–protein interaction (PPI) network. (A) A PPI network of differentially expressed genes in the ceRNA network. (B) The order of different edges of each gene in the PPI network.

Table 8

Top 10 correlating mRNAs presented in the protein–protein interaction network.

mRNA1	mRNA2	Combined score
JUN	MYB	0.973
MYB	CCNB1	0.964
CCNE1	CCNB1	0.949
JUN	SELE	0.949
E2F7	CCNE2	0.925
FGF9	FGF2	0.919
CCNE2	CCNE1	0.91
FGF2	THBS1	0.861
JUN	FGF2	0.811
TPM2	CFL2	0.793

Figure 6

The correlation plots of the top 10 predicted combined model.

Clinicopathological features related gene expression in ceRNA network

We further analyzed the relationship between the differentially expressed RNAs in the ceRNA network in relation to clinicopathological features. As a result, we found 5 mRNAs, 21 lnRNAs, and 6 miRNAs to be differentially expressed among high-grade and low-grade bladder cancer (Figure 7, Table 9). In addition, of these RNAs, 2 out of 5 mRNAs (HOXB5, MKX) and 10 out of 21 lncRNAs (AC112721.1, ADAMTS9-AS1, AP004609.1, LINC00460, AC110491.1, C20orf166-AS1, LINC00163, ARAP1-AS1, HCG22, LINC00473) were found to be significantly differentially expressed across stage I, stage II, stage III, and stage IV (Figure 8). The rest of the RNAs were not associated with tumor stage in bladder cancer (P>0.05) (Supplementary Figure 1). In relation to prognosis, 1 mRNA (HOXB5) and 6 lncRNAs (ADAMTS9-AS1, AC112721.1, LINC00460, AC110491.1, LINC00163, HCG22) were observed to be associated with overall survival (Figure 9).

Figure 7

A heatmap of differentially expressed genes in the ceRNA network associated with high-grade bladder cancer. (A) Differentially expressed mRNAs in ceRNA network associated with high-grade bladder cancer. (B) Differentially expressed lncRNAs in ceRNA network associated with high-grade bladder cancer. (C) Differentially expressed miRNAs in ceRNA network associated with high-grade bladder cancer.

Table 9

Differentially expressed RNA profiles correlating to high-grade bladder cancer from the ceRNA network.

RNAs	Log fold change	FDR
HOXB5(mRNA)	−2.944596082	2.39E-18
MKX(mRNA)	3.512064708	3.98E-08
SALL3(mRNA)	−3.009599773	6.86E-07
ELAVL2(mRNA)	3.104383407	0.000197897
SELE(mRNA)	2.275611899	0.000649315
LINC00336(lncRNA)	−2.274109048	1.50E-09
ARAP1-AS1(lncRNA)	−2.922663469	3.31E-09
LINC00460(lncRNA)	6.246084381	1.13E-06
AC112721.1(lncRNA)	6.126188656	1.90E-06
ADAMTS9-AS1(lncRNA)	4.959254878	8.37E-06
C2orf48(lncRNA)	2.14039587	2.37E-05
MIR137HG(lncRNA)	6.327534674	5.76E-05
LINC00473(lncRNA)	5.415263905	6.33E-05
LINC00163(lncRNA)	4.423919282	6.33E-05
AC110491.1(lncRNA)	4.746368747	0.000165643
LINC00402(lncRNA)	3.533793424	0.000229383
MYO16-AS1(lncRNA)	5.102890553	0.000231527
LINC00520(lncRNA)	4.493602761	0.000289874
HCG22(lncRNA)	6.362121452	0.000499845
TLR8-AS1(lncRNA)	4.019357381	0.000964633
C20orf166-AS1(lncRNA)	3.194349357	0.001177387
LINC00518(lncRNA)	2.965146588	0.004120812
AP004609.1(lncRNA)	2.067991973	0.004431844
SACS-AS1(lncRNA)	2.266344736	0.026097183
LINC00330(lncRNA)	3.348989332	0.026097183
AC127496.3(lncRNA)	2.19115354	0.034777432
hsa-mir-143(miRNA)	3.330527062	0.00099772
hsa-mir-215(miRNA)	2.871943848	0.002121278
hsa-mir-217(miRNA)	3.321184311	0.002354727
hsa-mir-372(miRNA)	4.593394591	0.00371048
hsa-mir-373(miRNA)	4.640595006	0.0118349
hsa-mir-383(miRNA)	2.381514994	0.024564256

Figure 8

Histograms of differentially expressed genes in the ceRNA network related to high-grade bladder cancer and disease stages. (A) Two significant mRNA expressions across different tumor stages. (B) Ten significant lncRNAs expressions across different tumor stages.

Figure 9

The Kaplan-Meier plots of differentially expressed RNAs in the ceRNA network which are also related to high-grade, stages, and overall survival.

Discussion

Molecular and genetic alterations in bladder cancer have been extensively studied in the past decades. Clinicopathological features, including tumor grade and stage have been shown to be linked to biological and genetic mechanisms [22-24]. However, clinical challenges still exist for bladder cancer making it crucial to focus on new prediction methods based on the molecular targets.

In this study, we have conducted an integrated characterization of RNA expression profiles and performed ceRNA network analysis in bladder cancer based on data from the TCGA database. As a result, 1819 mRNAs, 659 lncRNAs, and 160 miRNAs were found to be differentially expressed in bladder cancer, of which 52 mRNAs, 58 lncRNAs, and 22 miRNAs were predicted to be involved in the ceRNA network. Furthermore, high-grade bladder cancer related RNA profiles were identified consisting of 5 mRNAs, 21 lnRNAs, and 6 mRNAs, of which 2 mRNAs and 10 lncRNAs were found to be associated with stage as well.

Based on the current findings of 1819 differentially expressed RNAs, our functional enrichment analysis indicated that that these genes were mainly enriched in “extracellular region”. Similarly, GO analysis of prostate cancer conducted by Jiang et al. showed that most of the differentially expressed mRNAs in the urinary system were aggregated in the extracellular region [25]. In addition, we predicted that 61 genes were enriched in the neuroactive ligand-receptor interaction, including PDGFR-α. Activation of PDGFR-α have previously been implicated in bladder tumor progression by a ras- and Src-independent activation of MEK/ERK pathway [26]. Interestingly, systematic lupus erythematosus and alcoholism were identified as important components in the occurrence and development of bladder cancer as previous reported [27,28].

Based on the results from the Cox regression, we successfully established an overall survival prediction model based on 12 mRNA expression profiles. The discriminative value of these genes combined was an AUC at 0.735. Of these 12 genes, the tumor necrosis factor-α induced protein 8 TNFAIP8, has been shown to be strongly associated to cancer progression [29]. However, except for TNFAIP8 in the prediction model, the other 11 genes have not yet been elucidated in bladder cancer. Hence, the prognostic value of the prediction model should be further demonstrated in future studies.

Through our analysis of the PPI network, we identified that JUN, known as c-Jun, plays a vital role in the differentially expressed mRNAs in our ceRNA network. In the analysis, we found that JUN has the most edges and interacts with 28 other genes. JUN is an AP-1 transcription factor subunit, which exhibits proto-oncogenic functions, and several studies have previously validated that the tumorigenesis of bladder cancer has frequently occurred through c-jun relevant signaling pathways [30-34]. Moreover, Chen et al. reported that TMP2 expression was significantly decreased in bladder cancer [35], which was in accordance with our findings that demonstrated the strong positive correlation between TPM2 and CFL2 expressions, where CFL2 was also found to be downregulated in bladder cancer with the lowest FDR. CFL2 has also been reported to act as a tumor suppressor gene in nasopharyngeal carcinoma, pancreas cancer, and gastric cancer [36-38]. Although the functions of CFL2 in bladder cancer has not yet been studied comprehensively, our results indicated that it is a promising novel and potentially meaningful biomarker for exploring new mechanisms in bladder cancer.

As we know, clinicopathological features such as grading and staging are systematic and important prognostic factors in bladder cancer, but their relationship to gene expression profiles remains unclear. In our study, we constructed a bladder cancer specific ceRNA network and analyzed the differentially expressed mRNAs and endogenous RNAs in relation to tumor grades and stages. As a result, we identified that 2 mRNAs and 10 lncRNAs were significantly associated to high-grade tumors and different stages. However, in survival analyses, only HOXB5 and 6 lncRNAs (ADAMTS9-AS1, AC112721.1, LINC00460, AC110491.1, LINC00163, HCG22) were related to overall survival. Although HOXB5 is upregulated in bladder cancer correlating with high-grade and high-stage tumors and has a miRNA-7-binding site of HOXB5 3′-UTR SNP [39], the mechanisms of HOXB5 in bladder cancer remains to be further explored in future studies. Of the survival related lncRNAs in our study, Ye et al. previously demonstrated that LINC00460 could facilitate tumor progression by acting as a sponge to miR-302c-5p and thereby regulating FOXA1 signaling pathway in human lung adenocarcinoma; LINC00460 has also been implicated in promoting malignant biological behaviors in gastric cancer by regulating KDM2A expression through the targeting of miR-342-3p [40,41]. However, none of the 6 survival-related lncRNAs have been reported in bladder cancer. In summary, we identified that HOXB5, ADAMTS9-AS1, AC112721.1, LINC00460, AC110491.1, LINC00163, and HCG22 were related to overall survival in bladder cancer. These genes could be promising future biomarkers for diagnosis, prognostication, and also neo-therapeutic targets in bladder cancer.

Conclusions

We have identified a large number of differentially expressed mRNAs, lncRNAs, and miRNAs in bladder cancer which were strongly associated with oncogenesis and prognosis. Many of these RNAs have not been reported in the current literature as related to bladder cancer and represent novel and promising future targets. However, larger cohorts of patients and future mechanistic studies are necessary to validate our results, and investigate their functional roles in bladder cancer.

Histograms of differentially expressed genes in the ceRNA network unrelated to high-grade and stages. (A) Three mRNAs expressions were unrelated to tumor stages. (B) Eleven lncRNAs expressions were unrelated to tumor stages. (C) Six mRNAs expressions were unrelated to tumor stages.