Reconstructing the coding and non-coding RNA regulatory networks of miRNAs and mRNAs in breast cancer

Gene 548 (2014) 6–13

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Reconstructing the coding and non-coding RNA regulatory networks of miRNAs and mRNAs in breast cancer Sheng Yang a, Hui Zhang a, Li Guo a,b, Yang Zhao a, Feng Chen a,b,⁎ a b

Department of Epidemiology and Biostatistics, School of Public Health, Nanjing Medical University, Nanjing, China Ministry of Education Key Lab for Modern Toxicology, School of Public Health, Nanjing Medical University, Nanjing, China

a r t i c l e

i n f o

Article history: Received 19 March 2014 Received in revised form 27 May 2014 Accepted 5 June 2014 Available online 27 June 2014 Keywords: Integrated analysis miRNA–mRNA Breast cancer (BC)

a b s t r a c t microRNAs (miRNAs) are a class of small non-coding RNAs that deregulate and/or decrease the expression of target messenger RNAs (mRNAs), which specifically contribute to complex diseases. In our study, we reanalyzed an integrated data to promote classification performance by rebuilding miRNA–mRNA modules, in which a group of deregulated miRNAs cooperatively regulated a group of significant mRNAs. In five-fold cross validation, the multiple processes flow considered the biological and statistical significant correlations. First, of statistical significant miRNAs, 6 were identified as core miRNAs. Second, in the 13 significant pathways enriched by gene set enrichment analysis (GSEA), 705 deregulated mRNAs were found. Based on the union of predicted sets and correlation sets, 6 modules were built. Finally, after verified by test sets, three indexes, including area under the ROC curve (AUC), Accuracy and Matthews correlation coefficients (MCCs), indicated only 4 modules (miR-106b-CITKPNA2-miR-93, miR-106b-POLQ-miR-93, miR-107-BTRC-UBR3-miR-16 and miR-200c-miR-16-EIF2B5-miR-15b) had discriminated ability and their classification performance were prior to that of the single molecules. By applying this flow to different subtypes, Module 1 was the consistent module across subtypes, but some different modules were still specific to each subtype. Taken together, this method gives new insight to building modules related to complex diseases and simultaneously can give a supplement to explain the mechanism of breast cancer (BC). © 2014 Elsevier B.V. All rights reserved.

1. Introduction MicroRNAs (miRNAs) are a series of endogenous, conserved and small (19–23 nts) non-coding RNAs (ncRNAs), which decrease the expression of their target messenger RNAs (mRNAs) at the posttranscriptional level via miRNA–mRNA association (Bartel, 2009; Kong et al., 2013). The maturation process of these small miRNAs generates two intermediates, including primary miRNA (pri-miRNA) and precursor miRNA (pre-miRNA). Initially, pri-miRNA is transcribed from DNA and then it is excised to pre-miRNA by RNA polymerase II, such as Drosha (Davis et al., 2010; Takeshita et al., 2013). Pre-miRNA is then transported from the nucleus to the cytoplasm, and Dicer, a key enzyme, cleaves it to form a miRNA:miRNA* duplex (Ota et al., 2013). miRNAs play versatile

Abbreviations: AUC, area under ROC curve; BC, breast cancer; BL, basal-like breast cancer; CI, confidential interval; FC, fold change; FDR, false discovery rate; GEO, Gene Expression Omnibus; GOBP, Gene Ontology Biological Process; GSEA, gene set enrichment analysis; HGT, hypergeometric distribution test; LA, luminal A breast cancer; LB, luminal B breast cancer; MCC, Matthews correlation coefficients; miRNA, microRNA; mRNA, messenger RNA; OR, odds ratio; ROC, receiver operating characteristic curves; SAM, significant analysis of microarray. ⁎ Corresponding author at: 818 East Tianyuan Road, Nanjing, Jiangsu 211166, China. E-mail address: [email protected] (F. Chen).

http://dx.doi.org/10.1016/j.gene.2014.06.010 0378-1119/© 2014 Elsevier B.V. All rights reserved.

roles in many biological processes, such as cell profiling, development and apoptosis, and also contribute to occurrence and development of complex diseases, especially cancers (Che and Huang, 2012; Hwang and Mendell, 2007; Lussier et al., 2012; Pritchard et al., 2012; Volinia et al., 2012). For example, both miR-335 and miR-31 have been identified as robust inhibitors in breast cancer (BC), and they further influence its metastasis of BC (Png et al., 2011; Tavazoie et al., 2008; Valastyan et al., 2009). Nevertheless, miRNAs lack the specificity in diagnosing certain diseases. Moreover, it is important to explore the specific associations of miRNA–mRNA from the many associations identified to date (Guo et al., 2013; Miska et al., 2007; Xu et al., 2011). An increasing number of research programs have proposed a series of approaches to identify the regulatory network or modules of miRNA–mRNA associations by integrating their expression profiles and/or the predicting datasets (Enright et al., 2003; Lewis et al., 2005; Liu et al., 2010; Xiao et al., 2013). For example, the network between miRNA and mRNA related to HCV infection has been constructed by a theoretical approach (Peng et al., 2009). Recently, Xiao et al. have detected dysfunctional miRNA– mRNA regulatory modules (dMiMRMs), involving the miRNA interactions and protein interactions, to discriminate the different degrees of glioblastoma and to diagnose this disease (Xiao et al., 2013). Thus, a method integrating the information of miRNAs and mRNA, that reduces

S. Yang et al. / Gene 548 (2014) 6–13

the Type I error in the selection of significant mRNAs and considering the interactions between mRNAs, is both urgent and necessary. Herein, we identified miRNA–mRNA modules from miRNA and mRNA microarray data. Firstly, significant deregulated miRNA expression profiles were screened based on hypergeometric distribution test (HGT) and fold change value (FC), which is similar to significant analysis of microarray (SAM). Secondly, gene set enrichment analysis (GSEA) of mRNAs was used to detect the significant pathways relative to BC. Thirdly, if the mRNA was statistically and biologically related to miRNAs, these associations were selected. The module involved mRNAs from the same pathway and deregulated miRNAs. Finally, Matthews correlation coefficient (MCC), area under the ROC curve (AUC) and Accuracy were the proxy of the discrimination ability of the modules.

7

2. Material and methods 2.1. Data resources Data was downloaded from Gene Expression Omnibus (GEO), including miRNA data (GSE19536) and mRNA data (GSE19783), and they were collected from the same cohort study in Norway. After excluding the unexpressed miRNAs and mRNAs, 489 miRNAs and 19593 mRNAs (up-regulated 12533 and down-regulated 7060) were identified in 94 samples (83 cases and 11 controls). Cases consisted of four diverse subtypes: 15 ERBB2, 12 luminal B (LB), 15 basal-like (BL) and 41 luminal A (LA) (Enerly et al., 2011). When one mRNA was detected by different probes, the average expression level was used.

Fig. 1. The workflow of the 3 steps for constructing special modules of breast cancer (BC). Step 1: identifying significant miRNAs, including non-core miRNAs and core miRNAs. Common miRNAs are selected by SAM (FDR ≤ 0.1 and log2FC ≥ 0.5). Also, binding to the information of their targets, core miRNAs are identified by hypergenomic distribution test. Step 2: identifying significant pathways. The pathways with values of count ≥ 50, P ≤ 0.05 and FDR ≤ 0.1 are defined as significant. Step 3: calculating Spearman correlation and constructing modules. The Spearman correlations between common dysregulated miRNAs are calculated. The mRNAs shared by significant pathways and correlated to dysregulated miRNAs are defined as significant mRNAs. The other basis of the modules is biological and significant mRNAs, which are common between the predicted and significant mRNAs. A significant and biological module is comprised of two related miRNAs and their co-targets (see the Material and methods section for details).

8

S. Yang et al. / Gene 548 (2014) 6–13

2.2. Module assembly To build significant modules between the two molecules, we followed three steps (Fig. 1). In step 1, the core dysregulated miRNAs were defined. Two algorithms were performed: FC and HGT (Enerly et al., 2011; Xiao et al., 2013). Common dysregulated miRNAs had a log-fold change of (log2FC) ≥ ± 0.5 and a false discovery rate of (FDR) ≤ 0.1 (Tusher et al., 2001). A permutation test was applied to obtain FDR, because the distribution of expression level was unclear. 10,000 replicates were shuffled between case and control to acquire a robust result. All common deregulated miRNAs were involved to construct the module. Next, the HGT was used to define core deregulated miRNAs by the predicted targets. Three databases were involved, including TargetScan, miRanda and miRTarBase (Griffiths-Jones et al., 2006; Hsu et al., 2011; Lewis et al., 2005). To improve precision, the selected targets were predicted by at least 2 algorithms. If the targets were out of the range of detected mRNAs, they were removed from the analysis. For a special miRNA, the P-value was calculated as follows:  m X

P ¼ 1−

M k

k¼0



N−M n−k   N n

 ð1Þ

where N represented the number of total genes in this microarray, n indicated the number of predicted targets, M was the inversely differentially expressed genes and m was the inversely deregulated targets. When P was less than or equal to 0.05, this indicated that the miRNA at least regulated one target, therefore it was defined as a biological dysfunctional miRNA. In step 2, GSEA was applied to identify significant pathways specific to BC based on Gene Ontology Biological Process (GOBP) (Subramanian et al., 2005). This method specified whether the pathways were randomly distributed at the top or bottom of the detected genes. The coefficients of Spearman correlation between mRNAs and sample label were defined as the weight of mRNAs. To obtain a robust result, the process was repeated by calculating ES 10,000 permutations. We selected pathways that include at least 50 genes to avoid random signals. The significant level of pathways was considered with levels of FDR ≤ 0.1 and P ≤ 0.05. In step 3, the modules of miRNA–mRNA were constructed based on the biological and statistically significant relationships. As the unclear distribution of microarray data, Spearman correlation coefficients indicated the relationships between miRNAs and between miRNAs and mRNAs. When FDR was less than 0.2, the association was considered as significant in both miRNA–miRNA and miRNA– mRNA. Constructing a module, we identified two topological features: its mRNAs both were the intersection of biological and statistical targets and belonged to the same pathways; its miRNAs were at least one biologically and statistically significant ones. In total, one module was constructed by two related miRNAs and one mRNA regulated by them.

−1 to 1. 0 indicating a random prediction and 1 determined a perfect prediction. If the value was less than 0, it was worse than a random prediction. AUC, ranging from 0 to 1, estimated the classification performance. If it was prior to 0.5, the model can discriminate the disease. The three indexes were calculated by five-fold cross validation. Cases and controls were randomly divided into 5 groups, respectively. These case and control groups combined into 5 datasets. Each dataset was in turn regarded as test set, and the residual ones were consider as training sets. In our flow, training sets were used to construct the relationship between miRNAs and mRNAs; however, test set verified the performance of modules.

3. Results 3.1. Summary of differentially expressed miRNAs and mRNAs Here, in the five different training sets, 40 common differentially expressed miRNAs were selected through SAM algorithm (Fig. 2). Through HGT, only 6 were defined as core miRNAs, which were biologically and statistically significant deregulated miRNAs, including miR106b, miR-15b, miR-424, miR-107, miR-125 and miR-16 (Table 1). Of 481 significant pairs between miRNAs, 12 pairs were between core ones, 142 were between core and non-core ones and residuals were between non-core ones, which were the pre-condition of the modules (Fig. 3, Additional file 1). Interestingly, miR-106b and miR-93, which belonged to the same gene cluster, were both statistically significant. Although different pathways were enriched by training sets, 7 pathways were identified by four times, which perhaps had tight correlation

2.3. Validation Considering the epidemiological significance, logistic regression was applied to fit the model. MCC, Accuracy and AUC estimated classification performance of module. MCC balanced the specificity and sensitivity. MCC was calculated by the following formula: TP  TN−FP  FN MCC ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðTP þ FNÞ  ðTN þ FPÞ  ðTP þ FPÞ  ðTN þ FNÞ

ð2Þ

in which TP, TN, FP and FN represented the number of true cases, true controls, false cases and false controls, respectively. MCC ranged from

Fig. 2. Volcano plots of miRNAs and mRNAs. X axis and y axis are identified as log2FC and -FDR, respectively. The dots far from the center are the significant dysregulated miRNAs and mRNAs. (A) Volcano plot of miRNAs; (B) volcano plot of mRNAs.

S. Yang et al. / Gene 548 (2014) 6–13 Table 1 The summary of deregulated miRNAs in cross validation.

9

Table 2 The frequency of deregulated pathways by GSEA.

No.

miRNA

Cross

FC

FDR

Pa

No.

GOBPID

Term

Freq

1 2 3 4 5 6 7 8 9

miR-106b miR-106b miR-15b miR-424 miR-106b miR-107 miR-125 miR-15b miR-16

1 3 3 4 5 5 5 5 5

−1.80 −1.82 −2.75 −3.38 −1.85 −3.90 7.50 −2.75 −0.83

0.05 0.06 0.06 0.09 0.02 0.09 0.08 0.06 0.09

5.41E−07 2.63E−05 5.47E−09 3.42E−02 7.43E−04 1.53E−03 0.00E+00 8.80E−04 0.00E+00

1 2 3 4 5 6 7 8 9 10 11 12 13

GO:0000075 GO:0044257 GO:0006281 GO:0000279 GO:0000087 GO:0007067 GO:0016071 GO:0006397 GO:0001116 GO:0006974 GO:0022613 GO:0006396 GO:0008380

Cell cycle checkpoint Cellular protein catabolic process DNA repair M phase M Phase of mitotic cell cycle Mitosis mRNA metabolic process mRNA processing Protein RNA complex assembly Response to DNA damage stimulus Ribonucleoprotein complex biogenesis and assembly RNA processing RNA splicing

4 4 4 4 4 4 4 4 4 3 2 2 2

a

The results from HGT.

to BC (Table 2). In these different pathways, we identified 705 different mRNAs. 3.2. Identification of modules of BC In the first training set, three different mRNAs, including CIT, KPNA2 and POLQ, were significantly related to miR-106b and miR-93 (Table 3).

These mRNAs were involved in five pathways: M phase of mitotic cell cycle, M phase, mitosis, DNA repair and response to DNA damage stimulus. CIT and KPNA2 co-regulated M Phase pathways, however, DNA repair and response to DNA damage stimulus both included POLQ. According to the topological features, miR-106b (OR = 8.23, 95% CI =

Fig. 3. The statistical correlations between core and non-core miRNAs in cross validation. The blue squares represent core miRNAs, and the gray circles represent non-core miRNAs. The red solid lines suggest the relationship between core miRNAs, and the black long dot lines indicate the relationship between core and non-core miRNAs. (A) Cross 1; (B) Cross 3; (C) Cross 4; (D) Cross 5.

10

S. Yang et al. / Gene 548 (2014) 6–13

Table 3 The significant biological relationships between miRNAs and mRNAs.

Cross1

miRNA

mRNA

Coef.

FDR

miRanda

TargetScan

miRTarBase

miR-106b

CIT KPNA2 POLQ CIT KPNA2 POLQ CIT KPNA2 POLQ BTRC UBR3 CIT EIF2B5 BTRC EIF2B5 UBR3 CIT POLQ KPNA2 EIF2B5 CIT CIT KPNA2 CIT KPNA2 POLQ

0.54 0.53 0.57 0.54 0.52 0.53 0.54 0.52 0.49 0.36 0.29 0.49 0.37 0.34 0.35 0.34 0.52 0.43 0.42 0.30 0.48 0.41 −0.28 0.49 0.50 0.43

2.01E−05 2.71E−05 5.53E−06 1.77E−05 5.63E−05 3.77E−05 1.68E−04 4.00E−04 1.28E−03 3.55E−02 1.17E−01 1.22E−03 3.40E−02 5.69E−02 4.92E−02 5.82E−02 4.24E−04 7.96E−03 9.30E−03 1.00E−01 1.80E−03 1.25E−02 1.31E−01 9.88E−04 7.88E−04 6.99E−03

True True True True True True True True True True True True True True True True True True True True True True True True True True

True True True True True True True True True True True True True True True True True True True True True True True True True True

False False False False False False False False False False False False False False True False False False False False False False False False False False

miR-93

Cross5

miR-106b

miR-107 miR-130b miR-15b miR-16

miR-17

miR-200c miR-301a miR-454 miR-485-3p miR-93

2.47–40.42), miR-93 (OR = 7.94, 95% CI = 2.30–43.86), CIT (OR =2.84, 95% CI = 2.84–62.88) and KPNA2 (OR = 2.49, 95% CI = 1.10–6.41) were involved in Module 1. Module 2 included miR-93, miR-106b and POLQ (OR = 3.39, 95% CI = 0.82–15.45) (Table 4). From the fifth training sets, we identified 13 different miRNAs and 6 mRNAs, including CIT, KPNA2, POLQ, BTRC, UBR3 and EIF2B5 (Table 3). BTRC and UBR3 regulated cellular protein catabolic process and protein catabolic process, and EIF2B5 were involved in ribonucleoprotein complex biogenesis and assembly. Therefore, we construct two new modules: Module 3 included miR-107 (OR = 15.64,

95% CI = 2.62–17.93), miR-16 (OR = 3.81, 95% CI = 1.38–13.21), BTRC (OR = 2.93, 95% CI = 1.18–8.50) and UBR3 (OR = 3.25, 95% CI = 0.93–16.11); Module 4 consisted of miR-16, miR15b (OR = 4.86, 95% CI =1.74–19.39), miR-200c (OR = 3.48, 95% CI = 1.48–9.63) and EIF2B5 (OR = 7.47, 95% CI = 1.32–58.53). Then the modules by the first training set included new molecules: miR-17, miR-130b (OR = 27.59, 95% CI = 4.44–429.60), miR-301 (OR = 6.18, 95% CI = 1.86–3.85), miR-454 (OR = 4.81, 95% CI = 1.84–17.36) and miR-4853p (OR =0.36, 95% CI = 0.14–0.75) for Module 1 and miR-17 for Module 2 (Table 4, Fig. 4). 3.3. Validation

Table 4 The performance of classification of significant molecules.

Cross 1

Cross5

a

Molecule

OR

95% CI

P

AUCa

CIT KPNA2 POLQ miR-106b miR-93 CIT KPNA2 POLQ BTRC UBR3 EIF2B5 miR-106b miR-107 miR-130b miR-15b miR-16 miR-17 miR-200c miR-25 miR-301a miR-454 miR-485-3p miR-92a miR-93

10.79 2.49 3.39 8.32 7.94 0.96 12.39 2.13 2.93 3.25 7.47 11.73 15.64 27.59 4.86 3.81 5.94 3.48 5.58 6.18 4.81 0.36 5.98 9.75

2.84–62.88 1.10–6.41 0.82–15.45 2.47–40.42 2.30–43.68 0.39–2.75 2.65–116.29 0.65–6.62 1.18–8.50 0.93–16.11 1.32–58.53 2.86–81.22 2.62–179.93 4.44–429.60 1.74–19.39 1.38–13.21 1.89–25.18 1.48–9.63 1.73–26.88 1.86–32.85 1.84–17.36 0.14–0.75 1.71–27.62 2.52–66.71

2.02E−03 4.01E−02 9.70E−02 2.34E−03 4.98E−03 9.30E−01 6.94E−03 1.89E−01 3.09E−02 9.68E−02 3.09E−02 3.23E−03 8.75E−03 3.44E−03 7.23E−03 1.63E−02 5.66E−03 7.34E−03 1.16E−02 1.14E−02 4.02E−03 1.22E−02 9.52E−03 5.42E−03

0.63 0.72 0.71 0.75 0.69 0.31 0.75 0.80 0.61 0.82 0.58 0.35 0.59 0.45 0.73 0.31 0.69 0.88 0.27 0.51 0.12 0.76 0.53 0.36

The AUC from testing sets.

In the five-fold cross validation process, AUC was applied to testify classification performance of single molecule and modules. Most molecules also discriminated the labels, however, miR-106 was unable to classify cases and controls (AUC = 0.35). Classification performance of modules was prior to that of the involved molecules (AUC = 0.84, 0.78, 0.86 and 0.88 for Module 1, Module 2, Module 3 and Module 4, respectively). The classification performance of Module 5 and Module 6 was inferior to a random prediction. MCC and Accuracy perhaps were unsuitable for this datasets, because number of controls was not enough (Table 5, Fig. 5). 3.4. Identification of modules from three subtypes BC was a highly heterogeneous cancer that was found in multiple subtypes, including LA, LB and BL (Nguyen et al., 2008). Using the subtype data, we further studied the modules of different subtypes by the same steps. For the different subtypes, we identified 12, 14 and 7 deregulated miRNAs, respectively. Interestingly, both miR-106b and miR-93 were identified as the statistically and biologically significant miRNAs in all three subtypes (Additional file 2). Meanwhile, we also identified multiple similar significant pathways and mRNAs to those of whole data, for example DNA repair. Interestingly, Module 1 perhaps acted as a non-specific module, which was involved in each subtype. However, each subtype also had the specific modules (Additional file 3).

S. Yang et al. / Gene 548 (2014) 6–13

11

Fig. 4. The regulatory relationships of four modules. For a module, the gray circles and triangles are core miRNAs and non-core miRNAs, respectively. Their relationships are represented by red. The blue squares show deregulated mRNAs in the same pathways, and their relationship between miRNAs are the black solid lines. (A) Module 1 consists of miR-106b, miR-93 and CIT and KPNA2; (B) Module 2 includes miR-106b, miR-93 and POLQ. (C) Module 3 involves miR-107, miR-16 and BTRC and UBR3. (D) Module 4 is assembled by miR-16, miR-15b, miR-200c and EIF2B5.

4. Discussion miRNAs are involved in many biological processes related to tumorigenesis, such as BC (Iorio et al., 2005; Ma et al., 2007). The multiple associations between miRNAs and mRNAs may cooperatively play essential roles in the generation and development of complex diseases. BC is the most common and most aggressive cancer in females whose generation includes the deregulation of different molecules, especially

Table 5 Association of AUCs, Accuracy and MCCs of modules. AUC

Accuracy

MCC

Cross 1 Module 1 Module 2

0.84 0.78

0.89 0.89

NaN NaN

Cross 5 Module 3 Module 4 Module 5 Module 6

0.86 0.88 0.47 0.45

0.90 0.95 0.75 0.65

NaN NaN 0.40 0.31

miRNAs and mRNAs (Stephens et al., 2012). Therefore, a method is needed that integrates the information of these molecules to explain the mechanism of complex biological phenomenon. In this study, the expression profiles of miRNAs and mRNAs, the predicted miRNA targets, and the statistical correlations between them were summarized and GSEA was used to promote the diagnostic ability of BC. The method has reconstructed four significant modules: miR-106bCIT-KPNA2-miR-93, miR-106b-POLQ-miR-93, miR-107-BTRC-UBR3miR-16 and miR-200c-miR-16-EIF2B5-miR-15b by reanalyzing microarray data (Fig. 4 and Table 3). Moreover, the modules are verified by the indexes of AUC, Accuracy and MCC, and they show that the classification performance of module is prior to that of the involving single molecules (Table 3 and 4). On the other hand, miR-106b-25 cluster, consisting of miR-106b, miR-93 and miR-25, is related to TGF-β in BC (Gong et al., 2013; Smith et al., 2012). miR-16 and miR-200c also contribute to BC by progestin-mediated oncogenic signaling and EGF, respectively (Rivas et al., 2012; Uhlmann et al., 2010). KPNA2 is a novel oncogenic factor in human BC (Noetzel et al., 2012). For the BL subtype, E2F1 is one of deregulated mRNAs, whose expression was different in luminal and basal-like (Xu et al., 2012). Although some evidence suggests that the single molecule is tight relative to BC, experiments can confirm the significance of modules.

12

S. Yang et al. / Gene 548 (2014) 6–13

Fig. 5. ROCs of four modules from testing sets. (A) Module 1; (B) Module 2; (C) Module 3; (D) Module 4.

Furthermore, compared to previous algorithms or workflows presented in the literature, our method has some advantages. Firstly, permutation test selects deregulated miRNAs considering sampling errors and unclear distribution of expression values of miRNAs. With the reciprocal regulation between miRNA and mRNA, HGT is further performed to identify significant deregulated miRNAs (Guo et al., 2010). Therefore, the profiles of differentially expressed miRNAs may be the causal miRNAs and the specific BC biomarkers (Guo et al., 2012b). Secondly, we also consider the statistical correlations between miRNAs and between miRNAs and mRNAs. The associations of miRNAs are also attenuated, because different miRNAs can synergistically or antagonistically regulate the same target roles (Guo et al., 2012a; Xu et al., 2011). Spearman correlation between miRNAs and mRNAs can identify specific modules from the network. Considering these correlations, the range of targets is narrowed and the accuracy of prediction is higher. Thirdly, GSEA is performed to select the significant and specific pathways because the association of mRNAs may contribute to the development and occurrence of disease prior to a single mRNA. In addition, it also reduces the false positive rate and decreases the multi-comparison of the coefficients of correlation. Finally, this method of constructing modules can apply to many complex diseases, because it only utilizes the expression levels of miRNAs and mRNAs. In summary, our method can detect not only the specific and significant miRNAs and mRNAs but also their interaction and the significant pathways they are involved in. With the development of microarray and high-throughput sequencing, the biological data (miRNA, lncRNA, methylation, SNP, etc.) is increasing in enormity and heterogeneity (Berger et al., 2013). The methods to mine these data presented a huge challenge. Moreover,

considering the increasing number of methods available, this one has some weaknesses. Firstly, other molecules were not considered, such as methylation. It is an epigenetic biomarker whose expression level is low; therefore it may provide more information about BC. Secondly, with the development of sequencing technologies this method only fits for microarray data. Therefore, further methods with integration of more molecular information and analysis of different data sources are warranted. In conclusion, we have proposed a new analysis flow to define the modules of miRNA–mRNA under a specific disease conditions by combining the biology and statistics. This method can apply to other complex diseases and supply new insight of the construction of modules. In a systematic view, the modules may contribute to the mechanism, diagnosis and prognosis of BC. Acknowledgments This work was supported by the National Natural Science Foundation of China (Nos. 61301251, 81072389 and 81373102), the Research Fund for the Doctoral Program of Higher Education of China (Nos. 211323411002 and 20133234120009), the China Postdoctoral Science Foundation funded project (No. 2012M521100), the key grant of the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (No. 10KJA33034), the National Natural Science Foundation of Jiangsu (No. BK20130885), the Natural Science Foundation of the Jiangsu Higher Education Institutions (Nos. 12KJB310003 and 13KJB330003), the Jiangsu Planned Projects for Postdoctoral Research Funds (No. 1201022B), the Science and Technology Development Fund Key Project

S. Yang et al. / Gene 548 (2014) 6–13

of Nanjing Medical University (No. 2012NJMU001), the Research and Innovation Projectfor College Graduates of Jiangsu Province (No. 944) and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.06.010. References Bartel, D.P., 2009. MicroRNAs: target recognition and regulatory functions. Cell 136, 215–233. Berger, B., Peng, J., Singh, M., 2013. Computational solutions for omics data. Nat. Rev. Genet. 14, 333–346. Che, X., Huang, C., 2012. microRNA, cancer and cancer chemoprevention. Curr. Mol. Pharmacol. 5, 362–371. Davis, B.N., Hilyard, A.C., Nguyen, P.H., Lagna, G., Hata, A., 2010. Smad proteins bind a conserved RNA sequence to promote microRNA maturation by Drosha. Mol. Cell 39, 373–384. Enerly, E., Steinfeld, I., Kleivi, K., Leivonen, S.K., Aure, M.R., Russnes, H.G., Ronneberg, J.A., Johnsen, H., Navon, R., Rodland, E., Makela, R., Naume, B., Perala, M., Kallioniemi, O., Kristensen, V.N., Yakhini, Z., Borresen-Dale, A.L., 2011. miRNA–mRNA integrated analysis reveals roles for miRNAs in primary breast tumors. PLoS One 6, e16915. Enright, A.J., John, B., Gaul, U., Tuschl, T., Sander, C., Marks, D.S., 2003. MicroRNA targets in Drosophila. Genome Biol. 5, R1. Gong, C., Qu, S., Liu, B., Pan, S., Jiao, Y., Nie, Y., Su, F., Liu, Q., Song, E., 2013. MiR-106b expression determines the proliferation paradox of TGF-β in breast cancer cells. Oncogene. http://dx.doi.org/10.1038/onc.2013.525. Griffiths-Jones, S., Grocock, R.J., van Dongen, S., Bateman, A., Enright, A.J., 2006. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 34, D140–D144. Guo, H., Ingolia, N.T., Weissman, J.S., Bartel, D.P., 2010. Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature 466, 835–840. Guo, L., Sun, B., Wu, Q., Yang, S., Chen, F., 2012a. miRNA–miRNA interaction implicates for potential mutual regulatory pattern. Gene 511, 187–194. Guo, L., Zhao, Y., Yang, S., Cai, M., Wu, Q., Chen, F., 2012b. Genome-wide screen for aberrantly expressed miRNAs reveals miRNA profile signature in breast cancer. Mol. Biol. Rep. 40, 2175–2186. Guo, L., Ji, X., Yang, S., Hou, Z., Luo, C., Fan, J., Ni, C., Chen, F., 2013. Genome-wide analysis of aberrantly expressed circulating miRNAs in patients with coal workers' pneumoconiosis. Mol. Biol. Rep. 40, 3739–3747. Hsu, S.-D., Lin, F.-M., Wu, W.-Y., Liang, C., Huang, W.-C., Chan, W.-L., Tsai, W.-T., Chen, G.-Z., Lee, C.-J., Chiu, C.-M., 2011. miRTarBase: a database curates experimentally validated microRNA–target interactions. Nucleic Acids Res. 39, D163–D169. Hwang, H.W., Mendell, J.T., 2007. MicroRNAs in cell proliferation, cell death, and tumorigenesis. Br. J. Cancer (96 Suppl.), R40–R44. Iorio, M.V., Ferracin, M., Liu, C.-G., Veronese, A., Spizzo, R., Sabbioni, S., Magri, E., Pedriali, M., Fabbri, M., Campiglio, M., 2005. MicroRNA gene expression deregulation in human breast cancer. Cancer Res. 65, 7065–7070. Kong, W., He, L., Richards, E.J., Challa, S., Xu, C.X., Permuth-Wey, J., Lancaster, J.M., Coppola, D., Sellers, T.A., Djeu, J.Y., Cheng, J.Q., 2013. Upregulation of miRNA-155 promotes tumour angiogenesis by targeting VHL and is associated with poor prognosis and triple-negative breast cancer. Oncogene 33, 679–689. Lewis, B.P., Burge, C.B., Bartel, D.P., 2005. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell 120, 15–20. Liu, B., Liu, L., Tsykin, A., Goodall, G.J., Green, J.E., Zhu, M., Kim, C.H., Li, J., 2010. Identifying functional miRNA–mRNA regulatory modules with correspondence latent dirichlet allocation. Bioinformatics 26, 3105–3111. Lussier, Y.A., Stadler, W.M., Chen, J.L., 2012. Advantages of genomic complexity: bioinformatics opportunities in microRNA cancer signatures. J. Am. Med. Inform. Assoc. 19, 156–160. Ma, L., Teruya-Feldstein, J., Weinberg, R.A., 2007. Tumour invasion and metastasis initiated by microRNA-10b in breast cancer. Nature 449, 682–688. Miska, E.A., Alvarez-Saavedra, E., Abbott, A.L., Lau, N.C., Hellman, A.B., McGonagle, S.M., Bartel, D.P., Ambros, V.R., Horvitz, H.R., 2007. Most Caenorhabditis elegans microRNAs are individually not essential for development or viability. PLoS Genet. 3, e215. Nguyen, P.L., Taghian, A.G., Katz, M.S., Niemierko, A., Raad, R.F.A., Boon, W.L., Bellon, J.R., Wong, J.S., Smith, B.L., Harris, J.R., 2008. Breast cancer subtype approximated by

13

estrogen receptor, progesterone receptor, and HER-2 is associated with local and distant recurrence after breast-conserving therapy. Proc. Am. Soc. Clin. Oncol. 2373–2378. Noetzel, E., Rose, M., Bornemann, J., Gajewski, M., Knüchel, R., Dahl, E., 2012. Nuclear transport receptor karyopherin-α2 promotes malignant breast cancer phenotypes in vitro. Oncogene 31, 2101–2114. Ota, H., Sakurai, M., Gupta, R., Valente, L., Wulff, B.E., Ariyoshi, K., Iizasa, H., Davuluri, R.V., Nishikura, K., 2013. ADAR1 forms a complex with Dicer to promote microRNA processing and RNA-induced gene silencing. Cell 153, 575–589. Peng, X., Li, Y., Walters, K.-A., Rosenzweig, E.R., Lederer, S.L., Aicher, L.D., Proll, S., Katze, M.G., 2009. Computational identification of hepatitis C virus associated microRNA–mRNA regulatory modules in human livers. BMC Genomics 10, 373. Png, K.J., Yoshida, M., Zhang, X.H., Shu, W., Lee, H., Rimner, A., Chan, T.A., Comen, E., Andrade, V.P., Kim, S.W., King, T.A., Hudis, C.A., Norton, L., Hicks, J., Massague, J., Tavazoie, S.F., 2011. MicroRNA-335 inhibits tumor reinitiation and is silenced through genetic and epigenetic mechanisms in human breast cancer. Genes Dev. 25, 226–231. Pritchard, C.C., Kroh, E., Wood, B., Arroyo, J.D., Dougherty, K.J., Miyaji, M.M., Tait, J.F., Tewari, M., 2012. Blood cell origin of circulating microRNAs: a cautionary note for cancer biomarker studies. Cancer Prev. Res. (Phila.) 5, 492–497. Rivas, M.A., Venturutti, L., Huang, Y.-W., Schillaci, R., Huang, T.H., Elizalde, P.V., 2012. Downregulation of the tumor-suppressor miR-16 via progestin-mediated oncogenic signaling contributes to breast cancer development. Breast Cancer Res. 14, R77. Smith, A.L., Iwanaga, R., Drasin, D.J., Micalizzi, D.S., Vartuli, R.L., Ford, H.L., 2012. Abstract A20: the Six1-regulated miR-106b-25 cluster is a mediator of the tumor promotional effects of TGF-β signaling in human breast cancer. Cancer Res. 72, A20-A20. Stephens, P.J., Tarpey, P.S., Davies, H., Van Loo, P., Greenman, C., Wedge, D.C., Nik-Zainal, S., Martin, S., Varela, I., Bignell, G.R., Yates, L.R., Papaemmanuil, E., Beare, D., Butler, A., Cheverton, A., Gamble, J., Hinton, J., Jia, M., Jayakumar, A., Jones, D., Latimer, C., Lau, K.W., McLaren, S., McBride, D.J., Menzies, A., Mudie, L., Raine, K., Rad, R., Chapman, M.S., Teague, J., Easton, D., Langerod, A., Lee, M.T., Shen, C.Y., Tee, B.T., Huimin, B.W., Broeks, A., Vargas, A.C., Turashvili, G., Martens, J., Fatima, A., Miron, P., Chin, S.F., Thomas, G., Boyault, S., Mariani, O., Lakhani, S.R., van de Vijver, M., van 't Veer, L., Foekens, J., Desmedt, C., Sotiriou, C., Tutt, A., Caldas, C., Reis-Filho, J.S., Aparicio, S.A., Salomon, A.V., Borresen-Dale, A.L., Richardson, A.L., Campbell, P.J., Futreal, P.A., Stratton, M.R., 2012. The landscape of cancer genes and mutational processes in breast cancer. Nature 486, 400–404. Subramanian, A., Tamayo, P., Mootha, V.K., Mukherjee, S., Ebert, B.L., Gillette, M.A., Paulovich, A., Pomeroy, S.L., Golub, T.R., Lander, E.S., Mesirov, J.P., 2005. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550. Takeshita, N., Hoshino, I., Mori, M., Akutsu, Y., Hanari, N., Yoneyama, Y., Ikeda, N., Isozaki, Y., Maruyama, T., Akanuma, N., Komatsu, A., Jitsukawa, M., Matsubara, H., 2013. Serum microRNA expression profile: miR-1246 as a novel diagnostic and prognostic biomarker for oesophageal squamous cell carcinoma. Br. J. Cancer 108, 644–652. Tavazoie, S.F., Alarcon, C., Oskarsson, T., Padua, D., Wang, Q., Bos, P.D., Gerald, W.L., Massague, J., 2008. Endogenous human microRNAs that suppress breast cancer metastasis. Nature 451, 147–152. Tusher, V.G., Tibshirani, R., Chu, G., 2001. Significance analysis of microarrays applied to the ionizing radiation response. Proc. Natl. Acad. Sci. U. S. A. 98, 5116–5121. Uhlmann, S., Zhang, J., Schwäger, A., Mannsperger, H., Riazalhosseini, Y., Burmester, S., Ward, A., Korf, U., Wiemann, S., Sahin, Ö., 2010. mir-200bc/429 cluster targets plcγ1 and differentially regulates proliferation and egf-driven invasion than mir200a/141 in breast cancer. Oncogene 29, 4297–4306. Valastyan, S., Reinhardt, F., Benaich, N., Calogrias, D., Szasz, A.M., Wang, Z.C., Brock, J.E., Richardson, A.L., Weinberg, R.A., 2009. A pleiotropically acting microRNA, miR-31, inhibits breast cancer metastasis. Cell 137, 1032–1046. Volinia, S., Galasso, M., Sana, M.E., Wise, T.F., Palatini, J., Huebner, K., Croce, C.M., 2012. Breast cancer signatures for invasiveness and prognosis defined by deep sequencing of microRNA. Proc. Natl. Acad. Sci. U. S. A. 109, 3024–3029. Xiao, Y., Ping, Y., Fan, H., Xu, C., Guan, J., Zhao, H., Li, Y., Lv, Y., Jin, Y., Wang, L., Li, X., 2013. Identifying dysfunctional miRNA–mRNA regulatory modules by inverse activation, cofunction, and high interconnection of target genes: a case study of glioblastoma. Neuro Oncol. 15, 818–828. Xu, J., Li, C.X., Li, Y.S., Lv, J.Y., Ma, Y., Shao, T.T., Xu, L.D., Wang, Y.Y., Du, L., Zhang, Y.P., Jiang, W., Li, C.Q., Xiao, Y., Li, X., 2011. miRNA–miRNA synergistic network: construction via co-regulating functional modules and disease miRNA topological features. Nucleic Acids Res. 39, 825–836. Xu, K., Usary, J., Kousis, P.C., Prat, A., Wang, D.-Y., Adams, J.R., Wang, W., Loch, A.J., Deng, T., Zhao, W., 2012. Lunatic fringe deficiency cooperates with the Met/Caveolin gene amplicon to induce basal-like breast cancer. Cancer Cell 21, 626–641.