Identification of the Key Factors Related to Bladder Cancer by lncRNA-miRNA-mRNA Three-Layer Network

Bladder cancer is the most common malignant tumor of the urinary system, and it has high incidence, high degree of malignancy, and easy recurrence after surgery. The etiology and pathogenesis of bladder cancer are not fully understood, but more and more studies have shown that its development may be regulated by some core molecules. To identify key molecules in bladder cancer, we constructed a three-layer network by merging lncRNA-miRNA regulatory network, miRNA-mRNA regulatory network, and lncRNA-mRNA coexpression network, and further analyzed the topology attributes of the network including the degree, betweenness centrality and closeness centrality of nodes. We found that miRNA-93 and miRNA-195 are controllers for a three-layer network and regulators of numerous target genes associated with bladder cancer. Functional enrichment analysis of their target mRNAs revealed that miRNA-93 and miRNA-195 may be closely related to bladder cancer by disturbing the homeostasis of the cell cycle or HTLV-I infection. In addition, since E2F1 and E2F2 are enriched in various KEGG signaling pathways, we conclude that they are important target genes of miRNA-93, and participate in the apoptotic process by forming a complex with a certain protein or transcription factor activity, sequence-specific DNA binding in bladder cancer. Similarly, AKT3 is an important target gene of miRNA-195, its expression is associated with PI3K-Akt-mTOR signaling pathway and AMPK-mTOR signaling pathway. Therefore, we speculate that AKT3 may participate in proliferation and apoptosis of bladder cancer cells through these pathways, and ultimately affect the biological behavior of tumor cells. Furthermore, through survival analysis, we found that miRNA-195 and miRNA-93 are associated with poor prognosis of bladder cancer. And the Kaplan-Meier curve showed that 24 mRNAs and nine lncRNAs are closely related to overall survival of bladder cancer.


INTRODUCTION
Bladder cancer (BC) is a common malignant tumor with the tenth incidence rate and the ninth mortality rate in the world (Bray et al., 2018). As the most common malignant tumor of the urinary system, its etiology and pathogenesis are still not very clear. It has some biological characteristics such as high incidence, high degree of malignancy and easy recurrence after surgery, and it shows a clear upward trend, which is a serious threat to human health and social development (Knollman et al., 2015;Tian et al., 2015). The etiology of BC is complex, including both intrinsic genetic factors and external environmental factors. At present, it has been clarified that the occurrence and development of BC are related to smoking, environmental pollution, and other factors (Kiriluk et al., 2012;Van Roekel et al., 2013;Smith et al., 2016;Wilcox et al., 2016). But its molecular mechanism is still not completely clear. Therefore, the research and clarification of its cause of disease and the specific molecular mechanism will effectively promote the diagnosis and treatment of BC.
In recent years, a large number of studies on single molecules of BC have found that dysregulated mRNA contributes to the occurrence of BC tumors by participating in abnormal activation of many pathways such as Akt signal pathway and NF-kB pathway Tsui et al., 2018). Abnormal expression of lncRNA UCA1 and H19 is closely related to the initiation, development, and metastasis of BC (Zhang et al., 2013). Studies on the relationship between urinary tumors and miR-145, miR-205 indicate that a large number of miRNAs have been involved in tumor progression (Catto et al., 2011;Simonato et al., 2013;Braicu et al., 2015;Chen et al., 2017). Some people also analyzed the differentially expressed miRNAs and their target genes in BC from the perspective of the network, and it is pointed out that miRNAs participate in the occurrence and development of BC by regulating their target genes, and are related to BC recurrence (Guancial et al., 2014;Jin et al., 2015;Li et al., 2017;Zhao et al., 2017). LncRNAs and mRNAs can form complex coexpression networks in various diseases including BC, and promote the development of BC by regulating proteincoding genes (Cao et al., 2018b;Li et al., 2018;Chen et al., 2019a;Chen et al., 2019b). In addition, the ceRNA hypothesis suggests that lncRNA can not only directly participate in the expression of target regulatory genes, but also may adsorb target miRNAs, affect the binding of miRNAs to mRNAs, and further intervene in mRNA expression Lv et al., 2017;Wang et al., 2017;Miao et al., 2019;Teng et al., 2019).
Since the mRNA-lncRNA coexpression network and the miRNA-related regulatory network are in the common system, it is one-sided to study the key molecules of BC by the singlelayer network method. In contrast, converged network analysis can reflect the importance of molecules in the network more accurately and more directly. Therefore, it is considered as a feasible method to further analyze the topology of the network according to the degree, betweenness centrality and closeness centrality of the nodes in the network to find and predict the core molecules potentially related to the development of BC (Greenbury et al., 2010). In this paper, we constructed a three-layer network by merging lncRNA-miRNA regulatory network, miRNA-mRNA regulatory network and lncRNA-mRNA coexpression network, and further analyzed the topology attributes of nodes in the network.

MATERIALS AND METHODS
The main steps of identification of the key factors related to BC by lncRNA-miRNA-mRNA three-layer network include: data collection and pretreatment, construct three-layer network, calculate topology properties, identify the key nodes, functional analysis, and survival analysis of the key nodes. The flowchart of data collection and method implementation was shown in Figure 1. Next, we will describe a more detailed process of identification of the key factors related to BC by lncRNA-miRNA-mRNA three-layer network.
Using R software, the DEmRNAs and DElncRNAs coexpression data were integrated into a matrix. Then, the lncRNA-mRNA interactions were predicted according to the method of weighted gene coexpression network by using the WGCNA package in R and the correlation coefficient is set to 0.8.
Finally, we merged all miRNA-mRNA pairs and miRNA-lncRNA pairs and mRNA-lncRNA pairs to construct a converged network. We defined the miRNA, mRNA, and lncRNA as the nodes of the network. If they interacted each other, we added an edge between them, then the converged three-layer network was constructed. Cytoscape software (https://cytoscape.org/) was used to visualize the network.

Calculation of Topology Attributes
Further calculate the topology attribute of the network according to the degree, betweenness centrality, and closeness centrality of the nodes in the three-layer network. For a given graph, G = (V, E), V and E represent a set of nodes and edges, respectively, and degree (D) is usually used to reflect the basic topology attributes of the network, indicating the number of nodes connected to neighboring nodes, assuming K is the number of edges connected to node v, and described as follows: The betweenness centrality is one of the most effective methods for evaluating the importance of a node in the network. It is defined as the proportion of the shortest path through a node in the network to the total number of shortest paths, reflecting the influence of the node in the network. It plays an important role in the process of information transfer. And it can be seen from the definition that if a node is the only way for other nodes in the network to communicate, the node has an important position in the network. The higher the node betweenness centrality value, the greater its influence, the expression is as follows: where 's st (v)' represents the number of the shortest path between nodes 's' and 't' passing through node 'v,' and 's st ' represents the number of shortest paths between node pairs 's' and 't.' Closeness centrality reflects the proximity of a node in the network to other nodes. Closeness centrality of a node 'v' is the reciprocal of the average shortest path distance from node 'v' to all other nodes. That is, the closer a node is to other nodes, the greater its closeness centrality. Described as follows: where 'd iv ' is the shortest-path distance between 'i' and 'v,' and '|V|' is the number of nodes in the network.

Extract a Subnetwork From Three-Layer Network
Find important nodes whose topological attributes including degree, betweenness centrality, and closeness centrality ranked in the top ten in the network, namely hub nodes. Then extract the hub nodes and their directly connected RNAs in a three-layer network, and construct a subnetwork related to hub nodes.

Enrichment Analysis
To achieve enrichment analysis of candidate genes in BC, we used DAVID 6.8(https://david.ncifcrf.gov/) for gene ontology (GO) analysis, including biological process (BP) analysis, cellular component (CC) analysis, and molecular function (MF) analysis, with p < 0.05 as a screening condition to find significantly enriched GO terms. Then, in order to more fully understand the situation of pathway enrichment, KOBAS 3.0 (http://kobas. cbi.pku.edu.cn/) was used as the metabolic pathway KEGG enrichment analysis, at P < 0.01, the number of genes enriched in the pathway was greater than or equal to 2 as a screening condition. By analyzing the significantly enriched GO biological process classification and KEGG pathways, we can further predict the biological function of DEmRNAs from gene function and pathway function, respectively.

Survival Analysis
According to the clinical information downloaded from TCGA database and the expression level of differentially expressed RNAs, the median was used as the cut-off point, and the RNAs expression was divided into two groups: high expression and low expression. Then, we performed Kaplan-Meier survival analysis on the differentially expressed RNAs, the RNAs in the three-layer converged network, and the RNAs in the subnetwork. The "survival" Package in R software was used to plot the survival curves and view RNAs that is significantly associated with the patient's survival prognosis. P-value < 0.05 was considered statistically significant.

Construction of a Three-Layer Network
The lncRNA-mRNA coexpression network included 75 lncRNAs, 169 mRNAs, and 380 edges. The lncRNA-miRNA regulatory network constructed in the present study exhibited 129 lncRNAs and 35 miRNAs and 882 edges. And the regulatory network between miRNAs and target genes was established, which included 201 mRNAs, 56 miRNAs, and 487 edges.
Then, we constructed a three-layer network by merging lncRNA-miRNA regulatory network, miRNA-mRNA regulatory network, and lncRNA-mRNA coexpression network, which included 367 mRNAs, 65 miRNAs, 197 lncRNAs, and 1749 edges. The three-layer network was visualized with Cytoscape software and was shown in Figure 3.

Analysis of Topology Attributes
Further analyze the topology of the network and find out the important nodes in the network according to the degree, closeness centrality, and betweenness centrality of the nodes. There are two miRNAs (miRNA-93 and miRNA-195) appeared in each dimension, indicating that they have higher degree, betweenness centrality, and closeness centrality and played crucial roles in the three-layer network. The topology attributes of the top 10 genes were shown in Table 1.

Functional Analysis
We extracted two key miRNAs and their linked mRNAs and lncRNAs from a three-layer network. Then rebuilt a new subnetwork, which included 71 mRNAs, 2 miRNAs, 73 lncRNAs, and 144 edges. The subnetwork was visualized with Cytoscape software and was shown in Figure 4.
Through GO analysis of 35 target genes of miRNA-93 ( Figure  5A), we found that the dysregulated mRNAs associated with several -BP-, including lens fiber cell apoptotic process and plusend-directed vesicle transport along microtubule. And they are also related to nucleus and nucleoplasm in -CC-. Furthermore, based on -MF-analysis, we conclude that miRNA-93 may be involved in the biological processes of BC by regulating the expression of its target genes. The functions of its target gene are mainly glucose homeostasis, protein binding and MAP kinase tyrosine/serine/threonine phosphatase activity and transcription factor activity, sequence-specific DNA binding. Among them, the discovery of MAP kinase serine phosphatase activity in molecular function further confirms the previously reported results that miRNA-93 promotes BC cell proliferation and invasion by targeting PEDF, a member of the serine protease inhibitor superfamily (Jiang et al., 2019). The details of the GO term associated with 35 DEmRNAs are shown in Table 2.
We also performed GO analysis on 36 target genes of miRNA-195 ( Figure 5B), and functional enrichment results showed that these differentially expressed genes mainly (30/36) focused on protein binding in molecular function. At the same time, we found that the -CC-associated with mRNA is nucleoplasm (p = 0.012459971), and the -BP-most strongly associated with mRNA is the regulation of cell cycle (p = 0.002164798). In addition, the most relevant and only relevant -MF-associated with mRNA is the protein binding (P = 0.000145982). Based on these data, it can be speculated that the target genes of miRNA-195 can form complexes with different proteins to participate in regulating the cell cycle and positive regulation of endothelial cell migration, somatic stem cell population maintenance, stem cell development, accordingly, have certain effects on the occurrence and development of BC. The details of the GO term associated with 36 DEmRNAs are displayed in Table 3.
KEGG enrichment analysis of metabolic pathways for miRNA-93 and miRNA-195 targeted mRNAs showed that a total of 28 genes are enriched in 39 pathways, of which 10 target genes of miRNA-93 are mainly involved in the regulation of 16 KEGG signaling pathways including cancer pathways and Hepatitis B and Axon guidance and other pathways, the 19 target genes of miRNA195 are mainly involved in the regulation of 31 KEGG signaling pathways including many cancer pathways and Cell cycle, HTLV-I infection ( Figure 6). Demonstrating that miRNAs can play a role in a variety of BC signaling pathways. In addition, we found that the signaling pathways "MicroRNAs in cancer," "HTLV-I infection," and "cell cycle" are significantly enriched KEGG pathways, indicating that these three signaling pathways are closely related to BC. The details of the KEGG pathway enrichment are displayed in Table 4.
Since E2F1 and E2F2 are enriched in various KEGG signaling pathways, they may be important target genes of miRNA-93. Similarly, AKT3 may be an important target gene of miRNA-195. From the functions annotated in the genes E2F1 and E2F2 in GO, we can speculate that the genes E2F1 and E2F2 are present on the nucleoplasm, and participate in the lens fiber cell apoptosis process or intrinsic apoptotic signaling pathway by p53 class mediator by forming a complex with a certain protein or transcription factor activity, sequence-specific DNA binding, which is consistent with the results reported previously (Castillo et al., 2015). In addition, PI3K-Akt and AMPK are reported to be positive and negative regulators of mTOR activity, respectively, and mTOR is a well-known sensor for extracellular nutrients and growth factors, which converges many key signals and is widely involved in cell cycle, cell proliferation, and cellular metabolic processes, and has important biological roles in regulating apoptosis, physiological, and pathological processes (Ciuffreda et al., 2010). The PI3K-Akt-mTOR signaling pathway has been shown to play a central role in several cellular functions, including survival, proliferation, adhesion, migration, differentiation, and metabolism (Gao et al., 2003). Activated AMPK inhibits mTOR and further inhibits protein synthesis, cell proliferation, cell cycle progression, and angiogenesis (Jalving et al., 2010).This study found that the expression of AKT3 is associated with these pathways, suggesting that AKT3 may participate in important processes such as cell proliferation and apoptosis through these pathways, ultimately affecting the biological behavior of tumor cells.
KEGG analysis of the downregulated target genes of miRNA-93 revealed that EPHA7, CFL2, and DPYSL2 are significantly enriched in axon guidance pathway (P = 0.000129528), and PPP1R3B and SLC2A4 are significantly attached to insulin resistance pathway (p = 0.001718367), we hypothesized that the activated axon guidance pathway and the insulin resistance pathway may be responsible for the downregulation of miRNA-93 in BC. KEGG analysis of the upregulated target genes of miRNA-195 revealed that CHEK1, CCNE1, and CDC25A are significantly enriched in the cell cycle pathway (p = 0.00000888), and CHEK1, MYB, and WNT7A are enriched in HTLV-I infection pathway (p = 0.000077), CCNE1, KIF23, and CDC25A enriches microRNAs in cancer (p = 0.000117168), CHEK1 and CCNE1 enriches p53 signaling pathway (p = 0.000241969). Therefore, we believe that "cell cycle, HTLV-I infection, MicroRNAs in cancer, p53 signaling pathway" may be related to downregulation of miRNA-195 in BC.

Survival Analysis
We performed Kaplan-Meier survival analysis to predict overall survival (OS) of DERNAs in BC. As a result, we found that 253 of the 1457 differentially expressed lncRNAs were related to survival, 38 of the 247 differentially expressed miRNAs were associated with survival, and 625 of the 3047 differentially expressed mRNAs were related to survival (Supplementary  Material Tables S1-S3). Similarly, Kaplan-Meier survival analysis was performed on miRNAs/lncRNAs and mRNAs in the three-layer converged network, and we found that 39 of the   Tables S4-S6).
In addition, we also performed Kaplan-Meier survival analysis on miRNAs/lncRNAs and mRNAs in the subnetwork. Among them, miRNA-93 is highly correlated with overall survival of BC ( Figure 7A). A large number of studies have shown that miRNA-93 can be widely involved in various biological behaviors of tumor cells and affect the development of tumors (Li et al., 2014). Therefore, in BC, the high expression of miRNA-93 may promote the proliferation of cancer cells and inhibit its apoptosis, which can lead to excessive proliferation and deterioration of cancer cells. miRNA-195 also has a correlation with the overall survival of BC ( Figure 7B). Related studies (Fei et al., 2012) showed that miRNA-195 is a tumor suppressor gene, which inhibits cell proliferation and promotes apoptosis when normal expression, and downregulation leads to cancer, while miRNA -195 is downregulated in BC, so it may be closely related to the occurrence and development of BC. And we also found that nine lncRNAs (AC127496.3, ADAMTS9-AS2, C7orf65, LINC00536, MAGI2-AS3, SACS-AS1, AC016773.1, AC112721.1, and HCG22) were significantly associated with survival outcomes in patients with BC, among which MAGI2-AS3, C7orf65, and LINC00536 are lncRNAs linked to both miRNAs in the network. Twenty-four mRNAs (KATNAL1,  ANKRD29, CFL2, DPYSL2, DUSP2, EGR2, FAM129A,  KPNA2, LDLR, RRM2, AKT3, AMOTL1, BTG2, CPEB2,  FASN, FGF2, LSM11, PTPRD, RAB23, RS1, RUNX1T1, TMEM100, VCL, and ZFHX4) are closely related to the survival prognosis of patients, 10 of which are target genes of miRNA-93, and 15 are target genes of miRNA-195. And KATNAL1 is a mRNA that is simultaneously regulated by two miRNAs in a network related to BC (Supplementary Material Figure S1).

DISCUSSION
A large number of studies have shown that the development of BC is highly correlated with the abnormal expression of noncoding RNA and protein-coding genes. However, there are still few studies on the key molecules of BC and abnormal gene expression patterns and pathways in BC. In this study, we constructed a three-layer network by merging lncRNA-miRNA regulatory network, miRNA-mRNA regulatory network, and lncRNA-mRNA coexpression network, and calculated the topology attributes of each node in the network, including degree, closeness centrality, and betweenness centrality. Our study suggested that miRNA-93 and miRNA-195 are controllers of the three-layer network associated with BC and a regulator of numerous target genes, and their dysregulation may be closely related to the pathogenesis of BC. A recent study also showed that miRNA-93 is upregulated in BC and is associated with tumor stage and lymph node metastasis (Jiang et al., 2019). And studies have reported that miRNA-195 inhibits glucose uptake and proliferation of BC T24 cells by inhibiting GLUT3 expression (Fei et al., 2012), another study found that miRNA-195 inhibits BC cell proliferation by inhibiting Cdc42 / STAT3 signaling .
To clarify the significance of miRNA-93 and miRNA-195 in BC as accurately as possible, we extracted mRNAs and lncRNAs directly linked to miRNA-195 and miRNA-93 from a three-layer network, and built a new subnetwork. Then, we assigned miRNA-targeted mRNA to terms in GO and known pathways FIGURE 6 | miRNA-targeted mRNA enrichment pathway. The x-axis represents the genes in the network, in which the red solid line is the target genes of miRNA-93, and the black dotted line is the target genes of miRNA-195. The y-axis represents the KEGG pathway, the red solid box is the enrichment pathway of the miRNA-93 target genes, and the black dotted box is the enrichment pathway of the miRNA-195 target genes.
in the KEGG database. Through GO analysis, we found that DEmRNAs in the network is mainly concentrated in the protein binding in molecular function, suggesting that the proteincoding genes play vital roles in the development of BC. What's more, the result from KEGG pathway analysis showed that dysregulated mRNAs contribute to the development of BC through abnormal activation of many pathways. Therefore, we concluded that miRNA-93 and miRNA-195 may participate in various biological processes by regulating the expression of target genes with different functions in BC, which may have an important impact on the occurrence and development of BC. In addition, in recent years, a large number of studies have shown that miRNA-93 and miRNA-195 play an important role in the development of tumors. miRNA-93 overexpression was associated with tumor progression, metastasis and poor prognosis in patients with head and neck squamous cell carcinoma .And inhibiting the expression of miRNA-93 can inhibit the proliferation and migration of gastric cancer and liver cancer cells (Xu et al., 2012;Zhang et al., 2016). In addition, studies have shown that miR-195 has tumor suppressive properties in gastric cancer and glioma (Deng et al., 2013;Hui et al., 2013). Some evidences indicate that miRNA-195 can inhibit the ability of liver cancer cells to promote endothelial cell migration and angiogenesis, and can directly inhibit the migration and invasion of liver cancer cells, and inhibit liver cancer. Some evidences indicate that miRNA-195 can inhibit the ability of hepatoma cells to promote endothelial cell migration and angiogenesis, and can directly inhibit the migration and invasion of hepatoma cells, thereby inhibiting liver cancer .In osteosarcoma, miRNA-195 inhibits tumor migration and invasion by acting on fatty acid synthase (Mao et al., 2012).Based on these reports and current findings, we believe that the study of the mechanisms of action of miRNA-93 and miRNA-195 has profound implications for the clinical treatment of tumors. Through survival analysis, we identified several RNAs associated with the prognosis of BC, some of which have not been reported. It is worth noting that the miRNAs of forming the hub nodes are related to survival. Previous studies have shown that the abnormal expression of miRNA-93 and miRNA-195 can promote the proliferation of cancer cells and inhibit their apoptosis (Fei et al., 2012;Li et al., 2014). Therefore, they may be involved in and influence the occurrence and development of BC. In addition, some of the genes in the subnetwork are often present in different cancers. For example, studies have found that ADAMTS9-AS2 is significantly downregulated in glioma tissues compared with normal tissues, and negatively correlated with prognosis (Yao et al., 2014). Downregulation of ADAMTS9-AS2 is associated with poor prognosis in patients with gastric cancer (Cao et al., 2018a). Low expression of ADAMTS9-AS2 is associated with overall survival in patients with ovarian cancer (Wang et al., 2018a). And ADAMTS9-AS2 is significantly upregulated in TSCC tissues of patients with lymph node metastasis and is closely related to poor prognosis (Li et al., 2019b). Besides, there are reports that the downregulation of MAGI2-AS3 in BC patients is associated with poor prognosis (Wang et al., 2018b) and the high expression of LINC00536 is negatively correlated with the survival rate of patients with BC (Li et al., 2019a), this is consistent with our results. These reports help identify and improve the accuracy of prognostic predictions.
In this work, we considered the relations of differentially expressed genes to construct the networks. As we know, the human body is a complex system. Maybe, there are some genes that do not differentially express have effects on the oncogenesis and progression of BC. In addition, we analyzed how the key genes affect other genes. It is meaningful to study the relation between key genes with proteins or other chemical molecules from the perspective of pathway function. We are currently working to integrate various valid data sets and try to develop more reliable prediction methods to provide reasonable and valuable candidate factors for the early diagnosis and treatment of BC.

CONCLUSION
Taken together, we found that miRNA-93 and miRNA-195 are key factors for the development of BC by calculating the topology attributes of the nodes in the three-layer network. In addition, many mRNAs with different functions were found from KEGG and GO analysis, of which E2F1 and E2F2 are important target genes of miRNA-93, and AKT3 is an important target gene of miRNA-195, their dysregulation may be closely related to cell proliferation and apoptosis. Further the Kaplan-Meier curve analysis revealed that miRNA-93 and miRNA-195 are associated with poor prognosis of BC. In summary, our results have great significance in clinical practice, which can provide a reference for a better understanding of the occurrence and development of BC.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study can be found in the Supplementary Material.

AUTHOR CONTRIBUTIONS
XW and YD designed the research ideas. XW collected and organized the data. XW, YD, JW, and YW analyzed the data. XW wrote the manuscript.

FUNDING
This work was funded by The National Natural Science Foundation of China (grant numbers 21541006, 61772241).