Molecular network-based identification of competing endogenous RNAs in bladder cancer

Background Circular RNAs (circRNAs) have been shown to interact with microRNAs (miRNA) as competitive endogenous RNAs (ceRNAs) to regulate target gene expression and participate in tumorigenesis. However, the role of circRNA-mediated ceRNAs in bladder cancer (BC) remains unknown. Accordingly, the aim of this study was to elucidate the regulatory mechanisms in BC based on construction of the ceRNA network. Methods The RNA expression profiles were obtained from public datasets in the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) database, and were used to establish a circRNA-miRNA-mRNA network. The interactions among proteins were analyzed using the STRING database and hubgenes were extracted using the cytoHubba application. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of differentially expressed mRNAs in BC and normal tissue samples were performed to determine the functions of the intersecting mRNAs. Results A total of 27 circRNAs, 76 miRNAs, and 4744 mRNAs were found to be differentially expressed between BC and normal tissues. The circRNA-miRNA-mRNA ceRNA network was established based on 21 circRNAs, 14 miRNAs, and 150 mRNAs differentially expressed in BC. We also established a protein-protein interaction network and identified 10 hubgenes, which were used to construct circRNA-miRNA-hubgene regulatory modules. The most enriched biological process GO term was strand displacement (P<0.05), and the homologous recombination and Fanconi anemia pathways were significantly enriched (P<0.05) for the differentially expressed genes in BC. Conclusions We screened several dysregulated circRNAs and established a circRNA-associated ceRNA network by bioinformatics analysis. The identified ceRNAs are likely critical in the pathogenesis of BC and may serve as future therapeutic biomarkers.

Introduction were obtained from TCGA database. No ethical approval nor informed consent was required in this study due to the public-available data of GEO. No ethical approval or informed consent was required in this study owing to the use of publicly available data.

Identification of DEGs
We applied the Limma package to identify differentially expressed circRNAs (DEcircRNAs) with a threshold of |log2 fold change (FC)|> 3 and adjusted P-value < 0.01. The differentially expressed mRNAs (DEmRNAs) and miRNAs (DEmiRNAs) were analyzed using the edgeR package with thresholds of |log2 FC|> 1 and an adjusted P-value < 0.05.

Construction of the ceRNA network
The circRNA-miRNA interactions were predicted using the Circular RNA Interactome (Cir-cInteractome) (https://circinteractome.nia.nih.gov/) databases. These target miRNAs were further screened according to the DEmiRNAs obtained from TCGA database. In addition, miRNA-targeted mRNAs were retrieved from the miRTarBase and TargetScan databases [17,18]. Only mRNAs recognized by both databases were considered to be candidate targets, and were intersected with the identified DEmRNAs to screen out the DEmRNAs targeted by the DEmiRNAs. Based on these DEmiRNA-DEcircRNA and DEmiRNA-DEmRNA interactions, we constructed circRNA-miRNA-mRNA regulatory network, which was visualized using Cytoscape 3.7.0 software.

PPI network construction and analysis
The Search Tool for the Retrieval of Interacting Genes database (Version 10.0, http://stringdb.org) was used to predict potential interactions among DEmRNAs. A combined score of > 0.4 was considered significant. Cytoscape 3.7.0 was used for visualization. We used the cytoHubba application to explore the hub genes of the obtained PPI network [19].

GO annotation and KEGG pathway analysis
The Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the clusterProfiler package of R software [20]. P value less than 0.05 were considered to represent statistically significant enrichment of DEGs in pathways or GO terms.

Construction of the ceRNA network
To further examine the underlying mechanism of circRNA in mediating mRNA based on miRNA, a ceRNA network was established. We predicted the miRNAs targeted by the 27 DEc-ircRNAs using the CircInteractome online database. A total of 660 circRNA-miRNA pairs were identified. After intersecting with the DEmiRNAs, 37 circRNA-miRNA pairs, including 21 circRNAs and 14 DEmiRNAs, remained. We further searched for mRNAs targeted by these 14 DEmiRNAs from the miRTarBase and TargetScan databases, and selected those overlapping with the identified DEmRNAs. Ultimately, a total of 150 DEmRNAs were involved in the ceRNA network, along with 21 circRNAs, and 14 miRNAs (Fig 3).

Construction of the PPI network
After removing unconnected nodes, we established a PPI network that contained 46 nodes and 53 edges ( Fig 4A). To explore the hubgenes in the network, indicating a critical role in the process of BC carcinogenesis, the degree and betweenness centrality were evaluated, and the following top 10 hub_genes were extracted using the cytoHubba app: CDK1, CENPM, CENPF, KNTC1, DSN1, HIST1H2BJ, RAD51, EZH2, EXO1, and BRCA1 ( Fig 4B). Based on this result, we established a circRNA-miRNA-hub_gene subnetwork, including 31 ceRNA regulatory modules. After excluding modules with inconsistent expression of circRNAs and mRNAs, 29 modules remained ( Fig 5).

Functional assessment
To provide insight into the underlying biological processes and pathways related to the DEmRNAs in the ceRNA network, we performed the GO annotation and KEGG pathway analyses. GO analysis showed that the DEmRNAs were, most significantly enriched in the biological process, cellular components, and molecular function terms strand displacement, condensed chromosome, and endodeoxyribonuclease activity, respectively (P<0.05). The top five GO terms e indicated in Table 1. In addition, the DEmRNAs were significantly enriched in two KEGG pathways: homologous recombination and Fanconi anemia pathway.

Discussion
In recent years, several computational models have been developed to identify cancer-related non-coding RNAs, and some of them perform well, such as Inductive Matrix Completion for MiRNA-Disease Association prediction (IMCMDA), Matrix Decomposition and Heterogeneous Graph Inference (MDHGI), and Laplacian Regularized Sparse Subspace Learning for MiRNA-Disease Association prediction (LRSSLMDA) [21][22][23][24]. However, few computational models predict potential associations between circRNAs and bladder cancer. Our study aims to establish a circRNA-miRNA-mRNA regulatory network. CircRNAs were first identified in viruses in the 1970s, and then subsequently discovered to be present in human cell lines and the human body [25,26]. Since circRNAs lack 5 0 or 3 0 polarities or polyadenylated tails, they are stable and cannot be degraded by RNase-R enzyme [7].
CircRNAs are abundant in eukaryotic cells and show a high degree of conservation, along with structural stability, and a certain degree of organization, timing and disease-specific activity [27,28]. Based on these features, circRNAs have become research hot_spots, especially with respect to cancer research. Recent, studies have revealed the abundance and function of cir-cRNAs in tumorigenesis [29,30] and the mechanisms by which circRNAs participate in regulating malignant biological processes [31,32], demonstrating their potential to serve as biomarkers of malignancies [33][34][35]. However, the exact role of circRNAs in BC remains largely unclear. This study represents the first attempt at integrating the differentially expressed circRNAs, miRNAs, and mRNAs in BC from public databases to provide a cir-cRNA-miRNA-mRNA regulatory network.
Indeed, accumulating evidence indicates a role of circRNAs in the initiation and progression of BC [36,37]. Xu et al. [36] analyzed 40 pairs of BC tissue and blood samples and found that circPTK2 was highly expressed in BC. Moreover, they showed that elevated circPTK2 expression could promote the proliferation and migration of BC cells. Chen et al. [37] found that circPRMT5 expression was elevated in BC tissues and was associated with an advanced clinical stage and poor OS of BC patients. They further revealed that circPRMT5 promotes the metastasis of urothelial carcinoma of the bladder through sponging miR-30c to induce the epithelial-mesenchymal transition. In our study, we analyzed BC samples and normal renal samples and they found 21 circRNAs involved in the ceRNA network based on analysis of BC tissues and normal renal samples, including hsa_circ_0000144 and hsa_circ_0023642, whose dysregulated expression has been associated with the pathogenesis and prognosis of BC, indicating their potential as tumor-related biomarker [38,39]. However, none of the other 19 cir-cRNAs in our ceRNA network have been reported in BC to date.
MiRNAs are non-coding single-stranded RNA molecules that consist of approximately 22 nucleotides [40]. The aberrant expression of miRNAs, which regulate the expression of multiple oncogenes and tumor suppressors, has been widely associated with cancer development [41]. In this study, we identified 14 DEmiRNAs in the ceRNA network. Some researchers have studied the binding of circRNAs to miRNAs and their interactions in BC [14,42]. Li et al. https://www.ncbi.nlm.nih.gov/pubmed/30382592 [14] demonstrated that the circRNA Cdr1as expression inhibited BC cell proliferation, apoptosis, and invasion by sponging miR-135a. Yang et al. [42] indicated that circ-ITCH inhibits BC progression by sponging miR-17/miR-224 and regulating p21 and PTEN expression. Of the 14 miRNAs involved in our ceRNA network, four have been reported to play important roles in the initiation and development of BC, including miR-217, miR-375, and miR-431, and miR-935 [38,[43][44][45].
To further identify the key circRNAs participating in the regulatory network, we established the PPI network, and screening out 10 hub_genes, followed by construction of the circRNA-miRNA-hubgene network, including 29 ceRNA regulatory modules. To understand the underlying biological processes and pathways between DEmRNAs in the ceRNA network, we performed the functional enrichment analyses. The result indicated that the DEmRNAs were significantly enriched in homologous recombination and Fanconi anemia pathway. Qiao et al. [46] indicated that Imatinib can radiosensitize BC by targeting homologous recombination. In addition, Madubata et al. [47] found that BC patients with Fanconi anemia nonsense variants display a BRCA-deficiency somatic mutation signature. Therefore, the DEmRNAs are involved in many important BC-associated biological functions and metabolic pathways. In the future, we will do more analysis based on large samples. Second, the conclusions of our study are only based on the current existing tools and databases. These conclusions e further validated by real experiments. In addition, deep learning, such as graph convolutional neural networks (GCN), convolutional neural networks (CNN), has invaded the field of bioinformatics [48,49]. In the future, we will try some analysis based on deep learning. Third, the prognostic value of these DEcircRNAs in BC has not been evaluated. In future studies, we will collect more clinical samples validate our findings and further explore the function of these DEcircR-NAs using in vitro and in vivo experiments.

Conclusions
In this study, we constructed the first ceRNA network in bladder cancer, by identifying DEmiRNAs, DEmRNAs, and DEcircRNAs between BC and normal tissues from data in public databases. Based on these interactions, we identified miRNA and circRNA targets and constructed a protein interaction network to highlight hub genes with likely involvement in BC pathogenesis. We further conducted functional analyses of the DEGs to provide insight into the key biological processes involved in this regulation network. We believe that our study makes a significant contribution to the literature because the roles of circRNAs as potential miRNA "sponges" in several types of cancers are increasingly being emphasized, including BC. Our constructed networks and corresponding modules can serve as a useful guide for further targeted research into the molecular pathogenesis of BC to identify new therapeutic targets and/or biomarkers.
Supporting information S1