Identification and analysis of long non-coding RNA related miRNA sponge regulatory network in bladder urothelial carcinoma

Background The aim of this study was to investigate the regulatory network of lncRNAs as competing endogenous RNAs (ceRNA) in bladder urothelial carcinoma (BUC) based on gene expression data derived from The Cancer Genome Atlas (TCGA). Materials and methods RNA sequence profiles and clinical information from 414 BUC tissues and 19 non-tumor adjacent tissues were downloaded from TCGA. Differentially expressed RNAs derived from BUC and non-tumor adjacent samples were identified using the R package “edgeR”. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using the “clusterProfiler” package. Gene ontology and protein–protein interaction (PPI) networks were analyzed for the differentially expressed mRNAs using the “STRING” database. The network for the dysregulated lncRNA associated ceRNAs was then constructed for BUC using miRcode, miRTarBase, miRDB, and TargetScan. Cox regression analysis was performed to identify independent prognostic RNAs associated with BUC overall survival (OS). Survival analysis for the independent prognostic RNAs within the ceRNA network was calculated using Kaplan–Meier curves. Results Based on our analysis, a total of 666, 1819 and 157 differentially expressed lncRNAs, mRNAs and miRNAs were identified respectively. The ceRNA network was then constructed and contained 59 lncRNAs, 23 DEmiRNAs, and 52 DEmRNAs. In total, 5 lncRNAs (HCG22, ADAMTS9-AS1, ADAMTS9-AS2, AC078778.1, and AC112721.1), 2 miRNAs (hsa-mir-145 and hsa-mir-141) and 6 mRNAs (ZEB1, TMEM100, MAP1B, DUSP2, JUN, and AIFM3) were found to be related to OS. Two lncRNAs (ADAMTS9-AS1 and ADAMTS9-AS2) and 4 mRNA (DUSP2, JUN, MAP1B, and TMEM100) were validated using GEPIA. Thirty key hub genes were identified using the ranking method of degree. KEGG analysis demonstrated that the majority of the DEmRNAs were involved in pathways associated with cancer. Conclusion Our findings provide an understanding of the important role of lncRNA–related ceRNAs in BUC. Additional experimental and clinical validations are required to support our findings.

common type is BUC and accounts for more than 90% of all BCs [2]. Unfortunately, the high recurrence rate is a characteristic of BUC [3,4]. Currently, the principal treatment strategy for BUC consists of surgery and adjuvant combination chemotherapy. However, chemotherapy resistance reduces the sensitivity of BUC to chemotherapeutic drugs and frequently results in treatment failure resulting in BUC clinical management being a major challenge [5]. In addition, only a limited number of biomarkers are available for diagnosing BUC compared to other cancers. Hence, identifying sensitive and specific BUC biomarkers, as well as therapeutic targets for BUC are critically needed.
Long non-coding RNAs (lncRNAs) are a subtype of ncRNAs with transcript lengths over 200 nucleotides and have recently attracted increased attention [6]. lncRNAs were initially regarded as transcriptional noise without the capacity to encode proteins [7]. However, growing evidence has demonstrated that lncRNAs may play crucial biological roles in a variety of biological processes that are associated with carcinogenesis and cancer metastasis [8]. With regards to bladder cancer, several studies have suggested that lncRNAs may function as oncogenes or tumor suppressors and may affect overall patient survival and mortality [9,10]. To date, only a few lncRNAs have been verified experimentally, but their roles in regulating gene expression remains to be deciphered.
Considerable efforts have been made to demonstrate how lncRNAs exert their diverse biological functions in human malignant tumors. Rapid progression has been made to elucidate the role of lncRNAs in miRNA function. miRNAs are endogenous single-stranded RNA with lengths between 20 and 25 nucleotides that do not encode proteins. They repress gene expression by complement binding to their target mRNA sequences (i.e. microRNA response element, MRE) [11]. In 2011, Salmena et al. [12] proposed the competing endogenous RNA (ceRNA) hypothesis, which states that mRNAs, transcribed pseudogenes, and lncRNAs could act as natural miRNA "sponges" and inhibit miRNA function by competing with the binding of one or more MREs in complex and comprehensive regulatory networks, leading to pathogenic conditions. The ceRNA regulation theory has been proven to be involved in bladder cancer initiation and progression in several studies [13,14]. Similarly, lncRNAs acting as ceRNAs has also been reported in other cancers [15][16][17]. Recently, Kouhsar et al. [18] constructed a ceRNA network related to the staging of Non-Muscle Invasive Bladder Cancer (Ta and T1) derived from public data sources. They identified several biomarkers associated with tumor stage. We hypothesized that lncRNAs may function as ceRNAs during BUC initiation and progression. Understanding how lncRNAs function as ceRNAs will be important in deciphering BUC carcinogenesis.
In the present study, we aimed to decipher the regulatory ceRNA network of lncRNAs-miRNAs-mRNAs in BUC by analyzing gene expression data. This was performed using bioinformatics prediction and correlation analyses. In addition, using clinical trials and survival analyses, we identified potential prognostic genes.

Study cohort
RNA sequence data from 406 BUC patients were retrieved from The Cancer Genome Atlas (TCGA) database (https ://cance rgeno me.nih.gov/) in 2018. LncRNA, miRNA, and mRNAseq data were obtained using the Data Transfer Tool (provided by GDC Apps). Patient clinical information was also downloaded using the Data Transfer Tool. Sequencing data derived from the Illumina HiSeq RNAseq and Illumina HiSeq miRNAseq platforms were publicly available. This study met the publication guidelines stated by TCGA (https ://cance rgeno me.nih. Gov/publi catio ns/publi catio nguid eline s). All data used in the study were obtained from TCGA, and hence ethics approval and informed consent were not required.

Differential expression analysis
The mRNAseq and lncRNAseq data derived from 414 BUC tissue samples and 19 non-tumor adjacent tissue samples were downloaded from TCGA. The BUC miR-NAseq data were derived from 418 BUC tissue samples and 19 non-tumor adjacent tissue samples. For tumor and non-tumor group comparison, differentially expressed mRNAs (DEmRNAs), miRNAs (DEmiRNAs) and lncRNAs (DElncRNAs) were identified using the "edgeR" package (http://bioco nduct or.org/packa ges/relea se/bioc/html/edgeR .html) with a cut-off criteria of |log2 (fold change [FC])| > 2.0 and adjusted P value < 0.01 [19]. Differentially expressed lncRNAs (DElncRNAs) were defined and annotated using the Encyclopedia of DNA Elements (ENCODE), which included 15,877 human lncRNAs. All P-values used the False discovery rate (FDR) to correct for statistical significance of multiple testing (Benjamini-Hochberg method) [20]. FDR significance level was set at 0.05.

KEGG enrichment analysis of DEmRNAs
KEGG enrichment analysis was performed using the "clusterProfiler" package in R software based on the retrieved DEmRNAs and visualized using the Cytoscape v 3.5.1 software.

Protein-protein interaction (PPI) network
To understand the interactions of the DEmRNAs, we constructed a PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING, http:// strin g.embl.de/). Combined scores greater than 0.4 were considered statistically significant. The PPI network was visualized using the Cytoscape v 3.5.1 software. Subsequently, the top 30 mRNAs were identified using the ranking method of degree. In addition, gene ontology (GO) enrichment analysis was performed using STRING to functionally annotate the DEmRNAs in BUC.

Independent prognostic factors for overall survival (OS)
Independent prognostic differentially expressed RNAs for OS were identified by univariate and multivariate cox regression analyses. The risk score (RS) was estimated using the following formula: where Exp(i) denotes the expression value for RNA; n represents the number of RNA modules, and R(i) denotes the estimated regression coefficient of RNA. If the RS for a given sample was greater than the mean RS of all samples, the latter was regarded as a high-risk sample, otherwise, it was regarded as a low-risk sample. Kaplan-Meier method was used to evaluate the survival curves in the high-and low-risk groups. Additionally, the sensitivity and specificity were accessed using the receiver operating

characteristic (ROC) curves and area under the ROC curves (AUC values).
Survival curves were plotted using the "survival" package in R for the independent prognostic RNAs that were identified. Long-rank test was used to evaluate statistical significance and P < 0.05 was considered statistically significant.
Validating the prognostic value of lncRNAs and mRNAs in the ceRNA network using GEPIA Independent prognostic lncRNAs and mRNAs were validated using Gene Expression Profiling Interactive Analysis (GEPIA,http://gepia .cance rpku.cn/index .htm), which was based on RNA sequencing data from 9736 tumors and 8587 normal samples in the Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) dataset project [21]. In addition, the correlation between lncRNAs and mRNAs was confirmed using the Pearson correlation statistic. The correlation of two RNAs was considered significant when R was greater than 0.4 and the P-value was less than 0.05.

Pathway enrichment analysis of DEmRNAs
To investigate the mechanisms associated with BUC tumorigenesis, 1819 DEmRNAs were used for KEGG enrichment analysis. The threshold was set at P < 0.01. The top 15 significantly enriched pathways are presented in Table 1. The results showed that the majority of DEmRNAs were enriched for "neuroactive ligandreceptor interaction", "viral carcinogenesis", "protein digestion and absorption", "ECM-receptor interaction", and "cAMP signaling pathway". In addition, a network of pathways and DEmRNAs was constructed based on the above KEGG analysis and visualized using Cytoscape v 3.5.1 (Additional file 2: Figure S2). Interestingly, several DEmRNAs were related to more than one pathway, such as SLC8A1, CAMK2B, and CACNA2D1. CAMK2B was enriched for Circadian entrainment, calcium signaling

Functional enrichment analysis of DEmRNAs
Functional enrichment of the top 600 DEmRNAs were analyzed using STRING (300 significantly up-regulated and 300 significantly down-regulated DEmRNAs based on log fold change (log FC)). The GO function for the various genes were divided into biological processes (BP), cellular function (CF) and molecular component (MC). The GO results are presented in Fig. 1. The results from the GO functional enrichment analysis suggested that DEmRNAs were significantly enriched in multicellular organisms (ontology: BP), binding (ontology: CF) and extracellular region (MC). These results demonstrated that the enriched DEmRNAs were associated with BUC proliferation and migration.

Protein-protein interaction (PPI) network
The interrelationship between the top 600 significant DEmRNAs (300 significantly up-regulated and 300 significantly down-regulated DEmRNAs based on log fold change (log FC)) were retrieved from the STRING database to construct the PPI network (Fig. 2). The PPI network consisted of 418 nodes and 1937 edges. The nodes denoted DEmRNAs, while the edges denoted interactions among the DEmRNAs. Additionally, the top 30 mRNAs in the PPI network were analyzed using the ranking method of degree (Additional file 3: Figure S3). KEGG pathway enrichment analysis for the 30 mRNAs identified with a high degree was performed using the "clusterProfiler" package in the R software with a P < 0.05 as the cut-off criteria. The results demonstrated that these mRNAs were enriched for pathways related to "Systemic lupus erythematosus", "Alcoholism", "viral carcinogenesis" and "transcriptional dysregulation in cancer" (Additional file 4: Figure S4A). Furthermore, we identified 21 (HIST2H2BF, HIST1H2BO, HIST1H2BH, HIST1H2BE, HIST1H2BF, HIST1H4E, HIST1H2BM, HIST1H4D, HIST1H2BI, HIST1H4C, HIST1H4B, HIST1H4A, HIST1H2BB, HIST1H4F, HIST1H4L, IL6, HIST1H3G, HIST1H3B, HIST1H3C, HIST1H3F and HIST1H3I) of these 30 DEmRNAs that were enriched for cancer-related pathways: "viral carcinogenesis" or "transcriptional dysregulation in cancer". These 21 DEmRNAs were classified as hub genes. In addition, a network linking the pathways and mRNAs were constructed and visualized using Cytoscape v 3.5.1 (Additional file 4: Figure  S4B).

Construction of the ceRNA network for BUC
To better understand the role of the identified differentially expressed RNAs in BUC, a dysregulated ceRNA network based on DElncRNA-DEmiRNA-DEmRNA interactions was constructed using Cytoscape v 3.5.1 (Fig. 3). The regulatory relationship between DEmiR-NAs and DElncRNAs pairs were retrieved from miRcode (http://www.mirco de.org/). All the DEmiRNAs and DEl-ncRNAs involved in the ceRNA network meet the cutoff criteria (|log2FC| > 2.0 and adjusted P-value < 0.01) mentioned above. We found 259 lncRNA-miRNA interaction pairs containing 59 lncRNAs and 23 DEmiRNAs from the miRcode database. Subsequently, we searched for mRNAs that were targeted by the 23 DEmiRNAs using miRTarBase (http://mirta rbase .mbc.nctu.edu.tw/),  and MEST) were closely related to OS in BUC patients (P < 0.05) ( Table 2 and heat maps in Additional file 5: Figure S5A-C). Based on risk scores (RS) of the independent prognostic RNAs identified above, patients included in the present study were assigned into high-and lowrisk groups. The low-risk group had a significantly better prognosis compared to the high-risk group (Fig. 4a for lncRNA, B for miRNA and C for mRNA). The results showed that the 5-year OS for the low-risk group was 56.4%, 55.0%, and 55.8% respectively, while it was 26.5%, 29.1% and 27.8% for the high-risk group, respectively. To determine the prognostic power of the 7 lncRNAs, 3 DEmiRNAs, and 7 DEmRNAs identified above, timedependent receiver operating characteristic (ROC) curve analysis was performed and the area under the curve (AUC values) was estimated. AUC values for the independent prognostic DElncRNAs, DEmiRNAs, and DEm-RNAs were 0.707, 0.624 and 0.681, respectively (Fig. 5, a  Kaplan-Meier curve analysis was performed to determine the OS for the independent prognostic RNAs. One patient was lost during follow-up and was excluded from the survival analysis. Five DElncRNAs were significantly related to OS, of which, four DElncR-NAs (HCG22, ADAMTS9-AS1, ADAMTS9-AS2, and AC112721.1) were negatively related to OS (Fig. 6b-e), while AC078778.1 was positively related to OS (logrank P < 0.05) (Fig. 6a). In addition, Kaplan-Meier curve analysis for the three DEmiRNAs and seven DEmRNAs showed that two DEmiRNAs (hsa-mir-141 and hsamir-145) and 6 DEmRNAs (ZEB1, TMEM100, MAP1B, DUSP2, JUN, and AIFM3) were significantly related to OS (log-rank P < 0.05) (Figs. 7a, b, 8a-f ).

Prognostic value of lncRNAs and mRNAs in the ceRNA network analyzed using GEPIA
GEPIA was used to validate the expression levels and prognostic value of the 5 independent lncRNAs. The expression levels of HCG22, ADAMTS9-AS1 and ADAMTS9-AS2 were negative or nearly negative, while AC112721.1 expression levels were positive in BUC tissues. These findings were concordant with previous results using TCGA analysis. However, only ADAMTS9-AS1 and ADAMTS9-AS2 expression levels were downregulated with statistical significance in BUC. The expression levels of these 4 lncRNAs are shown in Fig. 9. However, AC078778.1 was not found in GEPIA, while the overall survival of the 2 lncRNAs, ADAMTS9-AS1 and ADAMTS9-AS2, were estimated using GEPIA. The results are shown in Fig. 10 and were similar to our previous results. Lower expression of ADAMTS9-AS1 and ADAMTS9-AS2 were associated with a good prognosis in BUC patients. Similar to the lncRNAs, 6 independent prognostic mRNAs (ZEB1, TMEM100, MAP1B, DUSP2, JUN, and AIFM3) in the ceRNA network were also analyzed using GEPIA. The six mRNA expression levels are shown in Fig. 11a-f. ZEB1, TMEM100, MAP1B, DUSP2 and JUN were down-regulated with statistical significance, while AIFM3 was up-regulated without statistical significance in BUC. The correlation with

Discussion
In this study, five DElncRNA factors (HCG22, ADAMTS9-AS1, ADAMTS9-AS2, AC078778.1, and Fig. 12 Survival analysis of the five-independent prognostic lncRNAs (a for DUSP2, b for JUN, c for MAP1B, d for TMEM100 and e for ZEB1) in BUC patients analyzed using GEPIA. The median expression levels of the five-independent prognostic lncRNAs were set as the thresholds for demarcating high and low-expression cohorts. Log Rank was used and P < 0.05 was considered statistically significant. TPM transcripts per million. Horizontal axis: overall survival time: months, Vertical axis: survival Fig. 13 Correlation of the independent prognostic lncRNAs and mRNAs (a for ADAMTS9-AS2 and TMEM100, b for ADAMTS9-AS2 and ZEB1, c for ADAMTS9-AS1 and ZEB1). R values greater than 0.4 and P-values < 0.05 were considered statistically significant AC112721.1) in the ceRNAs network were identified as independent prognostic factors for OS in BUC patients. HCG22 expression levels have been reported to be down-regulated in oral cancer and its low expression was associated with poor survival in a recent study based on TCGA data analysis [22]. Lu et al. [23] investigated HCG22 expression levels in 20 oral cavity and oropharyngeal squamous cell carcinoma (OSCC) samples and 10 control samples by qRT-PCR. They demonstrated that HCG22 was downregulated in OSCC tissues compared to controls, while no association was observed between HCG22 expression levels and overall survival. ADAMTS9-AS1 (ADAMTS9 antisense transcript), is a novel lncRNA without any functional annotation but could interact with two RNA-binding proteins, DCGR8 and FUS [24]. lncRNAs play an important role in cancer mainly via their associations with RNA-binding g proteins, which include HOTTIP, MaLAT1, H19, and HOTAIR. They participate in several biological pathways involved in cell differentiation and proliferation, apoptosis and tumorigenesis by interacting with RNAbinding proteins in hepatocellular carcinoma [25]. We hypothesized that ADAMTS9-AS1 may play a role in the development of cancer. Wang et al. [26] found that ADAMTS9-1 and ADAMTS9-2 expression levels were decreased in malignant epithelial ovarian cancer tissues compared to normal ovary tissues and benign ovarian cysts using lncRNA and mRNA microarray analysis. These results were confirmed using 8 normal ovarian, 17 benign ovarian cysts and 15 malignant epithelial ovarian cancer samples by qPCR assays. Low ADAMTS9-AS2 levels were found to be a significant independent predictor of poor survival in glioma patients [27]. Liu et al. [28] suggested that lncRNA ADAMTS9-AS2 could suppress cancer progression by inhibiting miR-223-3p and activating TGFBR3. Additionally, increased ADAMTS9-AS2 levels could reduce lung cancer tumor size and improve OS. However, to date, no studies have been performed to determine the role of AC078778.1 and AC112721.1 in cancer. miRNAs are involved in multiple roles during carcinogenesis. In this study, we found two independent prognostic DEmiRNA factors (hsa-mir-141 and hsa-mir-145) that were involved in the ceRNA network. miR-141 was found to be up-regulated in malignant bladder tissue samples compared to healthy tissues and was a favorable prognostic biomarker [29]. microRNA-141 (hsa-mir-141) has been shown to exert a regulatory role during epithelial to mesenchymal transition process, and its expression levels have been associated with tumorigenicity and invasiveness in several human cancers. hsa-mir-141 has been shown to be associated with the development of certain epithelial cancer cell types, including prostate [30], colorectal [31] and breast cancer [32]. Huang et al. [33] demonstrated that miR-141 could inhibit gastric cancer cell proliferation and tumor growth, while low miR-141 levels were associated with poor prognosis. Another study on gastric cancer suggested that miR-141 could play an important anti-tumor role by interacting with MEG3 and targeting E2F3 during gastric cancer pathogenesis and may be a therapeutic target. miR-145 has been frequently observed to be down-regulated in cancers and restoration of miR-145 levels suppressed cancer cell invasion by reversing the EMT phenotype [34]. Tan et al. [35] demonstrated that TUG1 promoted bladder cancer cell metastasis and radio-resistance by negatively regulating miR-145 expression.
In the present study, 6 prognostic DEmRNA factors (ZEB1, TMEM100, MAP1B, DUSP2, JUN, and AIFM3) were involved in the ceRNAs network and functioned as independent prognostic factors for OS in BUC patients. Several studies have demonstrated that ZEB1 was significantly overexpressed in bladder cancer tissues compared to normal healthy adjacent tissues [36]. Li et al. [37] reported that ZEB1 was significantly overexpressed in bladder cancers compared to normal tissues, and played a crucial role during VM formation, and was closely associated with invasion, metastasis and poor prognosis of malignant tumors [38,39]. Transmembrane protein 100 (TMEM100), located at 17q32, was first identified as a transcript in the mouse genome. Han et al. [32] found that TMEM100 could function as a tumor suppressor by inhibiting the growth and metastasis of non-small-cell lung cancer via the inhibition of the TNF pathway. Low TMEM100 expression levels were associated with poor prognosis. Similar results have been reported for hepatocellular carcinoma [40]. MAP1B, which encodes for the microtubule-associated protein 1B (MAP1B), is one of the main cytoskeletal proteins. Several studies have demonstrated that MAP1B plays an important role in a number of cellular processes, including synaptic transmission, autophagy, and cancer [41][42][43]. DUSP2 is a member of the class 1 DUSP family of proteins and is localized in the nucleus. DUSP2 levels are significantly decreased in bladder cancer and low expression of DUSP2 is correlated with poor prognosis [44]. JUN (C-Jun, AP-1 transcription factor subunit) is specifically phosphorylated by JNK and plays a central role in the AP-1 complex. It is involved in cellular DNA damage response by regulating the expression of several genes [45,46]. c-Jun is a protooncogene and is involved in transformation and tumor development [47,48]. AIFM3 (apoptosis-inducing factor mitochondria associated 3) is a gene with homology to apoptosis-inducing factor (AIF). AIF induces apoptosis in a caspase-dependent manner. AIFM3 was observed to be highly expressed in breast cancer tissues and associated with shorter overall survival and disease-free survival [49].
BUC specific DEmRNA pathways were assessed using KEGG pathway analysis. Our results demonstrated that the majority of DEmRNAs were frequently enriched for cancer-related pathways. Of these, the neuroactive ligand-receptor interaction signaling pathway has been reported to be associated with the progression of renal cell carcinoma in bioinformatics studies [50]. In addition, previous studies have suggested that cAMP-related signaling could control apoptosis induction and cell growth [51,52], while another study demonstrated that cAMP was an inhibitor of cell cycle progression and apoptosis in gastric cancer cells [53]. Zhang et al. [54] showed that the ECM-receptor interaction pathway played a significant role in tumor progression and metastasis. We found several overlapping DEmRNAs that were involved in multiple pathways, such as SLC8A1, CAMK2B, and CACNA2D1, and have been demonstrated to play important roles in cancer pathogenesis. Muñoz reported that lower levels of SLC8A1, which was at least partly mediated by miR-223, was associated with reduced calcium and apoptosis levels in penile carcinoma [55]. Based on an integrative meta-analysis, CAMK2B was found to be associated with the development of cancer cachexia [56]. Recently, Feng et al. [57] assayed 15 colorectal cancer tissues and 10 paracancerous tissues using microarrays and found that CAMK2B was involved in the progression of Fusobacterium nucleatum-induced colorectal cancer. Another study also reported that CAMK2B played an important role in glioblastoma multiforme using bioinformatics analysis on publicly available datasets [58]. High expression levels of CACNA2D1 in epithelial ovarian cancers were significantly correlated with histological subtypes, advanced FIGO stages and tumor differentiation [59]. In addition, the majority of PLB4 hotspot mutations are gain-of function mutations that have been demonstrated to be involved in uveal melanoma tumorigenesis by activating the same signaling pathway [60].
In terms of GO functional enrichment analysis, we found that DEmRNAs were significantly enriched for the regulation of multiple processes, such as binding, protein binding, receptor binding, ion channel activity, and endopeptidase activity. Interestingly, several studies have reported that these functions overlapped in different cancers. Several studies have demonstrated that knocking down ezrin and P65 expression induces tumor metastasis in different cancers [61,62]. Furthermore, Tang et al. [63] demonstrated that ezrin and P65 were physically associated with one another. We hypothesize that the interaction between ezrin and P65 is associated with the activation of the NF-κB pathway leading to breast cancer metastasis. A previous study suggested that Estrogen receptor β could increase the levels of miR-92a by binding to the estrogen-response-element (ERE) leading to a decrease in DAB2IP tumor suppressor expression to ultimately promote bladder cancer growth and invasion [64]. Numerous studies have demonstrated that ion channels play an important role in tumorigenesis and progression, such as inducing neo-angiogenesis [65], apoptosis resistance [66], proliferative potential [67] as well as cell migration and invasiveness [68,69]. In terms of endopeptidase activity, Zhu et al. [70]. demonstrated that asparaginyl endopeptidase (AEP) was highly expressed in tissues and ascites of patients with epithelial ovarian cancer and promoted tumor growth and progression both in vivo and in vitro.
We then constructed a PPI network to identify hub DEmRNAs. Proteins that corresponded to genes were used to build the PPI network, and the top 30 DEmRNAs with a high degree were selected. KEGG pathway analysis of these 30 DEmRNAs was performed using the "cluster-Profiler" package. The results showed that 21 of these 30 DEmRNAs were enriched for "viral carcinogenesis" and "transcriptional dysregulation in cancer". The 21 identified DEmRNAs may play an important role in cancer. For example, IL6, HIST1H3C, and HIST1H3G, which were classified with high degrees are present in pathways associated with transcriptional dysregulation in cancer. In addition, these three genes have been previously reported to be closely associated with tumorigenesis and development [71][72][73].
In terms of correlation analysis, a positive correlation between ZEB1 and ADAMTS9-AS1-AS2 and TMEM100 was observed using the Pearson correlation statistic based on the expression levels of these genes. The correlation between these four genes may play an important role in the initiation and progression of BUC. Only a few studies on their interactions have been published in public databases, such as Pubmed and Embase.
Several limitations of the present study should be stated. First, the number of normal bladder tissues (19 samples) was limited and may have compromised the reliability of our results. Second, BUC patient information from TCGA was not validated using experimental procedures. Third, we only investigated the ceRNAs network associated with lncRNAs, miRNAs, and mRNAs, and did not include other regulatory models. Finally, several novel lncRNAs with significant clinical value needs to be investigated to determine their functional role and underlying mechanism during BUC carcinogenesis. We verified the independent prognostic lncR-NAs and mRNAs using GEPIA and this supported our findings. Our results provide a better understanding of lncRNA-related ceRNAs and its important role in BUC.