Modeling microRNA-mRNA Interactions Using PLS Regression in Human Colon Cancer

Outline of our analysis

We used a multi-step approach in our analysis. In STEP 1, an un-paired t-test was used with FDR correction methods to identify differentially expressed miRNAs and mRNAs in colon tumors and normal colon tissues. In STEP 2, the PLS regression method and statistical tests on the association measures were performed. Repeated bootstrap samples were used to identify significant associations/interactions between miRNA and mRNA using 31 down-regulated miRNAs and 71 up-regulated mRNAs. This PLS-bootstrap algorithm work-flow is outlined in the Figure 1. In STEP 3, we compared the predictions from our PLS regression method with those from TargetScan 5.1 and miRanda. In STEP 4, a gene ontology analysis was constructed with the aid of Pathway Studio software. In STEP 5, miRNA-mRNA interacting networks were developed with the aid of Cytoscape 2.7 software. In STEP 6 and 7, IPA and IPA pathway designer were used to create a functional analysis of biological processes related to miRNAs and their target mRNAs. Finally, in STEP 8, a simulation was performed to evaluate the false detection rate for PLS. Steps 1-3 can be considered as the discoveries after mining the real data and Steps 6-8 correspond to the validities of these mined results.

STEP 1. Differentially expressed miRNA and mRNA in human colon tissues

The miRNA and mRNA microarray expression profiles from the same tissues including 7 human colon tumor tissues and 4 normal tissues were analyzed with the aid of GeneSpring GX10 software (Agilent, CA). The Benjamini & Hochberg FDR (false discovery rate) was used to obtain adjusted p-values for unpaired t-tests [28]. Differentially expressed miRNAs and mRNAs with adjusted p-values less than or equal to 10% and the absolute fold changes greater than or equal to 1.2 were obtained. We found a total of 31 down-regulated miRNAs and 2 up-regulated miRNAs that fit these criteria (Additional file 1, table S1). In addition, we found 73 up-regulated mRNAs and 63 down-regulated mRNAs (Additional file 2, table S2). Some of the mRNAs did not have gene IDs. There were only two up-regulated miRNAs using the above mentioned cutoff. Therefore, in this study we selected the 31 down-regulated miRNAs and the 71 up-regulated mRNAs to find the most likely miRNA:mRNA negative/inverse associations with the aid of the PLS regression method. Among the differentially expressed miRNAs, only miR-135b and miR-221 are significantly up-regulated (4.3 and 1.7 fold, respectively) in the colon cancer tissue samples. The remaining 31 miRNAs (Additional file 1, table S1) were significantly down-regulated in these samples.

STEP 2. miRNA:mRNA Associations predicted by PLS

To identify the most likely prospective targets the 31 down-regulated miRNAs were used as independent variables and the 71 mRNAs were denoted as dependent variables in the PLS regression analysis. Subsequently, statistical tests were performed to identify significant negative associations between each mRNA with the 31 down-regulated miRNAs. The algorithm flow diagram for this multi-step procedure is shown in Figure 1.

STEP 3. Comparisons between PLS Regression, TargetScan 5.1 and miRanda

Comparing the PLS regression method with other computational miRNA-target prediction algorithms such as TargetScan 5.1 and miRanda we found some areas of agreement. Thus, some of the miRNA-targets predicted by PLS were also predicted by either TargetScan or miRanda or both of them (Table 2). This table provides the number of miRNA-targets for those miRNAs reportedly linked to cancer [29–42]. It shows the percentage of common targets predicted by our use of the PLS regression method, with TargetScan 5.1 and miRanda.

Using 15 down-regulated miRNAs (let-7 g, miR-101, miR-133a, miR-150, miR-15a, miR-16, miR-29b, miR-29c, miR-30a, miR-30b, miR-30c, miR-30d, miR-30e, miR-34b and miR-342), known to be associated with cancer, we found 16.5% and 11.0% of our PLS-predicted miRNA-targets, on average, were also predicted as targets for the corresponding miRNAs by TargetScan5.1 and miRanda, respectively (Table 2). Moreover, we found that five out of 25 of our PLS-predicted targets (20%) for miR-30d were also predicted as targets for this miRNA by both TargetScan5.1 and miRanda. These five target genes predicted by our PLS regression approach, TargetScan 5.1 and miRanda were BHLHE40, DDIT4, NEBL, SLC7A11 and UBE2I. In this regard, it is interesting to note that our method identified more common targets with TargetScan than miRanda. From Tables 1 and 2 it is evident that our analysis provides many novel targets that were not identified by these existing algorithms which may prove to be significant biomarkers for colon cancer. In order to further pursue the integrative analysis an association network for these significantly associated miRNAs and mRNAs was constructed (see STEP 5).

STEP 4. Gene ontology (GO) terms for miRNA-targets in the category of biological processes

STEP 5. Visualization of miRNA:mRNA association networks

Using the statistically significant β _ij -association scores between the i-th mRNA and j-th miRNA, derived from the PLS regression method, we constructed association networks including the significant miRNAs and mRNAs. The resulting miRNA:mRNA association network had 97 nodes and 155 connections between the 31 down-regulated miRNAs and the 71 up-regulated mRNAs (Figure 2A). We developed a network to demonstrate the overlapping miRNA targets for miRNAs in the miR-30 family (a-e) (Figure 2B) since these miRNAs have a large number of mRNA targets. The number of targets for miR-30 a-e was 22, 20, 15, 25 and 17, respectively (Table 3). This network (Figure 2B) indicated that two mRNAs (VKORC1 and CSTB; in red) were predicted to be targets of five miR-30 miRNAs, six mRNAs (ANXA2, GPATCH3, OAZ1, PSMB4, SLC7A11, ZNF37A; in turquoise) were targeted by four of the miRNAs, and ten mRNAs (CAPG, C7orf28A, DAP, GPX1, ITPA, MMP11, NEBL, POLR2H, S100A11 and GTF2IRD1; 3 in green and 7 in blue) were targeted by three of the miRNAs (Figure 2B). Networks were also developed for the seven miRNAs (let-7 g, miR-101, miR-133a, miR-15a, miR-16, miR-29b and miR-29c) closely related to cancer and their associated mRNAs (Figure 2C). Most of the targets mentioned here have been reported to be linked to cancer. For example, ANXA2 up-regulation was reported to be associated with human hepatocellular carcinoma, pancreatic adenocarcinoma, high-grade glioma, gastric carcinoma, acute promyelocytic leukemia and primary colorectal cancer [43]. Also, CapG message and protein are reported to be up-regulated in colorectal cancer [44].

STEP 6. Pathway analysis with the aid of IPA

To gain additional insight into prospective miRNA linked biological networks and canonical pathways, the 31 down-regulated miRNAs from the colon cancer microarrays were analyzed with the aid of Ingenuity Pathway Analysis (IPA) software. We focused on two miRNA-associated networks related to cancer. The most significant (top) functional pathways highlighted by IPA consisted of 1) cancer, reproductive system disease, and genetic disorders and 2) cell death, hematological system development and function, and cancer. We found that ten of the down-regulated miRNAs (miR101, miR26a, miR26b, miR30a, miR30b, miR30d, miR30e, miR34b, miR-let7 g and miRN140) were grouped together in a functional network (Figure 3A) and nine of the down-regulated miRNAs (miR-130a, miR-133a, miR-142, miR-150, miR15a, miR-16, miR-29b, miR-30c and miR-99a) were grouped together in a second network (Figure 3B). These are colored in green in Figures 3A and 3B.

STEP 7. Combinatorial analysis with the aid of the PLS method and IPA pathway designer

For the above pathway analysis results from 31 down-regulated miRNAs and 71 up-regulated mRNAs in the previous step, we were interested in knowing which miRNAs were associated with the mRNAs grouped in the above cancer-related networks, functions and pathways predicted by our PLS regression method. With the aid of IPA pathway designer we expanded the above networks, functions and pathways generated by IPA (Figure 3, Tables 4, 5 and 6) by combining the PLS regression method and IPA pathway.

At first, we examined the two networks that were mediated by miRNAs and that were generated from the 31 down-regulated miRNAs (Figure 3A and 3B). With the aid of the IPA pathway designer, we connected miRNA-targets predicted by the PLS regression method with each associated miRNA labeled as my list (ML in Figure 4A and 4B). These miRNAs in ML were connected to the two networks (Figure 3A and 3B) through their predicted targets by IPA. The new expanded networks from Figure 3A and 3B are shown in Figure 4A and 4B. In Figure 4A, we found that 27 of 31 down-regulated miRNAs were associated with the miRNA-associated network through 15 genes (ANXA2, ARPC1B, BHLHE40, GPX1, MIF, OAZ1, ODC1, PHB, PSMB4, RAN, RPL10, RPL34, RPS9, TUFM, UBE2I) predicted by the PLS regression method. These mRNAs are all included in the list of predicted target in Table 2. It was interesting to see that the four remaining miRNAs (miR-100, rno-miR-140, miR15a and miR-26a) were grouped in this network by IPA and not connected through the above cancer-related target mRNAs by the pathway designer and hence two of them (miR-100 and miR-15a) were not included in Figure 4A. In this network, many of the cancer-affected target mRNAs were mainly regulated by the transcription factors MYCN and FOXO1 which play a vital role in cell death, apoptosis, proliferation and survival. Particularly MYCN was involved in the ERK/MAPK signaling pathway and cancer.

In a similar manner in Figure 4B, we found that 12 of 31-miRNAs were grouped with the miRNA-associated network through their predicted targets (BHLHE40, GPX1, MDK, S100A9 and UBE2I). These cancer-associated target mRNAs were connected to this network mainly through the interaction with apoptosis and cell death genes such as BCL2, CASP9 and ELK1 (Figure 4B).

Next we examined the expanded IPA results from 71 up-regulated mRNAs which were shown in the last column of Tables 4, 5 and 6. At first we wanted to know which of the 31 down-regulated miRNAs were associated with the top four networks through their target mRNAs predicted by PLS regression method. With the aid of IPA pathway designer, we found that 27 of the 31 down-regulated miRNAs were linked to one or more mRNA networks and 20 of them (let-7 g, miR-101, miR-126, miR-133a, miR-142-5p, miR-150, miR-15a, miR-26b, miR-28, miR-29b, miR-30a, miR-30b, miR-30c, miR-30d, miR-30e, miR-34b, miR-99a, mmu-miR-151, mmu-miR-342 and rno-miR-151) were involved in all of the top 4 networks. By restricting our attention to only the 15 cancer-associated miRNAs in the top four mRNA networks, we found that all 15 miRNAs were involved in network 1, all but miR-16 were in network 2, and all but miR-29c were in network 3 and 4, as shown in the fourth column of Table 4.

We also checked which miRNAs associated with the mRNA targets predicted by PLS regression method were linked to the cancer-related function. Through their targets, we found that 26 out of the 31 down-regulated miRNAs were associated with these targets and were in some ways related to cancer, 24 miRNAs were in the pathway of polyamine regulation in colon cancer and 12 miRNAs were involved in angiogenesis pathway. Again, we were interested in examining the 15 miRNAs directly linked with cancer-and their biological functions. We found that all 15 miRNAs were involved in cancer and tumorigenesis, 12 of them (all except miR-101, miR-15a and miR-29c) were in carcinoma, malignant tumor and primary tumor and 8 of them (all except let-7 g, miR-101, miR-150, miR-15a, miR-16, miR-29c and miR-342) were in angiogenesis as were shown in the last column of Table 5.

Furthermore, we examined which associated miRNAs among the 15 cancer-related miRNAs were involved in the canonical pathways associated with cancer. We found that 14 miRNAs were associated with protein ubiquitination pathway, all 15 miRNAs were in polyamine regulation in colon cancer, 11 miRNAs were associated with purine metabolism, 2 miRNAs were associated with circadian rhythm signaling, 4 miRNAs each were associated with hypoxia signaling in cardiovascular system and cdc42 signaling, and 7 miRNAs each were associated with pyrimidine metabolism, mTOR signaling or NRF2-mediated oxidative stress response or xenobiotic metabolism signaling along and their targets are mentioned in the last column of Table 6.

Overall, it clearly showed that these cancer-related miRNAs may mediate many biological functions and pathways by regulating the target genes predicted by our integrative and computational method.

STEP 8. Estimating false detection rate for PLS based target prediction

In order to determine the false detection rate of our mRNA target prediction method using PLS regression we used a simulation experiment. As defined in [22], the false detection rate is the percentage of detected significant miRNA-mRNA pairs derived as output from the PLS method out of the total possible pairs that would have the same or better associations by chance. Based on this procedure, the estimated false detection rate for the threshold FDR of 0.3 was 2.50%. We then reported the estimated false detection rates for many values of the thresholds simultaneously, and the estimated false detection rate is plotted as a function of the threshold in Figure 5.

We also compared the PLS regression method with a simpler method based only on the correlations for all miRNA-mRNA pairs. A correlation threshold of -0.750 corresponded to an estimated false detection rate of 2.50%. The PLS based method was more sensitive in detecting targets when applied to the actual data set; only 33 miRNA-mRNA connections were declared to be significant based on the correlation analysis with a threshold -0.750, whereas 155 connections were significant for the FDR threshold of 0.3.

Integrative computational modeling of miRNA targets using PLS

Although there are many computational algorithms for predicting miRNA targets, they vary widely with regards to specific targets identified. The relative usefulness of these algorithms is strengthened if there are some agreements in their conclusions. However, that is usually not the case among the existing algorithms. Additionally, the existing algorithms depend on the sequence complementarity of the miRNA-mRNA pairs. So, more complex algorithms using joint expressions of miRNA and mRNA may be useful to predict miRNA targets. In this study we developed an approach for predicting miRNA-mRNA associations based upon a statistical analysis using unpaired t-tests, PLS regression based association scores and the known inverse relationship between certain mRNAs and miRNAs.

Our efforts focused on a set of 31 significantly down-regulated miRNAs and a set of 71 significantly up-regulated mRNAs identified by the t-test for colon cancer vs. the control samples after multiple hypotheses corrections. We selected these miRNAs and mRNAs based on their inverse relationship for our integrative analysis. Using miRNAs as independent variables and mRNAs as dependent variables in our PLS regression method we identified significant miRNA-mRNA associations and in turn predicted mRNA targets for individual miRNAs (Table 1). When comparing these to miRNA-targets predicted by TargetScan 5.1 and miRanda, we found on average about 16.5% and 11.0% of our predicted targets matched those predicted by these two programs, respectively. Furthermore, some of our predicted target genes lie within the union or intersection of the other two target prediction programs. It seems likely that the miRNA target genes that lie within such intersections could be more accurate, although this claim remains to be evaluated further, experimentally. Interestingly, many of the targets predicted by our methods are not predicted by the other prediction algorithms and thus our approach provides possible novel targets that need further investigation for relevance to colorectal cancer {Step 7}. We estimated the false detection rate of the proposed method through a simulation study by sampling from the real data under the null model. A limitation of the present study is that we are unable to compute statistical performance measures such as the false discovery rate since there does not exist any realistic simulator currently for generating miRNA and mRNA datasets simultaneously with controlled biological interactions between them.

We have also compared our computational method of target prediction with another computational method of target prediction by association resulting from pairwise simple correlations between all 31 the up-regulated miRNAs and 71 down-regulated mRNAs. For the same dataset using the same threshold we get many more targets and so the PLS based method proved to be more sensitive than just correlation based target prediction described in Step 8 in the Results section.

miRNA and mRNA microarray data and gene expression profile analysis

The miRNA and mRNA microarray data sets for both human colon tumors and normal tissues used here were obtained from the Broad Institute and downloaded from their web-accessible database: http://www.broadinstitute.org/cgi-bin/cancer/datasets.cgi[20]. The data sets originally consisted of five normal human tissues and 10 human colon tumors for both miRNA and mRNA. However, only four of the normal tissue samples and seven of the colon tumor samples passed quality control criteria [46]. The mRNA expression values were obtained using Affymetrix GENECHIP analysis software (after hybridization to Affymetrix oligonucleotide microarrays Hu6800 and Hu35KsubA) containing a total of 16,063 probes [46]. The miRNA median fluorescence intensity values were measured by using a new, bead-based flow cytometric miRNA expression profiling method. All data were filtered by minimum value 32 (values less than 32 were considered as a background value according to the reference) and log₂ transformed [20].

The miRNA and mRNA expression data sets for colon cancer patients were loaded into GeneSpring GX 10.0 (Agilent Technologies). An unpaired t-test was used to compare miRNA and mRNA expressions for normal and tumor tissue samples. The p-values obtained were then adjusted with the aid of the Benjamini-Hochberg false discovery rate (FDR) [28] correction method. The differentially expressed miRNAs or mRNAs with q-values less than or equal to 10% and the absolute fold changes greater than or equal to 1.2 were determined to be significant.

miRNA target prediction using partial least square (PLS) method and bootstrapping

An algorithmic version of the whole procedure is provided in Figure 1. A detailed description of the PLS regression method used in this paper was previously reported in [23] and [27]. Briefly, consider a data set with n observations of p miRNAs (predictor variables) and m mRNAs (response variables). Let Y = (Y _1i ,Y _2i , ...,Y _ni )' be the column vector of observed values for the ith mRNA;i = 1,2, ... ,m, and let X _j = (X _1j , X _2j , ..., X _nj )' be the column vector of observed values for the jth miRNA for j = 1,2, ...,p. For each mRNA Y _i , PLS fits a regression model of the form

where e _i is an n-dimensional vector of residuals and is the k-th latent variable amongst the set of all orthogonal latent variables (l<<<p), all of which are linear combinations of the covariates. All of the original variables are standardized before doing PLS.

Using ordinary least squares leads to

and then the coefficients associated with the original covariates can be obtained from the coefficients of the latent variables.

Specifically, after q latent factors are found by selecting the number of component as 3, we calculated the scaled score coefficient . Then the score coefficients and intercept can be obtained by transforming back to the original scale as follows ,

and

. These coefficients are called the PLS scores, and these scores are defined as the associations between each response variable (mRNA) and the miRNA covariate for a specific i and j.

Statistical test of the significance of the association between miRNA and mRNA (β _ij )

In order to test the significance of the association between miRNA and mRNA, we tested the null hypothesis is H₀: β _ij = 0, (which implies there is no significant association between a specified miRNA and mRNA) versus the alternative H_a: β _ij < 0( which implies there is an inverse association between them). As we do not know the asymptotic distribution of such associations we used a bootstrap re-sampling scheme to construct the sampling distribution of the score β _ij for testing the hypothesis.

We ran PLS after removing the jth miRNA gene, to generate new β _ij and we compute the residual of the model which is the estimate of the error e _i .

After the residuals are computed, bootstrapping residuals are generated by simple random sampling with replacement of the components e _i1,e _i2,...,e _in of the residual vector e _i . Then the bootstrapped data can be found by the original regression equation:

where the parameters hold their usual meanings. We then run PLS with the new and original full set of miRNAs, and store the bootstrapped PLS based score . We did this bootstrapping step 1000 times for each gene and calculate the proportion of times when the negative bootstrapped association scores were smaller than negative observed association score . These proportions are the estimated p-values of the tests. The Benjamini & Hochberg FDR [28] then is used to calculate q-values. A threshold FDR was arbitrarily chosen as 0.3 since the number of variables in this test was small after having identified the inverse relationship between miRNA and mRNA. The associations for which the q-values were less than equal to the FDR cutoff were deemed to be significant.

Simulation Experiment

A simulation method similar in concept to the one proposed in [22] was used to estimate the false detection rate of possible targets using PLS regression for a specific threshold of FDR. The definition of the false detection rate was provided before in Step 8 of the Results section, as described on page 4 of [22]. Note that the concept of false detection rate refers to the simulation-based estimate of the percentage of detected significant miRNA-mRNA pairs derived by chance, and when using this concept, we refer to it in full as false detection rate so as not to confuse it with the abbreviation FDR which we have already reserved for the more commonly-used false discovery rate [28]. While we still use the FDR to control false positives for our PLS method, we use the false detection rate for this simulation to determine equivalent thresholds for our method and for a simple correlation method which does not choose its threshold based on FDR.

To generate a simulated data set with n "unassociated" observations of p miRNAs and m mRNAs, we randomly permuted the order of the n observations for both the miRNA genes and the mRNA genes so that none of the miRNA arrays corresponded to its true mRNA array. Specifically, we randomly permuted the labels of both arrays and rejected any generated data sets for which the new labels matched. Therefore, since the miRNAs do not correspond to the mRNAs in the simulated data set, any low q-values only occur by chance. For the simulated data set, we then compute the number of false positives by counting the number of q-values which fell below the threshold. The process was repeated for 1000 simulated data sets and the estimated false detection rate was computed by dividing the total number of false positives for the 1000 data sets by the overall total number of miRNA-mRNA pairs for the 1000 data sets.

We also compared the PLS regression with a simpler method based only on the correlations for all miRNA-mRNA pairs. In the correlation method, we declare that a miRNA-mRNA pair is significant if the corresponding pairwise correlation fell below a certain threshold. To compare the sensitivity of the methods, we used the same simulation method described above, and we chose the threshold for which the estimated false detection rate was the same as that found for the PLS regression method.

Biological function and pathway analysis

To obtain gene ontology (GO) terms for each miRNA, we used Pathway Studio software (Ariadne Genomics Inc., Rockville, MD, USA). This software enabled an enrichment analysis of expressed genes for GO terms using the differentially expressed miRNA target genes with respect to the category of biological process. The GO terms were selected by significant p-value and cancer-associated biological functions.

Using Ingenuity Pathway Analysis System and pathway designer (IPA 8.5, Ingenuity Systems, Edwood, CA, USA), we performed more detailed functional analysis to identify miRNA-mediated, cancer-related and statistically significant networks, biological functions and canonical signaling pathways for both differentially expressed miRNAs and mRNAs.