To construct a ceRNA regulatory network as prognostic biomarkers for bladder cancer.

Emerging evidence demonstrates that competing endogenous RNA (ceRNA) hypothesis has played a role in molecular biological mechanisms of cancer occurrence and development. But the effect of ceRNA network in bladder cancer (BC), especially lncRNA-miRNA-mRNA regulatory network of BC, was not completely expounded. By means of The Cancer Genome Atlas (TCGA) database, we compared the expression of RNA sequencing (RNA-Seq) data between 19 normal bladder tissue and 414 primary bladder tumours. Then, weighted gene co-expression network analysis (WGCNA) was conducted to analyse the correlation between two sets of genes with traits. Interactions between miRNAs, lncRNAs and target mRNAs were predicted by MiRcode, miRDB, starBase, miRTarBase and TargetScan. Next, by univariate Cox regression and LASSO regression analysis, the 86 mRNAs obtained by prediction were used to construct a prognostic model which contained 4 mRNAs (ACTC1 + FAM129A + OSBPL10 + EPHA2). Then, by the 4 mRNAs in the prognostic model, a ceRNA regulatory network with 48 lncRNAs, 14 miRNAs and 4 mRNAs was constructed. To sum up, the ceRNA network can further explore gene regulation and predict the prognosis of BC patients. ©2020 Zhongnan Hospital of Wuhan University. Journal of Cellular and Molecular Medicine published by Foundation for Cellular and Molecular Medicine and John Wiley & Sons Ltd.

INTRODUCTION

Bladder cancer (BC) is among the leading cancers worldwide. For 30 years, the prognosis and treatment for BC had not made great progress. , The poor prognosis and the difficulty in treatment are bound to the alteration and regulation of genes. An alteration in genes that play a key role in cell cycle regulation can contribute to high‐grade tumours and may alter the sensitivity to drugs in BC. Recent studies demonstrate that the regulation of ceRNA network plays a crucial role in the occurrence and development of cancers. ,

The ceRNA hypothesis postulates that in addition to the conventional function, miRNAs targeting RNAs, a reversed logic exists. Any RNA transcript with miRNA response elements (MREs) could bind to miRNA, thereby regulating the expression of RNAs with the same MREs. Significantly, because the two transcripts share the same MREs, we can predict the ability of two transcripts to regulate each other based on MRE overlap. , There is increasing evidence that lncRNA, miRNA and mRNA are involved in the biological processes of cancer. Among the ceRNA, lncRNAs and mRNAs might harbour the same MREs. Therefore, when miRNA binds to MRE on lncRNA, mRNA expression is not inhibited which may contribute to disease pathogenesis. Currently, the ceRNA network has been shown to play an important role in the occurrence and progression of different types of cancers. But the roles of ceRNA network in BC, especially the analysis of ceRNA regulatory network based on large sample size and next‐generation sequencing, were not completely expounded. So, this research aims to construct a ceRNA regulatory network in BC.

METHODS

RNA‐Seq data and clinical characteristics of patients from TCGA

Figure 1 shows the flow diagram of our study. We obtained RNA‐Seq data of BC from the TCGA database (https://portal.gdc.cancer.gov/), which contained 414 bladder urothelial carcinoma samples and 19 normal bladder samples. Besides, TCGA also included the corresponding clinical characteristics (406 bladder urothelial carcinoma patients) including overall survival and other clinicopathological characteristics of BC patients.

Figure 1

The flow diagram of this study

Identification of differentially expressed mRNAs and lncRNAs

We transformed the ensemble IDs of BC samples into gene symbols based on GENCODE (https://www.gencodegenes.org/human/). We identified differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) between bladder urothelial carcinoma and normal bladder tissue by means of EdgeR R package. And when they met these criteria (Adjust P < .05, and |log2 FC| ≥ 1), we considered mRNAs and lncRNAs as DEmRNAs and DElncRNAs.

Analysis of gene function and pathway enrichment

We performed Gene Ontology (GO) and Kyoto Encyclopedia of genes and Genomes (KEGG) and Gene Set Enrichment Analysis (GSEA) for DEmRNAs by means of clusterProfiler R package. We take advantage of GO to describe gene functions, including three aspects: biological process (BP), molecular function (MF) and cellular component (CC). KEGG‐GSEA was used to get possible pathways, and the cut‐off criteria for gene sets were set at P < .05.

Construction of a co‐expression network

WGCNA R package was conducted to construct a co‐expression network. First, unqualified genes were removed by goodSamplesgenes. Then, construct a sample network based on squared Euclidean distance, and samples with Z.k < −2.5 which were considered as outlying samples were excluded from subsequent studies. Then, the variances of mRNAs and lncRNAs across the left samples were calculated, and we chose the top 60% mRNAs and lncRNAs with highest variances for WGCNA. Next based on the scale‐free topology criterion, we chose an appropriate soft‐thresholding power (β) by pickSoftThreshold to set up a weighted adjacency matrix (WAM). Given the function of the topological overlap matrix (TOM), we transformed adjacency into TOM. After calculating module eigengenes, densely correlated and highly co‐expressed genes with BC patients are grouped into the same module.

Construction of the ceRNA network

LncRNAs in the key modules most relevant with BC were picked to predict miRNAs by the means of miRcode (http://www.mircode.org). Then, we explore target mRNAs of miRNAs by the MiRDB (http://www.mirdb.org/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw//), StarBase (http://starbase.sysu.edu.cn/) and TargetScan (http://www.targetscan.org//). Besides, Cytoscape software (https://cytoscape.org/) was used to visualize the relationship of ceRNA network.

Construction of BC prognostic signatures

We collected the clinical data of 186 BC patients with both OS and DFS. At the same time, for verification, we randomly divided patients into two groups based on the ratio of 5:4 (the training set = 103 and internal test set = 83). We conducted the chi‐square test to analyse the differences in clinical features including age, gender, tumour subtype, differentiation grade, T stage, TNM stage and survival status between training set and internal testing set(the N stage and M stage are discarded because of patients with DFS data only being N0 or M0 status). By means of the ‘survival’ package in R, univariate Cox regression method was conducted to explore the correlation between the mRNAs obtained by prediction and overall survival (OS) of BC patients. By means of LASSO regression analysis, we set up a survival model. And we used the following formula to calculate the prognostic risk score of patients: Risk score = βmRNA1 * expmRNA1 + βmRNA2 * expmRNA2 + … + βmRNAn * expmRNAn (the lower the risk score, the lower the risk of recurrence, where ‘β’ is the regression coefficient of mRNAs obtained by LASSO regression analysis, and ‘exp’ is the expression of corresponding mRNAs). Based on the survivalROC R package, a  receiver operating characteristic (ROC) curve was drawn with risk score against survival status. At the same time, the cut‐off point of risk score was picked by survminer R package which divided patients into high‐ and low‐risk groups. Then, chi‐square test was employed to investigate the correlations between the risk score of the 4‐gene‐based classifier with OS and clinicopathological characteristics.

RESULTS

Analysis of DEmRNAs of TCGA data

We identified the DEmRNAs in the TCGA database of BC which contained 841 up‐regulated and 1408 down‐regulated mRNAs. The distribution of the DEmRNAs is shown in Figure 2A. To explore the function of DEmRNAs and their molecular mechanisms in cell cycle progression, we performed GO enrichment analysis on DEmRNAs. GO analysis results showed that DEmRNAs were enriched in muscle system process, muscle contraction and muscle organ development in BP (Figure 2B). Figure 2C depicts the linkages of DEmRNAs and biological terms as a network. In addition, KEGG‐GSEA was used to visualize the distribution of the gene set and the enrichment score. Figure 2D shows DEmRNAs might act through adrenergic signalling in cardiomyocytes, calcium signalling pathway and hypertrophic cardiomyopathy in the development of BC.

Figure 2

Identification of DEmRNAs from TCGA data. A, Volcano map of DEmRNAs displaying 841 up‐regulated and 1408 down‐regulated mRNAs. B, GO analysis of DEmRNAs showed they may act in muscle system process, muscle contraction and muscle organ development. C, The linkages of DEmRNAs and biological processes. D, KEGG‐GSEA pathway enrichment of DEmRNAs

Identification of key modules by WGCNA

Then, we used WGCNA to identify gene co‐expression network modules. After removing unqualified genes and outlying samples, a total of 385 BC and 17 normal samples with 7836 genes (top 60% variances) were left for subsequent analysis. Soft power 8 was picked as the soft thresholding to set up a weighted adjacency matrix (scale‐free R2 = 0.86). After the co‐expressed genes were included in a module and the highly correlated modules were merged (the relevance of modules > 0.75), 25 modules were generated (Figure 3A). Among the 25 modules, the module eigengenes (ME) of the yellow and lightcyan modules are most positively related to traits (BC and normal) (Figure 3B). So, the yellow and lightcyan modules were considered as the key module because of high correlation with trait. Figure 3C‐D shows the correlation between gene significance (GS) and module membership (MM) in the lightcyan and yellow modules. Besides, the heatmap plot of selected genes is shown in Figure 3E.

Figure 3

Analysis of key modules highly related to traits. A, Cluster dendrogram of genes with top 60% variances in the co‐expression network. B, The correlation between modules and traits was displayed, and the ME of the yellow and lightcyan modules are most positively related to traits. C‐D, The correlation between GS and MM in the lightcyan and yellow modules. E, TOM plot of selected genes. Modules are sorted by cluster dendrogram. Red colour represents high topological overlap and light colour represents lower overlap. This TOM plot showed that there is no high topological overlap among 25 modules, indicating that the modules have a high degree of scale independence. F, GO‐GSEA showed the linkages of genes and BP. G, KEGG showed the pathway enrichment of genes in yellow and lightcyan modules

As shown in Figure 3F, we performed GO‐GSEA on the genes in yellow and lightcyan modules to explore their functions in BP. The genes in yellow and lightcyan modules were most related to extracellular matrix and structure organization, and muscle contraction. In addition, genes were highly enriched in focal adhesion, PI3K‐Akt and MAPK signalling pathway, and proteoglycans in cancer by KEGG analysis (Figure 3G), indicating that they might be involved in the process of BC.

Construction of a lncRNA co‐expression network

Then, the construction of a lncRNA co‐expression network was aimed to identify the key modules. The top 60% variance lncRNAs (8604 lncRNAs) were selected for WGCNA. Soft power 4 was selected as the soft thresholding to set up a weighted adjacency matrix (scale‐free R 2 = .89). And a total of 33 modules were generated (Figure 4A). Furthermore, Figure 4B‐C shows that pink module which contained 357 lncRNAs was regarded as key module based on correlation analysis. The network heatmap is shown in Figure 4D.

Figure 4

Identification of lncRNAs modules highly related to traits. A, Cluster dendrogram of lncRNAs with top 60% variances in the co‐expression network. B, The correlation between modules and traits were displayed. C, The correlation between GS and MM in the pink module. D, TOM plot of selected lncRNAs

Overlapped mRNAs of predicted mRNAs, WGCNA mRNAs and DEmRNAs

By construction of lncRNA co‐expression network, we obtained 357 lncRNAs in the pink module. Then, by miRcode online tool, the 357 lncRNAs were used to predict miRNAs which lncRNAs can bind to them through MREs. At the same time, we obtained the top 400 miRNAs with the highest expression through TCGA miRNA‐Seq. After that, we selected the overlapped miRNAs (55) among the predicted miRNAs (208) and the highest expressed miRNAs (400). And by the means of miRDB, TargetScan, miRTarBase and the starBase dataset, we took advantage of these 55 overlapped miRNAs for the second prediction and obtained 2266 mRNAs. Figure 5A‐B shows that we obtained 86 overlapped mRNAs by analysing the 2266 predicted mRNAs, 2,249 differentially expressed mRNAs (841 up‐regulated and 1408 down‐regulated mRNAs) and the 415 mRNAs in the yellow and lightcyan modules. The expression of 86 overlapped mRNAs, to our surprise, was all down‐regulated. Figure 5C shows the expression heatmap of the 86 genes in 433 samples. The GO analysis was further used to display biological function of 86 mRNAs. Figure 5D shows that the mRNAs were mostly related to ameboidal‐type cell migration, muscle organ development and muscle tissue development. In addition, genes were highly enriched in MAPK signalling pathway, proteoglycans in cancer, parathyroid hormone synthesis, secretion and action by KEGG analysis (Figure 5E).

Figure 5

86 mRNAs were displayed by GO and KEGG. A‐B, Common 86 mRNAs among the predicted mRNAs, differentially expressed mRNAs and the mRNAs in the yellow and lightcyan modules. C, The expression heatmap showed the 86 genes were down‐regulated in the patients with tumour. D, The linkages of 86 mRNAs and biological terms in BP were shown. E, KEGG analysis showed the pathway enrichment of 86 mRNAs

Construction of the 4‐gene‐based prognostic model

A total of 186 patients from TCGA were randomly divided into training set (n = 103) and internal testing set (n = 83) separately. Demographic and clinical data for the training and internal testing set are summarized in Table 1. As shown in Table 1, there is no significant difference in the clinicopathological characteristics. Then, we conducted univariate regression analysis for the 86 genes to explore the genes relevant to the OS of BC patients from the training set; then, a 7‐gene‐based model was constructed with the cut‐off for relevance set at P < .05. All the seven genes were then performed with the LASSO regression, and a prognostic model for 1,3,5 years of OS with 4 genes was constructed: ACTC1 + FAM129A + OSBPL10 + EPHA2 (Table 2). Figure 6A shows LASSO coefficient of 4‐gene. After the risk score was calculated for each BC patient, according to the risk scores, we classified the patients into the high‐risk group (n = 23) and low‐risk group (n = 80) with the cut‐off for risk score set at ‘0.3073’. (Figure 6B). Besides, Figure 6C shows the heatmap of the 4 genes expression between the high‐risk and low‐risk group. Kaplan‐Meier survival curves of OS and DFS in the whole sets (n = 186) showed that compared to the low‐risk group, the high‐risk group was tended to lower OS survival rate (P < .001), and DFS survival rate (P = .04) (Figure 6D,E).

Table 1

Clinical features for the BC patients in the training set, testing set

Parameters	Training set n = 103	Testing set n = 83	χ 2	P value
Age
≤65	41	43	2.2105	.1371
>65	62	40
Gender
Male	80	59	0.73567	.3911
Female	23	24
Subtype
Papillary	48	35	0.19702	.6571
Non‐Papillary	54	47
Grade
Low	12	7	0.20189	.6532
High	91	75
T stage
T1‐2	62	48	0.065705	.7977
T3‐4	40	35
TNM stage
I‐II	56	39	0.4756	.4904
III‐IV	47	42
Survival status
Alive	80	68	0.28411	.594
Dead	23	15

Table 2

Univariate Cox proportional hazard regression analysis and LASSO of 7 genes

Gene	Univariate cox regression analysis			LASSO coefficient
	HR (95%CI)	95% CI	P value
ACTC1	2.766	1.145‐6.679	.024	0.025202
EPHA2	2.388	1.042‐5.477	.040	0.118639
OSBPL10	2.442	1.012‐5.891	.047	0.271164
FAM129A	2.560	1.087‐6.025	.031	0.114802
EGR1	3.226	1.332‐7.814	.010
PALLD	2.602	1.074‐6.301	.034

Figure 6

Construction of the 4‐gene‐based prognostic model. A, LASSO coefficient profiles of the 4 survival‐related RNAs. B, BC patients were divided into high‐risk and low‐risk groups. C, The expression heatmap of the 4 genes between high‐risk and low‐risk group. D‐E, Kaplan‐Meier survival curve of OS and DFS in the whole set

Validation of survival prediction in the integrated test sets

To validate the ability of this prognostic model to predict survival of BC patients, we draw Kaplan‐Meier survival curves of OS and DFS in the training and integrated test sets. In the training set, compared to the low‐risk group, the high‐risk group was tended to lower OS survival rate (P < .001) and DFS survival rate (P = .045) (Figure 7A,B). ROC curve was conducted to analyse the false‐positive rate and true‐positive rate of the 4‐gene‐based prognostic model, and the area under ROC curve based on OS at 1,3,5 years was 0.665, 0.734 and 0.735 (Figure 7C). In addition, we did the same analyses in the internal testing set (Figure 7D,E,F). By the validation, we showed the high accuracy prediction ability of this prognostic model.

Figure 7

Validation of survival prediction. A‐B, Kaplan‐Meier survival curve of OS and DFS in the training set. C, Time‐dependent ROC curves in the training set. D‐E, Kaplan‐Meier survival curve of OS and DFS in the integrated test set. F, Time‐dependent ROC curves in the integrated test set

Correlation between 4 gene classifier and clinicopathologic characteristics in the training set

Based on TCGA data, we obtained the clinical information of 186 BC patients. Chi‐square test was conducted on 6 clinical factors between the high‐risk and low‐risk group. As shown in Table 3, the 2 clinical characteristics (T stage and TNM stage) showed significant differences between the two groups. The patients in the high‐risk group tended to have higher staging.

Table 3

Correlations between risk score of the 4‐gene‐based classifier with OS and clinicopathological characteristics in the training set

Parameters	High risk	Low risk	χ 2	P Value
Age
≤65	6	35	1.6472	.1993
>65	17	45
Gender
Male	14	66	3.6526	.05598
Female	9	14
Subtype
Papillary	7	41	1.8934	.1688
Non‐Papillary	15	39
Grade
Low	1	11	0.2877169	.2911
High	22	69
T stage
T1‐2	9	53	4.7274	.02969
T3‐4	14	26
TNM stage
I‐II	7	49	5.6519	.01744
III‐IV	16	31

Construction of ceRNA network

First, the relationship among the 4 mRNAs and the 14 miRNAs which might target the 4 mRNAs were displayed. Among the network, miR‐124 could target OSBPL10 and EPHA2, while EPHA2 was simultaneously regulated by miR‐124, miR‐141, miR‐200a, miR‐26a, miR‐26b, miR‐302b, miR‐506, miR‐520d and miR‐95(Figure 8A). Then, we showed the correlation of mRNAs and miRNAs in Figure 8B‐C. By means of edgeR, the differentially expressed lncRNAs were analysed based on lncRNA expression data from TCGA. A total of 788 differentially expressed lncRNAs (249 up‐regulated and 539 down‐regulated) were identified (Figure 8D). Next, 48 overlapped lncRNAs were obtained between the 788 differentially expressed lncRNAs and the 649 lncRNAs predicted from 14 miRNAs. Finally, Figure 8E shows that we constructed a lncRNA‐miRNA‐mRNA ceRNA network including 4 mRNAs, 14 miRNAs and 48 lncRNAs.

Figure 8

Construction of ceRNA network. A, The linkage among the 4 mRNAs and the 14 miRNAs. B‐C, The correlation of the mRNAs and miRNAs. D, Volcano map of DElncRNAs. E, The construction of a ceRNA network including 4 mRNAs, 14 miRNAs and 48 lncRNAs. The purple spots represent the lncRNAs. The light green spots represent the miRNAs. The pink spots represent the mRNAs

DISCUSSION

BC, one of the leading causes of cancer death worldwide, is rarely diagnosed in individuals aged < 40 years. But the patient remains at risk for developing a new lesion probably with a more advanced stage when a tumour is diagnosed anywhere within the urothelial tract. Not surprisingly, after diagnosis of BC, most patients were reported significant decline in health‐related quality of life scores over time. The histological grade of malignant tumours shows the tumour development and is often used to evaluate tumour differentiation and patients’ prognosis, as well as bladder cancer. But, given the difficulty in obtaining the histological samples of BC patients, we want to explore some novel markers to screen BC earlier and predict the prognosis of patients with BC.

In recent years, increasing ceRNA networks, including mRNAs and lncRNAs, are involved in the process of diseases, especially human tumours. Molecular interaction in ceRNA network could promote co‐ordination of BP, and if the balance is compromised, it can lead to cancer. As a miRNA sponge, lncRNA HULC may be related to tumorigenesis. Once HULC is up‐regulated, it absorbs miR‐372 to the putative binding site and PRKACB is released from the binding of the miR‐372, while PRKACB could phosphorylate CREB, and thus controls HULC transcription 5. What's more, lncRNA‐miRNA‐mRNA ceRNA network might play a role for prediction prognosis of patients. A study, based on TCGA data, constructed a ceRNA network including 148 mRNAs, 60 miRNAs and 37 lncRNAs. Then by means of multivariate logistic regression, they established a 5‐lncRNA‐based survival model to predict the prognosis of patients with colon adenocarcinoma. In the hepatocellular carcinoma (HCC), BCYRN1 might regulate some cancer‐related pathways to promote HCC initiation via an lncRNA‐miRNA‐mRNA network. Through combing plasma BCYRN1 with alpha‐fetoprotein, the diagnosis of HCC was remarkably improved. Therefore, the study of RNA interactions at the molecular level will lead to a better understanding of the gene regulatory network in cancer.

In the present study, EdgeR was conducted to identify the DEmRNAs between the BC and normal bladder samples. Further study, we obtained 86 overlapped mRNAs among the predicted mRNAs, the DEmRNAs and the mRNAs in the yellow and lightcyan modules. And GO analysis showed that the 86 mRNAs were involved in the amoeboid cell migration that was accomplished by extension and retraction of a pseudopodium. Considering that some tumour cell lines have shown inherent amoeboid appearance, the ameboidal‐type cell migration in which 86 mRNAs participate may be the true mechanism of bladder cancer cell migration. Then, a survival model with 4 genes was constructed by univariate Cox regression and LASSO regression analysis: ACTC1 + FAM129A + OSBPL10 + EPHA2. Based on this model, we predicted risk score of BC patients and identified the high accuracy prediction ability of this prognostic model. At last, 4 mRNAs, 14 miRNAs and 48 lncRNAs were used to construct a ceRNA network. In this ceRNA network, we found that the patients in the high‐risk group tended to have higher staging. Molecularly, mRNA expression is regulated by miRNA and lncRNA. Based on this, miRNAs and lncRNAs might regulate mRNAs and thus affect the biological process and staging of BC.

EphA2, a member of the family of Ephrin receptor tyrosine kinases, is a functional signalling receptor for progranulin. Interaction of progranulin with EphA2 may act in cancer progression and tumour angiogenesis. Angiogenesis plays an important role in BC. A study found the immunosuppressive drug leflunomide inhibited angiogenesis in bladder carcinogenesis model via significant inhibition of the sEphrin‐A1/EphA2 system and might be a potential biomarker for treating BC beyond immunosuppressive therapy. In addition, EphA2 had been demonstrated to be essential in ceRNA network. In gastric cancer, miR‐302b exhibited anti‐tumour activity by reversing EphA2 regulation. And this modulation of EphA2 also had distinct effects on cell proliferation and migration in GC. Besides, ephA2/ephirnA1 signalling regulates miR‐302b expression and attenuates malignant pleural mesothelioma cell growth. Further study had shown that MIAT/miR‐520d‐3p/ EphA2 might be a new target for HCC therapy. Other studies showed the interaction between miR‐520d and EphA2 could enhance tumour suppression in ovarian cancer and gastric cancer. ,

MiR‐26a and miR‐26b play important roles in BC as tumour suppressors. Overexpression of miR‐26a/miR‐26b inhibited the proliferation and migration of BC cells. , Importantly, a study found that blood miR‐26b detected the presence of invasive bladder tumours with 94% specificity and 65% sensitivity. And with increasing tumour‐nodes‐metastasis staging in bladder cancer, miR‐26b showed a moderate decreasing trend. As for interaction between EphA2 and miR‐26b in cancers, EphA2 was verified as the target of miR‐26b. MiR‐26b could enhance radiosensitivity of hepatocellular carcinoma cells by targeting EphA2. , And overexpression of miR‐26b dramatically inhibited the proliferation, invasion and migration of HCC cells by targeting EphA2. So, among the network we constructed, we wonder if the interaction between EphA2 and miR‐26b could play roles in tumour angiogenesis and proliferation in BC. We will do further research to validation our hypothesis.

Besides, FAM129A/miR‐195, ACTC1, XIST/ miR‐124/miR‐133a, miR‐506, miR‐200a and miR‐141 were all related to the development of BC. , , , , , Thus, a lot of previous studies have demonstrated the molecular functions of the ceRNA network in the occurrence and development of BC. Among the ceRNA network we constructed, we took advantage of the 4 genes to explore the prognosis of BC patients which indicated the higher risk group was tended to a poor prognosis. And based on the 4 genes, we constructed a ceRNA network as a biomarker for BC patients. However, we do not know how mRNAs, miRNAs and lncRNAs interact in this network. For example, we predicted that lncRNAs could bind miRNAs through MRE and indirectly regulate mRNAs expression. However, as for the specific signalling pathway in which this regulatory network influences specific biological processes, we still do not know which needs further research.

Inevitably, our research had some innate limitations which need to be addressed. Firstly, although this bioinformatics study was designed well, and a strict threshold for WGCNA analysis was set up, the major drawback of our study was the lack of in vivo and vitro validation. Secondly, the current study was of a retrospective nature, as it was based on data from TCGA dataset without validating it in a prospective clinical trial. Despite these drawbacks, the results of our study demonstrate that our survival model with 4 genes could be used as reliable prognostic predictors of BC, as well as the ceRNA network plays a role in the molecular regulatory network of BC.

CONFLICTS OF INTEREST

The authors declare that there are no conflicts of interest.

AUTHOR CONTRIBUTIONS

S. L. conceived and designed the study, J. J., Y. B. and X. L. performed the analysis procedures, D. Y., X. Y., J. Y, T. L. and J. J. analysed the results, T. L., and S. L. contributed analysis tools, J. J., and S. L. contributed to the writing of the manuscript. All authors reviewed the manuscript.