A new method to study the change of miRNA-mRNA interactions due to environmental exposures.

MOTIVATION: Integrative approaches characterizing the interactions among different types of biological molecules have been demonstrated to be useful for revealing informative biological mechanisms. One such example is the interaction between microRNA (miRNA) and messenger RNA (mRNA), whose deregulation may be sensitive to environmental insult leading to altered phenotypes. The goal of this work is to develop an effective data integration method to characterize deregulation between miRNA and mRNA due to environmental toxicant exposures. We will use data from an animal experiment designed to investigate the effect of low-dose environmental chemical exposure on normal mammary gland development in rats to motivate and evaluate the proposed method.
RESULTS: We propose a new network approach-integrative Joint Random Forest (iJRF), which characterizes the regulatory system between miRNAs and mRNAs using a network model. iJRF is designed to work under the high-dimension low-sample-size regime, and can borrow information across different treatment conditions to achieve more accurate network inference. It also effectively takes into account prior information of miRNA-mRNA regulatory relationships from existing databases. When iJRF is applied to the data from the environmental chemical exposure study, we detected a few important miRNAs that regulated a large number of mRNAs in the control group but not in the exposed groups, suggesting the disruption of miRNA activity due to chemical exposure. Effects of chemical exposure on two affected miRNAs were further validated using breast cancer human cell lines.
AVAILABILITY AND IMPLEMENTATION: R package iJRF is available at CRAN. CONTACTS: pei.wang@mssm.edu or susan.teitelbaum@mssm.edu. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

With the past decades rapid technological advances in genomic profiling, multiple types of genomic profiles can be collected on the same set of samples. Recently, we carried out an animal study to assess the effect of three environmental chemicals on miRNA and mRNA activity in mammary tissue (Gopalakrishnan ). In this and similar studies, there is increasing interest to characterize changes in the regulatory patterns among different molecular types across different experimental conditions. Compared to commonly used marginal analyses, integrative approaches examining interactions often help to reveal more subtle yet biologically important mechanisms. For instance, it is well known that miRNAs drive the development of many diseases via the regulation of post-transcriptional gene expression (Jansson and Lund, 2012; Nogales-Cadenas ). Thus it will be more powerful to characterize miRNA activities through monitoring the global regulatory system between miRNA and mRNA (Arner and Kulyté, 2015).

Multiple challenges arose during the construction of the high dimensional miRNA–mRNA interaction networks. First, such analysis involves thousands or tens of thousands of genes but a much smaller sample size (n = 20 in the motivating chemical exposure study). This problem, well known in statistics as the ‘large p, small n’ paradigm (Bernardo ), arises in many biological applications (Kosorok ). Second, samples under different conditions in the same study often share common properties, and borrowing information across conditions is crucial to maximize the power of the estimation process (Flutre ; Li ; Petretto ). Third, properly incorporating information from existing databases in the analysis could greatly enhance the accuracy of the inference (Bernard and Hartemink, 2004; Petralia ; Werhli ; Yip ; Zhu ).

However, recent methodologies for integrating different genomic profiles (for example, Mo ; Hwang ; Ebrahim ; Ritchie ) either were not designed for characterizing interaction networks, or do not fully address the above challenges, especially the latter two. Recently, random forest based methods have been utilized to estimate multiple networks simultaneously (Petralia ). In particular, we demonstrated the advantage of joint learning and showed the high performance of random-forest compared to Gaussian graphical models such as Danaher . In this work, we will extend these ideas to model the relationship of dependence among miRNAs and mRNAs. Moreover, in Petralia , we introduced a probability sampling scheme to effectively incorporate prior information in building random forest models. This framework can be effectively modified to incorporate existing knowledge on miRNA–mRNA regulatory relationships documented in miRNA–mRNA databases (Agarwal ; Betel ; Hsu ; Kertesz ).

Therefore, in this paper, we propose a new method—integrative Joint Random Forest (iJRF), which borrows information across multiple chemical exposure conditions and takes into account prior information from existing databases when inferring miRNA–mRNA interactions. iJRF is built upon our two previous random-forest based algorithms (Petralia , 2016). The advantages of our integrative framework are multiple. First, its ensemble nature allows the delivery of excellent performance with moderate sample size requirements. Second, treatment-specific tree ensembles are designed to share common structures, so that miRNAs regulations playing a crucial role in multiple conditions will be detected more accurately. Third, existing databases are utilized in order to prioritize miRNA–mRNA interactions. Specifically, TargetScan (Agarwal ) was chosen over other databases (Betel ; Kertesz ) for its comprehensive list of predicted miRNA–mRNA interactions and, as mentioned in Lee , its consistency with databases containing experimentally validated targets (Farazi ; Helwak ).

Applying iJRF to the chemical exposure study, we simultaneously estimated miRNA–mRNA interaction networks for the control group and three chemical exposed groups: diethyl phthalate (DEP), methyl-paraben (MPB) and triclosan (TCS). We found that the interaction among miRNAs and mRNAs was greatly reduced in the chemical exposed groups compared to the control group. Among all chemicals, DEP exposure was associated with the highest loss of connectivity. iJRF also detected two important miRNAs: miR-200a and miR-375, that played crucial roles in the inferred regulatory network of the control group, but lost more than 90% connectivity in the network corresponding to the DEP exposure group. mRNAs connecting with miR-200a and miR-375 in control network only, were enriched for the ‘Mammary Gland Development’ and ‘Gland Morphogenesis’ pathways. This suggests a mediating role for these miRNAs associated with chemical exposure in mammary gland development. We then confirmed the effect of DEP on miR-375 and miR-200a using human breast cancer cell lines.

2 Materials and methods

2.1 Random forest for network construction

Random forest is a non-linear algorithm that models the response variable via a series of decision trees where each tree is constructed based on a random subset of samples (Breiman, 2001). At each node, a random subset of predictors is considered and the predictor maximizing a certain utility function (i.e. decrease in node impurity) is chosen to split observations into two subsets. Recently, Huynh-Thu introduced GENIE3, a random forest based model for inferring gene regulatory networks (GRNs). In GENIE3, first, the expression of each target gene k is modeled as a function of the expression of all other genes via random forest, then, the regulatory events are ranked based on random forest importance scores. Specifically, the importance score is defined as the total decrease in node impurities from splitting on the jth predictor, averaged over all trees. Recently, Petralia proposed iRafNet—a new random-forest based algorithm for network construction which can integrate prior information from database and independent datasets. According to iRafNet, potential regulators considered important by other datasets are prioritized and sampled more often within the random-forest framework. Petralia extended the original random-forest algorithm to estimate multiple related networks (JRF). As shown by Petralia , borrowing information across multiple networks is crucial to accurately detect common mechanisms. In particular, information across different class of data is borrowed by using the same splitting variables for the tree construction. In this paper, JRF and iRafNet are combined to jointly estimate miRNA–mRNA interactions from different exposure conditions while integrating information from existing databases.

2.2 iJRF: integrative joint random forest

We are interested in inferring miRNA - mRNA interactions in tissues samples for control and three common environmental chemicals: diethyl phthalate (DEP), methyl paraben (MPB) and triclosan (TCS). Denote as the index over different treatment conditions. For each treatment condition g, we observe the expression of M miRNAs and p mRNAs for n samples. Denote and as the expression of the kth mRNA and the jth miRNA for the ith sample exposed to the gth treatment.

An overview of the proposed algorithm is shown in Figure 1. Specifically, for each treatment condition g, the expression of the kth mRNA is modeled as a function of the expression of miRNAs via random forest, i.e. . Therefore, for each target mRNA, G random forest models corresponding to G treatment conditions are constructed. The key idea of iJRF is to build G tree ensembles simultaneously. Let τ denote the current node in the gth tree model corresponding to the gth condition. As illustrated in Figure 1, the allocation processes for are performed simultaneously through the following steps:

Fig. 1

Joint Random Forest with iRafNet sampling scheme. For each exposure condition, model the expression of mRNAs as function of the expression of miRNAs via random forest. At each node, sample miRNAs prioritizing those present in TargetScan (Agarwal ). Following JRF model (Petralia ), the four random forest tree ensembles (Control, DEP, MPB and TCS) use the same splitting variables (miRNAs) to build trees. In this way we achieve borrowing information across them. This procedure is repeated for each mRNA and, then, interactions are ranked based on random forest importance scores

iRafNet Step (Petralia et al., : For different tree ensembles corresponding to different treatments, the same set of predictors (miRNAs) are proposed for the splitting rule of nodes . This subset of predictors is selected by prioritizing miRNAs that have similar sequences to that of the target mRNA and, thus, are more likely to bind to the target mRNA. Specifically, we sample a set containing N miRNAs from the entire set of miRNAs with probabilities

where scores are derived based on prior information on miRNA–mRNA regulations from existing databases. In this study, TargetScan (Agarwal ) was chosen over other databases given its comprehensive list of predicted miRNA–mRNA interactions. Also, as mentioned in Lee , TargetScan is consistent with databases containing experimentally validated targets (Farazi ; Helwak ). For each interaction , TargetScan provides a context score based on sequence similarity (Garcia ). Context scores are non-positive with more negative values corresponding to more favorable sites. We then calculate: with and being the minimum context score. Scores take values in the interval [1, 2] and the probability to sample the miRNA with the most similar sequence to that of the target mRNA will be twice the probability of the least similar miRNA (context score equal to zero). This transformation was considered in order to not excessively penalize interactions that are not contained in prior databases. with being the decrease in node impurity observed in the gth tree after splitting τ based on the jth predictor. Specifically, the decrease in node impurity is defined as where is the variance of observations allocated to set is the set of samples allocated to node τ in the gth tree ensemble; while sets and are respectively the sets of samples allocated to the left-child and right-child of node τ according to a splitting rule based on the jth predictor.

2. JRF Step (Petralia et al., : Among predictors contained in subset , the optimal splitting variable of nodes is the predictor maximizing the summation of the decrease in node impurity across different treatment conditions, i.e.,

Once G random forest models are constructed for the kth mRNA, interactions between miRNAs and the kth mRNA under different conditions are ranked based on random forest importance scores. In order to derive the final unweighted networks, a proper cut-off value for importance scores needs to be chosen. Specifically, we follow the same permutation based procedure described in Petralia , whose details are provided in Supplementary Section 1.

The computational complexity of the proposed algorithm is in the same order as the complexity of JRF (Petralia ), i.e. with T being the number of random forest trees. The computational burden can be greatly reduced by estimating in parallel random forest models corresponding to different target mRNAs (further information can be found in Supplementary Section 2.1). In practice, the number of trees (T) and the number of potential regulators to be sampled at each node (N) are parameters to be specified by the users. In this paper, we used the conventional choice of T = 1000 and with M being the number of predictors (Breiman, 2001).

3 Data

3.1 Data overview

Exposure to environmental chemicals, especially during mammary gland development, has been linked to breast cancer risk using animal models (Manservisi ; Moral ; Rudel ; Russo ). Increased understanding of the biological and genomic mechanisms mediating the effects of chemical exposures can lead to better prevention and treatment of the disease. In a recent pioneering study using a Sprague-Dawley (SD) rat model, animals were exposed chronically from birth through adulthood (postnatal day 146) to three environmental chemicals present in daily personal care products: diethyl phthalate (DEP), methyl paraben (MPB) and triclosan (TCS) at an exposure level comparable to those observed in human population (Teitelbaum ). The study involved 20 rats in each of the following groups, control, DEP and MPB treatment groups, as well as 15 rats in the TCS treatment group. Gopalakrishnan provides a detailed description of the animal study and the chemical treatment for MPB and TCS exposure. Likewise other chemicals, diethyl phthalate (DEP) (CAS 84-66-2, lot STBB0862V, 99% purity) was supplied in plastic containers (Sigma Aldrich, Italy). The experimental oral dose of DEP was 0.1735 mg/Kg/day, which represented of no observed adverse effect levels (NOAEL) of DEP (Brown ; Moody and Reddy, 1978; Oishi and Hiraga, 1980).

3.2 Data processing

The expression matrices of both mRNAs and miRNAs can be found in the GEO database (ID: GSE72276) (Barrett ). Quality control of.CEL files and preprocessing based on the robust multiarray average method (RMA) were done using the expression console software (Affymetrix, CA). For mRNA data, batch effects were removed using the ComBat package (Leek ) available in R CRAN. We applied a signal intensity filter to retain only those probesets with high and stable expression (signal value > 30th percentile in at least 1 of the experimental groups). A variance-based filter was used to retain the top 50% of the probesets with high interquartile range. For miRNA data, 283 miRNAs were profiled and only miRNAs with variance different from zero were considered. The data were normalized using the package NanoStringNorm available in R CRAN (Waggott ). The filtered data contained 7546 genes and 272 miRNAs which were used for downstream analyses. For both miRNA and mRNA data matrices, quantile normalization across samples was performed.

4 Results

4.1 Network estimation

For ease of explanation, we will refer to the network from control rats as Control-Net and networks from different chemical exposed rats as DEP-Net, MPB-Net and TCS-Net. For the analyses, we considered 7546 messenger RNAs and 272 miRNAs. As mentioned in section 3, the sample size was 20 each for the control, DEP and MPB treatment groups, and 15 for TCS. In order to implement the proposed algorithm, 1000 trees were considered and a total of miRNAs (predictors) were sampled at each node (Breiman, 2001). The four networks were estimated using iJRF and mRNA-miRNA interactions were derived using permutation techniques considering an FDR cut-off of 0.001 (further information can be found in Supplementary Section 1).

Table 1 shows the total number of interactions inferred for each network as well as the number of interactions shared across networks. As shown, all three chemicals result in a loss of interaction compared to Control-Net which involved 6829 edges linking 47 miRNAs and 2270 mRNAs. In particular, the total number of interactions in DEP-Net, MPB-Net and TCS-Net were respectively 44%, 84% and 52% of the total number of interactions in Control-Net. Therefore, DEP was the chemical exposure resulting in the most dramatic loss of interaction compared to control.

Table 1

Number of interactions inferred in Control-Net, DEP-Net, MPB-Net and TCS-Net and number of interactions shared across networks

	Control	DEP	MPB	TCS
Control	6829	2079	3593	2374
DEP		3018	842	1630
MPB			5743	3491
TCS				3557

Figure 2(a) shows the top 10 hub-miRNAs in Control-Net, which were responsible for more than 85% of connecting edges in Control-Net. In particular, a substantial portion of those interactions () were not present in any of the chemical-networks. Figure 2(b) compares each chemical-Net and Control-Net showing, for each miRNA, the number of edges shared between chemical-Net and Control-Net (green bar), the number of control-specific edges (blue bar) and the number of chemical-specific edges (red bar). The three quantities have been normalized dividing them by the total number of connecting edges present in either Control-Net or chemical-networks. As shown, DEP-Net has miRNAs such as miR-375-3p, miR-200a-3p and miR-214-3p with remarkable loss in connectivity () compared to Control-Net. Given the dramatic loss in connectivity observed in DEP-Net, we decided to focus on this chemical and miRNAs miR-375-3p, miR-200a-3p and miR-214-3p for further investigation. Besides the consistent loss of connection in chemical-networks compared to Control-Net, differences among chemical-networks were also observed, as shown in Supplementary Section 2.2. Further investigation is needed to understand the biological implication of these differences across chemical-networks.

Fig. 2

(a) Plot of top ten hub-miRNAs in Control-Net with interactions detected only in Control-Net (Shannon ). These miRNAs were responsible for more than 85% of connecting edges in Control-Net. The total number of interactions detected in Control-Net was 6829. (b) For each miRNA, we show the number of edges shared by chemical and control (green bar), the number of control-specific edges (blue bar) and the number of chemical-specific edges (red bar). The three quantities have been normalized dividing them by the total number of connecting edges in either Control-Net or chemical-networks

It is important to note that, no significant results were detected by traditional univariate analysis (i.e. unpaired t-test, Wilcoxon test) when testing was conducted on one miRNA at-a-time, suggesting the advantage of the proposed network based approach (further information on univariate analysis based on the Wilcoxon test can be found in Supplementary Section 3.1). As a comparison, we also estimated miRNA–mRNA interactions under different treatment conditions via Pearson’s correlation test, a commonly used approach in the literature (Mertins ; Zhang , 2016). As shown in Supplementary Section 3.2, when considering the same FDR cut-off (fdr = 0.001) as in the iJRF analysis, the correlation test detects far fewer edges than iJRF, and fails to reveal any informative hub structure or pathway enriched network module. We then relaxed the FDR cut-off to 0.01 for the correlation tests to obtain more connected networks, and tested for enriched GO terms. As shown in Supplementary Figure S6, the correlation test resulted in far fewer enriched GO terms compared to iJRF. This suggests that iJRF is more effective to detect biologically relevant interactions. Moreover, the percentage of shared edges across networks based on correlation tests is much smaller than that based on iJRF, which makes the detection of treatment-specific interactions particularly vulnerable to high false positive rates. This result is expected since iJRF, through joint learning, is more effective in detecting common associations than algorithms handling different treatment conditions separately (Petralia ).

4.2 miR-375, miR-214 and miR-200a

Overview: One of the top hub miRNAs in Control-Net poorly connected in chemicals was miR-375-3p. In particular, the number of mRNAs interacting with miR-375-3p in DEP-Net, MPB-Net and TCS-Net was, respectively, 2%, 51% and 27% of the number of interacting mRNAs in Control-Net. Recent studies have investigated the role of miR-375 in breast cancer (Madhavan ; Ward ; Zehentmayr ). In particular, Ward reported a loss of miR-375 expression in drug-resistant breast cancer cells; while Madhavan showed that miR-375 was significantly associated with breast cancer survival. As shown in Figure 2(b), other miRNAs poorly connected in DEP-Net were miR-214 and miR-200a. Their role in breast cancer has been investigated in several papers (Kalniete ; Ming ; Penna ; Pieraccioli ; Wang ; Yao ; Yu ). In particular, miR-214 was shown to be associated with breast cancer survival (Kalniete ) and drug sensitivity (Yu ). On the other hand, miR-200a has been associated with survival in metastatic breast cancer (Madhavan ) and its role in cell proliferation inhibition has been demonstrated (Yao ).

Overlap with MirTarBase: Some Control-specific interactions detected by iJRF are contained in miRTarBase (Hsu )—the database of experimentally validated miRNA–mRNA interactions. Table 2 shows the list of interactions in Control-Net contained in miRTarBase for miR-375 and miR-200a (other interactions contained in miRTarBase can be found in Supplementary Table S2). As shown in Table 2, some of these interactions were contained in DEP-Net, MPB-Net and TCS-Net as well. Particularly interesting is the interaction of miR-375 with the human epidermal growth factor receptor 2 (HER2) (Pillai ; Shen ) only contained in Control-Net and MPB-Net. HER2 stimulates the growth of breast cancer cells and is one of the main targets for breast cancer survival and therapy. Another interesting interaction is miR-200a - ZEB2 which has been investigated in different ovarian and breast cancer studies (Ahmad ; Bracken ; Jang ; Park ; Truong ; Wu ). In particular, many articles have described ZEB2 as the crucial target of miR-200 family members (Burk ; Christoffersen ; Korpal ). Brabletz and Brabletz (2010) showed that ZEB2 and miR-200a are involved in a ‘feedback loop’ which drives the progression of metastasis in breast cancer. As shown in Table 2, the interaction between ZEB2 and miR-200a is only contained in Control-Net, suggesting the impact of all chemicals on the regulatory mechanism between miR-200a and ZEB2.

Table 2

List of interactions in Control-Net contained in miRTarBase for miR-200a and miR-375

miRNAs	mRNAs	DEP	MPB	TCS
miR-375	HER2, TMTC4, SFT2D2, KRT8		x
miR-375	PLAG1, CCDC88A, CELF2			x
miR-375	GATA6	x	x
miR-375	CMTM4, FOLR1, CTSC
miR-200a	ZEB2
miR-200a	HOXB5			x
miR-200a	DLC1	x

For each interaction, we indicate if it was contained in other networks such as DEP-Net, MPB-Net and TCS-Net.

Enrichment Analysis: Among hub-miRNAs in Control-Net, interesting pathways were obtained for miR-375-3p and miR-200a-3p. Figure 3(a) shows some interesting enriched categories for mRNAs connected to miR-375-3p and miR-200a-3p in Control-Net but not in DEP-Net (enriched pathways for other hub-miRNAs in Control-Net can be found in Supplementary Table S3). The enrichment analysis was performed using David Tool 6.7 (Huang ) and Benjamini adjusted p-values were reported in Figure 3(a). As shown, miR-375-3p and miR-200a-3p share some enriched categories such as ‘Mammary Gland Development’ with Benjamini adjusted p-values less than 0.1. This result is not surprising since the two miRNAs share some connections in Control-Net. Various studies have shown that chemical exposure can alter mammary gland development (Manservisi ; Mandrup ; Schwarzman ). In this context, DEP exposure might alter the regulatory mechanism of miR-375-3p and miR-200a-3p and affect mammary gland development. As shown in Figure 3(a), enriched pathways include genes such as ERBB2 (HER2), FOXA1 and SFRP1 which play a crucial role in breast cancer. To further demonstrate the loss in connectivity in DEP-Net, Figure 3(b) shows the correlation density between each miRNA and mRNAs connected in Control-Net but not in DEP-Net. As expected, the correlation density for DEP exposure is shifted to the left compared to that of control revealing a loss of correlation. Supplementary Figure S7 shows that the loss of correlation observed in DEP-Net for both miR-375-3p and miR-200a-3p is significant.

Fig. 3

(a) We consider genes connected to miR-375-3p and miR-200a-3p in Control-Net but not in DEP-Net and derived enriched categories using David Tools (Huang ). Pathways ‘Gland Development’, ‘Plasma Membrane’ and ‘Mammary Gland Development’ were enriched for both miR-375-3p and miR-200a-3p with Benjamini adjusted p-values smaller than 0.10. Pathway ‘Gland Morphogenesis’ was enriched only for miR-375-3p. (b) Density of absolute correlation between miR-375-3p and miR-200a-3p with mRNAs connected only in Control-Net for DEP exposed data (red) and control data (black)

4.3 Validation via cell line experiments

In Section 4.1,we showed that DEP-Net resulted in a dramatic loss of interaction compared to Control-Net. In particular, the three hub- miRNAs poorly connected in DEP-Net were miR-375, miR-214 and miR-200a. To validate the effect of DEP, in vitro experiments of a human breast cancer cell line are utilized.

Experiments: MCF-7 cells were maintained in phenol-red-free DMEM (Gibco 11054) containing 5% (v/v) dextran-charcoal-stripped fetal calf serum, with M diethyl phthalate, for 4 weeks. These cells were subcultured every 3-4 days with a confluence level less than 70%. The total RNA was isolated from 10 cm Petri plate using the Promega Maxwell simplyRNA kit according to the manufacturer’s instructions. Reverse transcription was done using the Exiqon Universal cDNA Synthesis kit and the detection of qPCR was performed with the Exiqon ExiLENT SYBR Green master mix on Roche LightCycler 480 machine.

Results: We quantified the expression of miR-200a, miR-375 and miR-214 in MCF-7 cells. Unfortunately, miR-214 was not expressed in this cell line, thus only results for miR-200a and miR-375 are reported. As shown in Figure 4, the expression levels of both miR-375 and miR-200a are significantly different between control and DEP-exposed cells. This result suggests that the two miRNAs are affected by DEP exposure in breast cancer human cell lines. Further analyses are necessary to elucidate the role of these two miRNAs on mammary gland development.

Fig. 4

Plot of cell line experiments. Expression levels of miR-375 and miR-200a for normal and DEP exposed cells. Red diamonds indicate average over the three replicates. Benjamini adjusted p-values from the unpaired t-test are reported

5 Discussion

In this paper, we focused on exposure to chemicals commonly used in personal care products during mammary gland development (birth to adulthood) and their effect on miRNA–mRNA interactions. Some challenges arose given the large number of miRNAs/genes involved and the limited sample size of our data. In order to overcome this problem, we proposed iJRF, an ensemble algorithm with small sample size requirement. The advantage of iJRF is dual. First, it is designed to borrow information across different treatment conditions so that associations shared across treatments can be detected more accurately. Second, iJRF can integrate information from existing miRNA–mRNA databases such as TargetScan.

Using our newly developed algorithm, we estimated miRNA–mRNA interaction in mammary tissues for control and chemical-exposed animals. All chemical-networks registered a loss in connectivity compared to control. In particular, DEP was the chemical registering the most dramatic loss of interaction compared to control. In fact, the total number of interactions in DEP-Net, MPB-Net and TCS-Net were respectively 44%, 84% and 52% of the total number of interactions in Control-Net.

Among the leading miRNAs in Control-Net, miR-200a, miR-214 and miR-375 lost more than 90% of connectivity in DEP-Net compared to Control-Net. Messenger RNAs connected to miR-200a and miR-375 in Control-Net but not in DEP-Net were enriched in ‘Gland Morphogenesis’ and ‘Mammary Gland Development’, indicating their potential involvement in mammary gland development mechanisms. Among genes in these pathways, we found targets in breast cancer such as ERBB2, FOXA1 and SFRP1. Recent studies have investigated the role of miR-375, miR-214 and miR-200a in breast cancer. For example, Ward reported a loss of miR-375 expression in tamoxifen-resistant breast cancer cells; while Madhavan reported a significant association between survival in metastatic breast cancer and both miR-375 and miR-200a. The expression of miR-214 was shown to be associated with breast cancer survival (Kalniete ) and drug sensitivity (Yu ); while Yao demonstrated the role of miR-200a in the inhibition of cell proliferation in breast cancer.

Given the dramatic loss of connectivity observed in DEP-Net, we validated the effect using cell line experiments. Specifically, we focused on miRNAs with the highest loss in connectivity: miR-375, miR-214 and miR-200a. Using MCF-7 cells, we showed that the expression levels of both miR-375 and miR-200a were significantly different between the control and DEP exposed groups. In this study, we hypothesized that chemical exposures first affect miRNA activity, which then changes the regulatory pattern among miRNAs and mRNAs. In our validation experiments, we demonstrated the effects of chemical exposure on miRNA activity. Future research is warranted to further validate the detected changing of regulatory relationships among miRNAs and mRNAs.

In this paper we examined miRNA–mRNA networks for single exposure conditions. On the other hand, humans are exposed to a mixture of different chemicals with the interaction and combination of chemicals playing a crucial role. In order to deal with such complex data, as future work, we will design a model to estimate networks which vary across exposure conditions in a dynamic way. Once chemical-induced mechanisms will be identified, their association to breast cancer phenotypes will be assessed.

Finally, the proposed algorithm can be utilized in different biological applications. As an example, eQTL analysis might be performed for different tissues simultaneously while borrowing information from existing databases. It is well known that borrowing information across tissues is crucial to detect shared eQTLs more accurately (Flutre ) and, iJRF is a non-parametric model that can be easily implemented for such analyses.

Click here for additional data file.