Cepred: Predicting the Co-Expression Patterns of the Human Intronic microRNAs with Their Host Genes

Identifying the tissues in which a microRNA is expressed could enhance the understanding of the functions, the biological processes, and the diseases associated with that microRNA. However, the mechanisms of microRNA biogenesis and expression remain largely unclear and the identification of the tissues in which a microRNA is expressed is limited. Here, we present a machine learning based approach to predict whether an intronic microRNA show high co-expression with its host gene, by doing so, we could infer the tissues in which a microRNA is high expressed through the expression profile of its host gene. Our approach is able to achieve an accuracy of 79% in the leave-one-out cross validation and 95% on an independent testing dataset. We further estimated our method through comparing the predicted tissue specific microRNAs and the tissue specific microRNAs identified by biological experiments. This study presented a valuable tool to predict the co-expression patterns between human intronic microRNAs and their host genes, which would also help to understand the microRNA expression and regulation mechanisms. Finally, this framework can be easily extended to other species.


Introduction
MicroRNA (miRNA) is a class of small non-coding RNAs (,22 nt) identified in recent years. They usually function as the negative regulators for certain genes at the post-transcriptional level by binding to the 39UTRs of the target mRNAs through base pairing, resulting in the target mRNAs cleavage or translation inhibition [1,2,3]. It has been newly shown that miRNAs may also function as positive regulators [4,5]. It is estimated that 1-4% of the genes in human genome are miRNAs, a single miRNA can regulate more than 200 targets [6], and the miRNAs preferentially regulate duplicated genes in mammals [7]. Increasing evidences suggest that the miRNAs play crucial roles in nearly all the important biological processes, such as cell growth, proliferation, differentiation, development, and apoptosis [6]. It has been reported that miRNAs also participate in the cellular signaling networks [8] and the gene regulatory networks [9], as well as the cross-species gene expression variation [10] and pathways [11]. Hereby, the miRNAs might be associated with various diseases [12].
To get a more comprehensive understanding of the miRNAs, it becomes important to identify the tissues in which a miRNA is expressed, which would promote not only the understanding of the miRNA biogenesis and expression mechanisms but also the prediction of their functions. However, it is difficult to identify the tissues in which a miRNA is expressed because great limitations exist in the current large-scale expression profile detecting techniques, such as the microarray [13]. For example, it is difficult to detect the expression levels of those miRNAs in rare cells or expressed at low levels [13]; the amount of miRNAs contained in a microarray chip is relatively limited comparing to the total amount of miRNAs in the human genome. Moreover, microarray experiments are both time and cost consuming. Most importantly, detecting the miRNA expressions of human directly is more difficult than other species since the lack of tissues. Therefore, it becomes critical to develop new approaches to decipher the expression of miRNAs. miRNAs show a very biased distribution in genomes. More than 36% (196/533) of the human miRNAs are within the introns (referred as the intronic miRNA here) of protein-coding genes (the host genes, miRBase Version 10, October 2007). The enrichment of the intronic miRNAs in human genome triggered researchers to wonder whether the intronic miRNAs and their host genes have special connections in terms of expression and regulation. Indeed, several groups have reported that many intronic miRNAs show significantly correlated expression profiles with their host genes [14,15,16,17]. One possible explanation for this phenomenon is that (1) the intronic miRNAs may be co-transcribed with their host genes by inclusion introns of their pre-mRNAs [17]; (2) the intronic miRNAs and their host genes may share common regulatory elements, such as common promoters [18]. By performing a comparison between the known miRNA expression profiles [19] and the RT-PCR expression profiles of the host genes, Rodriguez et al. found perfect correlations between the expression profiles of the investigated intronic miRNAs and their host genes [17]. These findings seemed to provide us a chance to uncover the tissues in which an intronic miRNA is expressed by associating with the expression profiles of its host gene. However, there are some intronic miRNAs which are low co-expressed with their host genes [14], which makes the problem much more complicated. Therefore, the first step to infer the expression tissues of the intronic miRNAs in this way is to identify which intronic miRNAs show high co-expression with the host genes and which do not.
In this study, we established a machine learning based approach to predict whether an intronic miRNA is high co-expressed with its host gene. We tested our method by the leave-one-out cross validation scheme, a novel independent dataset, and tissue specific miRNAs. Finally, we applied our approach to all the human intronic miRNAs.

Results
The co-expression dataset of Baskerville et al. [14] was employed as the training dataset. We first calculated the feature vectors for each sample in the training dataset, which was used to train a classifier. The performance of the classifier was then evaluated. Finally, we applied the classifier to all the human intronic miRNAs. As a result, all human intronic miRNAs showing high co-expression with host genes are identified and then the expression profiles of these miRNAs are inferred by identifying the expression profiles of the host genes. The complete flowchart of our framework is shown in Figure 1.

Genomic features of the human intronic miRNAs
In this study, the different co-expression patterns of the intronic miRNAs with their host genes triggered us to explore what the reasons are. Zhou et al. reported that the miRNA hosting introns have a 59-biased relative position distribution compared to all the other introns in human and mouse genomes, suggesting that the cis-signals within the 59UTRs of the host genes may interfere the transcription and regulation of the intronic miRNAs to certain extent [20]. Moreover, this finding suggests that the relative positions of the intronic miRNAs in their host genes is not completely random and may have functional significance. Therefore, here we present four feature vectors, which reflect the positional characteristics of the intronic miRNAs in their host genes to train the classifier (Table 1). These four features describe the position of miRNA hosting intron, the position of miRNA, the length of the miRNA hosting intron, and the ratio of the length of miRNA to that of the host gene.

Support vector machine training and predicting
A support vector machine (SVM) based classifier was implemented to predict the co-expression levels between the intronic miRNAs and their corresponding host genes. The dataset from Baskerville et . [14] was adopted as the training set. The feature vectors (Table 1) for each sample were computed and are listed in Table 2. As shown in Table 2, columns ''1'', ''2'', ''3'', and ''4'' list the values of these four features for each miRNA; column ''coexpressions'' lists the correlation coefficients of the miRNA and host gene expression profiles; column ''ER'' lists the co-expression patterns of miRNA and host gene identified by biological experiments, in which ''+'' represents high co-expression and ''2'' represents low co-expression. A leave-one-out cross validation scheme was employed to evaluate the performance of this SVM classifier. For each round of the leave-one-out experiments, one sample was eliminated from the training set and the remaining samples were used to train the SVM classifier. The prediction was made on the uncovered sample with the well-trained model. Each sample was left out for one time and the prediction results for all the samples are listed in the column ''PR'' of Table 2. The results of this experiment revealed that our approach achieves an accuracy of 79% and shows good sensitivity and specificity (Table 3), which indicates that the co-expression patterns are predictable using SVM and these four features. The accuracy is calculated as the ratio of the number of the samples that are predicted correctly to the number of the whole samples. The sensitivity of high/low co-expression class is calculated as the ratio of the number of the samples that are predicted as high/low coexpression correctly to the total number of true high/low coexpression samples, respectively. The specificity of high/low coexpression class is calculated as the ratio of the number of the samples that are true high/low co-expression to the number of the samples that are predicted as high/low co-expression, respectively.
In order to confirm the results, we further test our method on an independent dataset. This testing dataset is combined from one set of mRNA microarray data and two sets of miRNA microarray data from three separated labs (see Materials and Methods). There are 5 and 16 intronic miRNAs which showed high or low coexpression patterns with their host genes, respectively (Table 4), in which we evaluated the co-expression patterns of miRNA and host  More specifically, all the 5 experimental high co-expressed intronic miRNAs, as well as 15 out of the 16 low coexpressed intronic miRNAs were correctly predicted. The missing one was hsa-mir-149, which showed low co-expression in the experimental results but was predicted as high co-expressed by our classifier. Finally, we got a total accuracy of 95% in this validation, which indicates that our approach is robust (Table 3).
Identifying the tissues in which the human intronic miRNAs are expressed Using the 29 miRNA-host gene pairs of Baskerville et al., we have constructed a SVM based classifier to predict whether an intronic miRNA is high co-expressed with its host gene. Here we applied this classifier on the rest 178 human intronic miRNAs reported in miRBase database [21] to predict whether they are high co-expressed with their host genes. Finally, 64 and 112 of them were predicted to be high or low co-expressed with their host genes, respectively (Supplementary Text S1). Together with the 12 high co-expressed intronic miRNA-host gene pairs in the training dataset, we obtained 76 high co-expressed intronic miRNA-host gene pairs. The expression profiles of these intronic miRNAs could be estimated by associating with the expression profiles of their host genes. By examining the expression profiles of the host genes across 79 human tissues [22], we obtained the expression profiles of 71 miRNAs in these 79 human tissues (Supplementary Text S2).  intronic miRNAs and show high co-expression with host genes in our prediction. Table 5 lists the specific tissues identified by prediction and biological experiments for these 12 miRNAs, respectively. As a result, 11 of the 12 overlaps were successfully predicted (Table 5), in which mir-302 cluster contains mir-302a, mir-302b, mir-302c, mir-302d, and mir-367. For example, we predicted that hsa-mir-488 is high expressed in the amygdala, temporal lobe, globu spallidus, cerebellum peduncles, cerebellum, caudate nucleus, whole brain, parietal lobe, medulla oblongata, prefrontal cortex, occipital lobe, hypothalamus, thalamus, subthalamic nucleus, cingulated cortex, pons, spinal cord, fetal brain, adrenal gland, adrenal cortex and pituitary, which mainly belong to nervous systems except the adrenal gland and adrenal cortex ( Figure 2). As shown in the Figure 2, each bar represents one tissue; the height of the bar represents the expression level of hasmir-488 in that tissue. The tissues with high expression level are highlighted in red and blue color, in which red bars represent tissues of nervous systems. The exception is mir-1-1, which was reported to be heart specific by Landgraf et al. and was predicted to show high specificity in skeletal muscle in our study. However, mir-1 was also reported to show high specificity in both heart and skeletal in several other studies, which indicates that all the above 12 predictions in our study are right [24,25].

Validation through tissue-specific miRNAs
Furthermore, we also predicted that hsa-mir-208b shows high specificity in heart, skeletal muscle and tongue ( Figure 3). As shown in the Figure 3, each bar represents one tissue; the height of the bar represents the expression level of has-mir-208b in that tissue. The tissues with high expression level are highlighted in red. However, hsa-mir-208b was not identified as tissue specific miRNA by Landgraf et al. [23]. It was reported that hsa-mir-208b show high heart-specificity by Liang et al. [24] and is associated with various cardiovascular diseases [12]. These results indicate that our prediction can provide us valuable assistance in identifying the tissue specific miRNAs.

Discussion
In summary, we presented a machine learning based approach to predict the co-expression patterns of the human intronic miRNAs and their host genes, which show a high accuracy and The combined data is described as Table 3. The symbol ''+'' means high co-expression and the symbol ''2'' means low coexpression. doi:10.1371/journal.pone.0004421.t004 Table 5. Tissue-specific miRNAs that are also found in our predictions. validation and could be further extended to other species. We further applied this approach to all the human intronic miRNAs and predicted the tissues in which the miRNAs are expressed for the miRNAs showing high co-expression pattern with their host genes. We also tested the validation of our method by predicting tissue specific miRNAs and performed a comparison with tissue specific miRNAs identified by biological experiments. As a result, our predictions are confirmed to be right. However, limitations exist in our current approach, which could be improved in the future study. We know that the process of miRNA and gene transcription is very complex and could be affected by lots of factors. Therefore, in order to improve the accuracy of the classifier, it will be very important to find better features that are associated with the expression of the intronic miRNAs and their host genes, for example the trans-action elements. For the machine learning methodology, currently we assigned the intronic miRNAhost gene pairs into high co-expression ones and low co-expression ones based on the correlation coefficients between the profiles of the intronic miRNAs and their host genes. This treatment could miss some detailed information, which may be useful in accurately predicting the correlations of the expression profiles between novel intronic miRNAs and their host genes. A regression model may provide more details in this situation. For the sample dataset, the current training dataset is small, which may affect the reliability of the classifier. Therefore, the prediction accuracy could be improved when more training samples become available. Furthermore, the current approach can only deal with the miRNAs that are within the introns of protein-coding genes. Developing computational approaches that can be used to predict the co-expression patterns for other kinds of miRNAs is also important. In a conclusion, our method can provide help in not only the understanding of miRNA expression and regulation but also the function of miRNAs, which can be further used to infer the miRNA-associated diseases.

Materials and Methods
The genome coordinate data of human miRNAs We downloaded the genome coordinate data of the human miRNAs from the miRBase [21] on August 2007 (gff-version 2). This database records 528 items of genome coordinates of human miRNAs. The genome coordinates were used to discriminate  whether one miRNA locates within the intron of a protein-coding gene or not.

The genome coordinates of human protein-coding genes
We obtained the genome coordinates data of known human RefSeq genes through the UCSC (University of California Santa Cruz)  [26]. Each item in this table lists one reference gene and its supplementary information, such as name, chromosome number, strand, transcription start position, transcription end position, coding region start position, coding region end position, number of exons, exon start positions, exon end positions etc. The genome coordinate data of protein coding genes together with that of miRNAs were used to determine whether a miRNA locates within the intron of a protein-coding gene.

Calculating the genome features of human miRNAs
The four feature vectors for the SVM classifier were calculated from the above genome coordinate data. We first mapped the genome coordinates of human miRNAs and protein-coding genes using Galaxy tool (http://g2.trac.bx.psu.edu/). We then calculated the lengths of all intronic miRNAs, the lengths of miRNAhosting introns, the lengths of miRNA-hosting genes, the distance from the transcription start position of one host gene to its intronic miRNA, and the distance of the transcription start position of one host gene to the host intron. These parameters will be used as the features of the samples (the miRNA-host pairs).

The dataset used in training and leave-one-out cross validation
Baskerville et al. performed a human miRNA microarray experiment across 30 normal tissues and investigated the coexpression patterns of 29 pairs of intronic miRNAs and host genes. Some pairs show high co-expression but some others show low coexpression. In the dataset, we set a cutoff of correlation coefficient at 0.6 to discriminate high co-expression and low co-expression. We also performed a leave-one-out cross validation on this dataset.

The testing dataset
In order to test the generalization power of our approach, we validated the method on an independent testing dataset containing miRNA expression profiles as well as mRNA expression profiles. It has been reported that microarray data from different platforms or laboratories often show great differences. Therefore, to improve the reliability and the persuasion of the results, we adopted two sets of miRNA microarray data. One miRNA microarray dataset was presented by Baskerville et al. [14], the other was obtained from Barad et al. [27]. We downloaded the host gene expression profile from Su's mRNA microarray data which are across 79 normal human tissues [22]. The intronic miRNAs in Baskerville et al.'s dataset and Barad et al.'s dataset were mapped to the host genes in Su's dataset, respectively. We then calculated the Pearson's correlation coefficients between the intronic miRNA expression profiles and their host gene profiles for both miRNA microarray datasets. Only those intronic miRNA-host gene pairs that showed consistent correlation patterns in both two miRNA microarray datasets was reserved for further analysis. As a result, we obtained 21 intronic miRNA-host gene pairs, in which there are 5 high co-expression pairs and 16 low co-expression pairs, respectively (Table 5).

Support vector machine
Support vector machine (SVM) was introduced by Vapnik [28], and has a comprehensive applications in many classification and regression problems. The goal of SVM is to construct a classifier from well-labeled data (the training data) that can be used to classify the incoming unlabeled data (the testing data). For a binary classification problem, for a given dataset, x i represents the feature vector, while y i represents the class labels (i = 1,2,…N, where N is the number of samples), where Here, for the intronic miRNA and host gene co-expression prediction problem, the input vector dimension is 5 (d = 5). The class label +1 represents the high co-expression class, while 21 represents the low co-expression class. Once the SVM was trained, it can be used to predict the class label (high or low co-expression) for a new sample (a new intronic miRNA and host gene pair). In this study, the LibSVM package [29] was used to train the SVM intronic miRNA and its host gene co-expression classifier and make predictions in the test dataset. We chose the Radial Basic Function (RBF) as the kernel function and tuned the parameters using the grid search strategy in LibSVM.

Software availability
We implemented our approach as a web server which is free for the scientific and technical community (http://cmbi.bjmu.edu.cn/ cepred/). The ''cepred'' software in linux and windows are also available for the scientific and technical community when requesting.

Supporting Information
Text S1 the predicted co-expression patterns of intronic miRNA with their host genes