Competing endogenous RNA network analysis identifies critical genes among the different breast cancer subtypes

Although competing endogenous RNAs (ceRNAs) have been implicated in many solid tumors, their roles in breast cancer subtypes are not well understood. We therefore generated a ceRNA network for each subtype based on the significance of both, positive co-expression and the shared miRNAs, in the corresponding subtype miRNA dys-regulatory network, which was constructed based on negative regulations between differentially expressed miRNAs and targets. All four subtype ceRNA networks exhibited scale-free architecture and showed that the common ceRNAs were at the core of the networks. Furthermore, the common ceRNA hubs had greater connectivity than the subtype-specific hubs. Functional analysis of the common subtype ceRNA hubs highlighted factors involved in proliferation, MAPK signaling pathways and tube morphogenesis. Subtype-specific ceRNA hubs highlighted unique subtype-specific pathways, like the estrogen response and inflammatory pathways in the luminal subtypes or the factors involved in the coagulation process that participates in the basal-like subtype. Ultimately, we identified 29 critical subtype-specific ceRNA hubs that characterized the different breast cancer subtypes. Our study thus provides new insight into the common and specific subtype ceRNA interactions that define the different categories of breast cancer and enhances our understanding of the pathology underlying the different breast cancer subtypes, which can have prognostic and therapeutic implications in the future.


INTRODUCTION
Non-coding RNAs are a class of RNAs that regulate gene expression transcriptionally and posttranscriptionally [1]. The miRNAs are a class of ~22-nucleotide long single-stranded non-coding RNAs that regulate gene expression by binding to MREs (miRNA response elements) on the RNAs [2]. Aberrant miRNA expression patterns have been associated with many human cancers including breast cancer [3]. The lncRNAs (long non-coding RNAs) are another class of non-coding RNAs (longer than 200 nucleotides), that are also involved in the pathology of many complex human diseases including cancers [4]. The lncRNAs also harbor MREs and compete with other RNAs for miRNA binding, thus can function as competing endogenous RNAs (ceRNAs) to influence the post-transcriptional regulation in cancers [5,6].
Breast cancer is a progressive disease that is a leading cause of cancer-related mortality in women and composed of distinct subtypes that are highly heterogeneous [7][8][9][10]. Based on gene expression, breast cancer is classified into five major subtypes namely luminal A, luminal B, HER2-enriched, basal-like, and normal-like [11,12]. The luminal subtypes (A and B) are characterized by the expression of ER (estrogen receptor) and/or PR (progesterone receptor) related genes [13]. The HER2-enriched subtype is often associated with frequent HER2/ERBB2 amplification (80%) and therefore HER2/ ERBB2 is targeted for therapy [14]. The basal-like subtype is generally negative for ER, PR and HER2 and therefore also called as triple-negative breast cancer or TNBC, www.impactjournals.com/oncotarget/ Oncotarget, 2017, Vol. 8, (No. 6), pp: 10171-10184 that is characterized by enhanced invasiveness and poor clinical outcomes [15]. The etiology of different breast cancer subtypes shows variant risk factors that may result in different clinical outcomes [16,17]. Therefore, understanding the similarities and differences among the various subtypes is important from the clinical standpoint.
Although many breast cancer biomarkers have been reported [18,19], it has come to light that many of these biomarkers are coordinately regulated by similar miRNAs that implies competition among shared miRNAs to bind to different RNAs. For example, Li et al. found that the STARD13-correlated ceRNA network inhibits epithelialmesenchymal transition (EMT) and metastasis of breast cancer [20]. However, these studies focused on the clinical utility of individual ceRNA interactions and did not investigate the ceRNA interactions at a system level. Xu et al. conducted a systemic analysis of ceRNA interactions in breast cancer without considering the heterogeneity of breast cancer subtypes and the involvement of lncRNAs [21]. Also, studies investigating ceRNA interactions in human cancer have focused on dysregulated RNAs that are aberrantly expressed during cancer initiation and progression [22]. Due to the heterogeneity of breast cancers, the expression profiling is more complex [7,23]. Hence, a systemic analysis of ceRNA crosstalk among different breast cancer subtypes may yield better information regarding the interplay of various biological networks that are involved in breast carcinogenesis [24]. Furthermore, since the potential interactions between the subtype molecular biomarkers have not been considered in previous studies, we opined that the ceRNA network analysis would help dissect the pathway interactions at a molecular level.
Therefore, the aim of this study was to generate subtype specific dysregulated miRNA mediated ceRNA networks to discover new critical subtype related genes and to unravel the shared and specific pathogenic pathways of the breast cancer subtypes at a system level

The miRNA dys-regulatory networks in the breast cancer subtypes
To systematically identify disease markers in each of the four breast cancer subtypes, we performed Wilcoxon Rank Sum tests to identify differentially expressed protein coding and non-coding genes in each of the breast cancer subtypes ( Figure 1A, Supplementary Table 1). Each subtype exhibited about 40% differentially expressed coding and non-coding genes ( Figure 1A). Further, majority of the miRNAs were upregulated and most of the lncRNAs were downregulated in all the subtypes. About 70% of dys-regulated coding and non-coding RNAs are shared in at least two breast cancer subtypes ( Figure 1B). Also, about 95% of subtype shared dysregulated mRNAs, miRNAs and lncRNAs are concordantly changed across breast cancer subtypes, comparing the four breast cancer subtypes, namely luminal A, luminal B, HER2-enriched and basal-like ( Figure 1B). These data indicated that these shared dys-regulated RNAs may result in accordant direction change of cancer-related pathways in different breast cancer subtypes. Among the subtypes, basal-like exhibited the highest fraction of dysregulated coding and non-coding RNAs that lead to malignancy as shown in supplementary Table 2 [15,25]. Notably, several subtype marker genes were identified that were consistent with previous studies like ESR1, PGR (progesterone receptor) and FOXC1 [11]. For example, ER, ESR1 and PGR genes were downregulated in the basal-like subtype, but, upregulated in the luminal A/B subtypes. Similarly, FOXC1 was upregulated in the basal-like subtype and down-regulated in the luminal A/B and HER2-enriched subtypes.
Regarding miRNAs, previous studies have suggested that miRNAs usually function as negative regulators of gene expression [2,26]. Therefore, we constructed the subtype miRNA dys-regulatory networks by assuming that the miRNAs negatively regulated their target genes (Supplementary Table 3). Our data suggested that most of the miRNA targets were subtype specific. However, the subtypes with similar tissue-of-origin, such as the luminal A and the luminal B shared many nodes and edges in the networks ( Figure 1C-1E). Also, the network structure of the miRNA dys-regulatory networks was organized according to a core set of principles and not randomly. This suggested that their analysis would help understand the regulation between the miRNAs and their targets at a system level (Supplementary Figure 1A-1B).

CeRNA networks among the breast cancer subtypes
To investigate the ceRNA regulation in breast cancer subtypes, we generated a ceRNA network for each subtype by assembling all the significant ceRNA interactions that were based on the significance of both, their positive co-expression and the shared miRNAs, in the corresponding subtype miRNA dys-regulatory network (Table 1). Accordingly, we identified 84,026, 11,449, 4,982 and 8,600 ceRNA interactions in the networks of luminal A, luminal B, HER2-enriched and basal-like subtype, respectively. We found 26.88%, 11.79%, 8.07% and 12.54% of dysregulated mRNAs as well as 6.50%, 2.58%, 1.75% and 3.30% of dysregulated lncRNAs in the ceRNA networks of luminal A, luminal B, HER2enriched and basal-like subtype, respectively. This indicated that more dys-regulated mRNAs are influenced by the ceRNA regulation mechanism compared with lncRNAs. Further, we observed that many ceRNAs were well connected and formed a tight network in each subtype (Supplementary Figure 2A-2D). Similar to the miRNA dys-regulatory networks, the ceRNA networks showed that the ceRNA interactions are much more specific to each subtype in spite of overlapping nodes among the four different subtypes, suggesting unique subtype specific regulation (Figure 2A-2B). As expected, the luminal A and B subtypes that are similar in the tissue of origin and at the molecular level showed greater overlap of the ceRNA network (Figure 2A-2B). We also observed that the degree distribution of the ceRNA networks satisfied the power law model suggesting that the networks were scale-free and not randomly organized implying rapid communication among different ceRNAs (Supplementary Figure 2E).
Several marker genes of specific breast cancer subtypes have been reported to be differentially expressed [9,11]. We observed that these marker genes were highly connected in the ceRNA networks. For example, in luminal A breast cancer, EGFR and FOXC1 had 453 and 343 ceRNA connections, respectively and shared four miRNAs (miR-15b-3p, miR-141-3p, miR-200a-3p and miR-760). This suggested that the marker genes of breast cancer subtypes can co-regulate each other and have high connectivity in the ceRNA networks.

Subtype common ceRNAs tend to be at the core of ceRNA networks
A key network property is node connectivity that shows the degree of communication among the nodes.
Hub node is pivotal and shows high connectivity with other nodes. In our study, we defined the hubs as the top 10% of the nodes with the highest degree within the networks. The hub nodes were sub-divided into subtype common and subtype specific ceRNA nodes that would identify mechanisms that are either shared or specific among the different breast cancer subtypes. The subtype common ceRNA hubs were defined as ceRNA hubs that occurred in at least three of the subtype ceRNA networks and accounted for 12.65%, 29.07%, 39.34% and 19.39% of total ceRNA hubs in luminal A, luminal B, HER2enriched and basal-like, respectively ( Figure 3A). On the other hand, the subtype specific ceRNA hubs were defined as hubs that occurred in each of the specific subtype networks only. Also, the luminal A and basal-like had high proportion of subtype specific ceRNA hubs relative to common ceRNA hubs.
While the node degree of a network represents the local centrality, k-core decomposition reveals the modularity of central hub nodes and is used to identify the largest sub-graph where the vertices have at least 'k' connections within a subgraph called as the k-cores [27]. Therefore, hubs that are closely connected with other hub nodes are within k-cores have larger k values, whereas, hubs closely connected to non-hub nodes have smaller k values. When we conducted k-core decomposition analysis to analyze the localization of common and subtype specific nodes, we observed that the subtype common ceRNAs had larger k values compared to subtype-specific ceRNA nodes ( Figure 3B). This indicated greater centrality in the networks and suggested that all breast cancer subtypes shared pathogenic pathways. Also, the Wilcoxon Rank Sum test revealed that the subtype common ceRNA hubs had higher connectivity compared with subtype-specific ceRNA hubs in each network, which may provide some clues of common biological processes disturbed among the different breast cancer subtypes ( Figure 3C).

Subtype common ceRNA hubs are associated with proliferation in breast cancers
Next, we analyzed the relationship among common ceRNA hubs to understand their regulation of breast cancer progression. Towards this, we collected all the ceRNA interactions among the 28 common ceRNA hubs in the four breast cancer subtypes to form a subtype common sub-network ( Figure 4A-4B). Five of the twenty eight ceRNA hubs were common to all the four breast cancer subtype networks, including RECK, which controls breast cancer metastasis by regulating MMP-2 (matrix metalloproteinases-2) to inhibit tumor angiogenesis [28]. We also found that these ceRNA hubs were highly connected with each other and regulated by subtype common hub miRNAs in the corresponding miRNA dys-regulatory network ( Figure 4A). Furtherly, 147 of the 353 ceRNA interactions among the 28 subtype common ceRNA hubs occurred in at least three breast cancer subtypes. The subtype common ceRNA network was regulated by 13 common hub miRNAs in at least three miRNA dys-regulatory networks, including breast cancer related miR-141-3p, miR-33a-5p, miR-200c-3p, and miR-182-5p [29][30][31]. We also identified 25 protein coding genes among the 28 subtype common ceRNA hubs, of which 12 (RECK, LHFP, FSTL1, TGFBR2, EBF1, ERG, MRVI1, DKK3, GAS7, RBMS3, PDGFRA and CCDC80) were listed as cancer related genes in the Cancer Gene Census (CGC) of the COSMIC database or the Genetic Association Database (GAD) [32,33]. LHFP, a common ceRNA hub in all the four breast cancer subtypes is associated with mesenchymal differentiation in glioma [34]. Defective MSRB3 has been shown to inhibit cell proliferation through the activation of p53-p21 and p27 pathways [35]. Also, MSRB3 and LHFP coregulate each other by sharing common hub miRNAs, miR-141-3p and miR-33a-5p (    Furthermore, our data showed that the lncRNAs functioned as ceRNAs involved in the cancer related biological processes. We found RP11-276H19.1 (ENSG00000226237) as a common hub that participated in the ceRNA interaction networks in luminal A, luminal B and HER2-enriched subtypes. It also co-regulated with MSRB3 and LHFP by sharing miRNAs to form ceRNA-ceRNA interactions in the three subtypes ( Figure  4C-4D). Also, the expression of these three ceRNAs had similar trend in each breast cancer subtype ( Figure  4D). Functional enrichment analysis of the ceRNA hubs revealed their involvement in proliferation (FDR  4.80e-03, hypergeometric test, Bonferroni corrected), MAPK signaling pathway (FDR = 4.10e-03), tube morphogenesis (FDR = 4.37e-03) and some breast cancer associated processes (Table 2, Supplementary Table 4), demonstrating that the proliferation process was altered substantially in all breast cancer subtypes [36].

Subtype specific ceRNA hubs contribute to breast cancer subtype phenotype
To discover the pathogenic mechanisms of different breast cancer subtypes, we further focused on the subtype specific ceRNA hubs. We extracted four subtype-specific ceRNA hub sub-networks by collecting all subtypespecific ceRNA hubs and their ceRNA interactions within each of the four subtype ceRNA networks and found that the subtype-specific hubs were highly connected with each other and regulated by hub miRNAs in the miRNA dys-regulatory network of each subtype (Figure 5A-5D; Table 3; Supplementary Table 5). Several miRNA hubs were subtype specific and had specific roles in the corresponding breast cancer subtype. For instance, the three members of the miR-17-92 cluster (miR-17-3p, miR-20a-5p and miR-19b-3p) were part of an oncogenic cluster and associated with lymph node metastases in TNBC [37]. These three miRNAs were significantly upregulated in basal-like subtype compared with normal breast tissue. In addition, KDR or VEGF2, is also a basal-like subtype specific hub and a specific biomarker in TNBC based on a long-term follow-up dataset [38]. KDR acted as a ceRNA partner of PGR, another basallike specific hub and competing for the binding to the miR-17-92 cluster (including miR-17-3p, miR-19b-1-5p and miR-20a-5p). Interestingly, we found that these subtype-specific ceRNA hubs participated in breast or other hormone related cancer like prostate cancer, ovarian cancer and endometrial cancer ( Table 3, Supplementary  Table 5), consistent with previous findings that hormone related cancers had similar pathways [39]. In addition, some biological pathways were identified to be assigned to specific subtypes. For example, estradiol response and IL2 -STAT5 signaling (inflammatory response pathway) were luminal A and luminal B specific (Table 3), whereas, the coagulation pathway was specifically involved in the basal-like subtype (FDR = 1.16e-02).
Based on these data, we hypothesized that the subtype-specific ceRNA hubs that took part in subtype specific pathways could be used to classify different breast cancer subtypes. We also performed the principal component analysis (ggfortify R package) and found that the 29 subtype-specific ceRNA hubs that were dysregulated among all breast cancer subtypes (Wilcoxon Rank Sum test, FDR<0.05) could distinguish the basallike, HER2-enriched and luminal subtypes, whereas distinguishing the luminal A and B subtypes was not clearcut ( Figure 5E). Taken together, these results suggested that investigating the properties of ceRNA interaction networks may help strengthen our understanding of different subtype breast cancers, which may contribute to breast cancer subtype phenotype, thus further to improve clinical prevention and treatment.

DISCUSSION
In this study, we developed a computational framework to construct and analyze ceRNA networks using the transcriptional profiling data to understand the dynamic clinical behaviors and outcomes of different breast cancer subtypes.
In this study, we adopted a three-step design to construct the ceRNA interaction network in each subtype to eliminate false discoveries [40,41]. Thereafter, we used highly stringent threshold values to construct miRNA dys-regulatory networks and the identified common and subtype-specific ceRNA hubs were not influenced by different threshold effects (Supplementary Table 6). In addition, the maximal FDRs for luminal A, luminal B, HER2-enriched and basal-like subtypes were 0, 2.06e -12 , 1.18e -04 , 7.11e -09 , respectively, suggesting that the ceRNA interactions were reliable. We further validated the accuracy of predicting the ceRNA interactions globally by analyzing expression datasets for the four breast cancer subtypes from another array, GSE45827, and found that the ceRNA interaction pairs showed higher positive co-expression compared with the gene pairs that were randomly selected (Supplementary Figure 3).
Further, although majority of the ceRNA nodes and interactions were subtype specific, several common ceRNAs were also shared by the four breast cancer subtypes (Figure 2A and 2B). These common ceRNA hubs indicate shared dynamic changes in many biological processes like MAPK signaling pathway, proliferation process and tube morphogenesis during breast tumorigenesis (Table 2). Conversely, the subtype specific ceRNA hubs provide clues to specific changes during the transformation of the normal cells into the different breast cancer subtypes (Table 3). For example, our analysis showed that luminal A and luminal B were associated with estradiol response and inflammatory response pathway. This validated our analysis since the luminal subtypes are recognized as estrogen receptor positive breast cancer (ER+) and are clinically treated with hormone therapy [13]. Also, interleukin-6 that is secreted by the transformed fibroblasts induces tamoxifen resistance in luminal breast cancer [42]. On the other hand, we found that the hallmark of the basal-like breast cancer subtype was coagulation pathway, consistent with previous reports that basal breast cancer cells constitutively express coagulation initiator, an alternatively spliced tissue factor (asTF) www.impactjournals.com/oncotarget   that promotes breast cancer growth in a β1 integrindependent manner [43]. In regard to the tumor suppressive and the oncogenes reported [42], our study demonstrated that both the common and the subtype-specific ceRNA hubs were enriched in tumor suppressive genes and depleted in oncogenes (Supplementary Table 7). To explain this, it has to be noted that the ceRNA hubs are down-regulated in each of the breast cancer subtypes. Therefore, the dysfunction of the enriched tumor suppressor genes probably drives the initiation and progression of breast cancer. Several cancer associated lncRNAs including the common ceRNA hubs like HAND2-AS1 (DEIN) [44,45] and RP11-276H19.1 [46], as well as subtypespecific hubs like ADAMTS9-AS2, a luminal A specific hub lncRNA, acted as tumor suppressors and inhibited the migration of the glioma cells [47]. Similarly, functional analysis of the miRNA hubs by the TAM tool showed that the common and specific miRNA hubs play important roles in breast cancer subtypes (Supplementary Table 8) [48].
Finally, our analysis revealed 29 critical ceRNA hub genes that could distinguish the basal-like, HER2enriched and luminal breast cancer subtypes. The Luminal A/B subtypes were similar as they shared many expression markers within the luminal epithelial cell layers. HER2-enriched and luminal A/B subtypes could not be well distinguished as several HER2-enriched subtype showed the expression of ER or PR as well as a global gene expression profile that was similar to the luminal A and B subtypes [49,50]. Among the 29 subtype related genes, four were part of the 50-gene PAM50 model (FOXC1, SFRP1, PGR and EGFR) [12]. Although the other 25 genes were not part of the PAM50 model, they probably compete for miRNA binding and thus influence the subtype specificities. For example, DCN (decorin) regulates PGR by competing for miR-19b-1-5p and promotes the basal-like subtype (Pearson correlation coefficient =0.54). This supports previous findings that DCN is down-regulated and influences the breast cancer cell motility and invasion, and acts as one of the six-gene signature that adds prognostic value independent of the expression of ER, PR, and HER2 [51,52]. Therefore, subtype specific expression of 50-gene in PAM50 model can be regulated by other ceRNAs that may influence breast cancer subtype specificity. Interestingly, three of the 29 subtype-specific ceRNA hubs, namely EGFR, IL6ST and MET, are listed in the Cancer Gene Census of the COSMIC database [32]. Therefore, all these data put together identify these 29 genes as crucial for breast cancer subtype determination, maintenance and function and need to be further analyzed in detail.
Systematic studies have shown that somatic mutations are the driving force for human carcinogenesis [53]. In our study, among the 28 common ceRNA hubs, we found 18 mutant genes in at least one breast cancer subtype. Therefore, when we analyzed gene mutations in the general ceRNA networks, we observed that the mutated ceRNAs are highly connected than the other ceRNAs within the four ceRNA networks (Supplementary Figure 4A, Wilcoxon Rank Rum test). This also showed that mutated ceRNAs were enriched in the hubs across all the four breast cancer subtypes (Supplementary Figure 4B, hypergeometric test). Moreover, k-core subnetwork analysis indicated that the mutated ceRNAs were localized in the highly dense layers of each ceRNA network, especially in the HER2-enriched and the basallike subtypes (Supplementary Figure 4C). Although the effects of somatic mutation in human cells is largely unknown [54], since mutated genes had greater ceRNA interactions, we postulate that somatic mutation may drive breast cancer initiation and progression by dysregulating the normal ceRNA interactions.
Therefore, in conclusion, our study provides a framework to compare disease subtypes based on miRNA mediated ceRNA interactions. The systemic network analyses provide understanding of the biological mechanisms underlying the different breast cancer subtypes and help identify new critical genes for each of the breast cancer subtypes. We predict that systematic characterization of ceRNA-ceRNA interactions will contribute to improving the subtype classification and also identifying novel prognostic markers that would help selecting appropriate subtype specific therapies that would improve the survival outcomes for breast cancer patients.

Genome wide transcriptome and mutation datasets of the five breast cancer subtypes
The mRNA and miRNA expression datasets measured by RNA sequencing were downloaded from The Cancer Genome Atlas (TCGA) [23], whereas, the lncRNA sequencing datasets were obtained from The Atlas of ncRNA in Cancer (TANRIC) [55]. The molecular subtype information of breast cancer was downloaded from Synapse (https://www.synapse.org/, syn1461151) for the luminal A, luminal B, HER2enriched, basal-like and normal-like, with normal breast samples [23]. As the normal-like subtype has high percentage of normal sample "contamination" in the tumor specimen [12], it was excluded from further analysis (Supplementary Table 1). The mutation datasets of each tumor sample were also obtained from Synapse. To filter out lowly expressed RNA transcripts (mRNA, miRNA, lncRNA), the one with expression value of 0 in all of the samples was removed. Finally, 20,211 mRNAs, 12,697 lncRNAs and 1,026 miRNAs were considered for further analysis. All the expression profiles were log2 transformed.

The miRNA-target relationships
To obtain the regulatory relationships between miRNAs and mRNAs (or lncRNA), miRanda algorithm was used to predict the miRNA-target relationships. As a result, 1,129,584 unique miRNA-lncRNA and 2,927,669 unique miRNA-mRNA relationships were generated among 1,017 miRNAs, 7,731 lncRNAs and 17,215 mRNAs, respectively.

Functional annotation datasets
Functional annotation gene sets were downloaded from the Molecular Signatures Database (MsigDB) [56]. This included 50 hallmark gene sets, 825 Biological Process (BP) terms for Gene Ontology (GO), 3396 chemical and genetic perturbations (CGP) and 186 Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets.

Overview of the construction of subtype ceRNA networks
A three-step model was proposed to construct ceRNA networks. Firstly, we identified differentially expressed mRNA, miRNAs and lncRNAs for each of the breast cancer subtypes compared with normal breast samples. Secondly, four miRNA dys-regulatory networks were constructed by considering negative regulation between dysregulated miRNAs and their targets. Thirdly, based on the miRNA dysregulatory networks, ceRNA interactions were identified by considering significant sharing of the the regulatory miRNAs and their co-expression. All ceRNA interactions were assembled into a ceRNA network for each breast cancer subtype, where the nodes were referred to as ceRNAs and the edges referred to as ceRNA interactions.

Identification of differentially expressed mRNAs, miRNAs and lncRNAs
Wilcoxon Rank Sum tests were used to identify differentially expressed protein coding genes and noncoding genes between each subtype and normal breast samples. We applied the Bonferroni correction and the RNAs that showed a FDR less than 0.05 and the absolute log2 fold change (LFC) no less than 1 were considered significant.

Construction of subtype miRNA dys-regulatory networks
Differentially expressed target genes (mRNAs and lncRNAs) for each differentially expressed miRNA were collected and a miRNA dys-regulatory network was constructed for each subtype by considering negative regulation between miRNAs and their targets based on the expression profiles of each breast cancer subtype samples (Pearson Correlation Coefficient<0, P<0.01).

Construction of subtype ceRNA networks
The subtype ceRNA interaction networks were constructed based on the miRNA dys-regulatory networks. Two target genes were considered to be competing only if 1) they significantly shared regulatory miRNAs in the miRNA dys-regulatory network as calculated by the hypergeometric test (FDR<0.01, Bonferroni corrected) and 2) the expression of ceRNA-ceRNA interactions showed a statistically significant positive correlation with a Pearson correlation coefficient that was no less than 0.5 (FDR<0.01, Bonferroni corrected). All the ceRNA interactions were assembled into a ceRNA network for each breast cancer subtype and the weights of all edges were set to 1.

K-core decomposition analysis of the ceRNA networks
The k-core is the largest sub-graph where every node has at least k links. The k-core can be determined by pruning all the vertices from a network with degrees less than k and their incident links until the degree of all the remaining vertices in the sub-graph are larger than or equal to k. The k layers reflect the connectivity of the ceRNA partners such that the first peeled layer (k=1) contains ceRNAs with only a few connections whereas the last peeled layer is most connected. Each node was assigned to a maximal k core sub-network in each breast cancer subtype ceRNA network. For among different breast cancer subtypes, the layers of each ceRNA network were normalized to 0-1 divided by the number of largest k layer in each network.

Functional enrichment analysis for subtype common and specific ceRNA hub sub-networks
To dissect the biological pathways represented by the common or subtype-specific ceRNA hub sub-networks, functional enrichment analysis of the hub ceRNA genes was performed by hypergeometric distribution with the Bonferroni correction (FDR < 0.05).