Construction of microRNA functional families by a mixture model of position weight matrices

MicroRNA sequences and position weight matrix

We collected the mature miRNA sequences from miRBase (release 14) (Griffiths-Jones et al., 2008; Kozomara & Griffiths-Jones, 2011) and extracted the mammalian miRNAs from these. We selected for the analysis of the 6-mers at nucleotide positions 2–7, the seed regions, in the mature sequences. Then, we constructed the PWMs, which are scoring matrices weighted according to a specific position in the given seed sequences, with a computational learning approach. The PWMs for our experiments were initialized by randomly and repeatedly extracting sequences of six nucleotides from human genome sequences.

Mixture model of PWMs and the expectation–maximization (EM) algorithm

We developed a mixture model of position weight matrices to construct miRNA families, and estimated the model parameters with an EM algorithm to maximize the likelihood of the model.

Suppose that X = {X₁, X₂, …, X_i, …, X_N} is a dataset of N miRNA sequences, and each sequence length |X_i| is L. When the model is set as having k PWMs, the matrices are denoted by W = {W₁, W₂, …, W_N}, which is a set of 4 × L position weight matrices derived from miRNA seed sequences. The probability of sequence data X is then represented as: ${\displaystyle P\left(X\right|W)=\sum _{i=1}^k{\lambda }_{i}P(X|{\lambda }_i,W_i)=\sum _{i=1}^k{\lambda }_i\prod _{m=1}^N\prod _{v=1}^L{W}_i[u=X_{mv},v],}$ where λ_i is a weight value for the i-th PWM, and X_mv is the v-th base symbol in the m-th miRNA sequence. W_i[u, v] is the value of index (u, v) in the i-th PWM, and W_i[u = X_mv, v] is the probabilistic value matched to the base symbol of X_mv in the index (u, v) of the matrix W_i. The sum of λs should be 1, i.e.,: ${\displaystyle \sum _{i=1}^k{\lambda }_i=1({\lambda }_i\ge 0).}$ To show assignment of the data samples to the model, we introduce a hidden variable z_mi indicating the probability that an input sequence X_m is represented by a PWM W_i. Given these data, the expected log-likelihood to be maximized by the learning process in the PWM-mixture model is given by: ${\displaystyle logL(W,\lambda |X)=\sum _{m=1}^N\sum _{i=1}^{k}E(z_{mi})log(P\left(X_m\right|W_i\left){\lambda }_i\right),}$ L(⋅) is a likelihood function in the model, and E(⋅) is an expected value. The parameters of the model are estimated with the EM algorithm. The EM algorithm searches for the maximum value of the likelihood by iteratively repeating the E-step and M-step until convergence is achieved.

The E-step calculates the expected value of the hidden variable z_mi as follows: ${\displaystyle E\left(z_{mi}\right)=P(z_{mi}=1|X_m,W_i,{\lambda }_i)=\frac{P\left(X_m\right|W_i){\lambda }_i}{\sum _{j=1}^{k}P\left(X_m\right|W_j){\lambda }_j}.}$ E(z_mi) is the posterior probability that a miRNA seed sequence X_m is fitted in the position weight matrix W_i.

The M-step computes the parameters, W and λs, to maximize the log-likelihood by ${\displaystyle \widehat{{\lambda }_i}=\frac{1}{N}\sum _{m=1}^{N}E\left(z_{mi}\right),}$ ${\displaystyle W_i[\hat{u},v]=\frac{\sum _{m=1}^N{W}_i[u=X_{mv},v]E\left(z_{mi}\right)}{\sum _{b\in \{A,U,G,C\}}\sum _{m=1}^N{W}_i[u=b,v]E\left(z_{mi}\right)}.}$ Here, $\hat{\cdot }$ is the estimate of the parameters, and N is the sample size.

The algorithm is then run iteratively until the marginal likelihood is maximized. When the learning is finished, a miRNA sequence X_m is assigned to the i-th group among the k groups, according to the hidden variable z_mi that has the maximum value.

Score functions for validation

To interpret and validate our grouping results, we used two approaches to score the miRNA sequence sets in each cluster. The first scoring function, the match score, was originally developed to search transcription factor binding sites in DNA sequences (Kel et al., 2003). The function Score_M is calculated as follows: ${\displaystyle {\mathit{Score}}_M\left(X_M\right)=\frac{\sum _{v=1}^{L}I\left(v\right)f_{v,b}-\sum _{v=1}^{L}I\left(v\right)\underset{v}{\overset{\min }{f}}}{\sum _{v=1}^{L}I\left(v\right)\underset{v}{\overset{\max }{f}}-\sum _{v=1}^{L}I\left(v\right)\underset{v}{\overset{\min }{f}}},}$ where f_v,b is the frequency of nucleotide b (b∈{A, T, G, C}) at position v, and $f_v^{\min }$ is the frequency of the nucleotide that is the lowest frequency at position v, and $f_v^{\max }$ is the highest at the position. I(v) represents the conservation of position v in a matrix, defined as: ${\displaystyle I\left(v\right)=\sum _{b\in \{A,U,G,C\}}{f}_{v,b}log\left(4f_{v,b}\right).}$ The information value I(⋅) makes that mismatches in highly conserved regions are suppressed, whereas mismatches in less conserved regions are relatively acceptable.

The other score function is the silhouette measure, first described by Rousseeuw (1987). The score, Score_S, can be used to determine how tightly grouped all the datasets in the cluster are. The score function is defined as: ${\displaystyle {\mathit{Score}}_S\left(X_m\right)=\frac{\beta \left(X_m\right)-\alpha \left(X_m\right)}{max\left\{\alpha \right(\underset{m}{X}),\beta (X_m\left)\right\}},}$ where α(X_m) is the average distance of X_m to all the other sequences in the same cluster, and β(X_m) is the minimum value for the average distance of X_m to every other single cluster. The score value varies between −1 and 1. If Score_S(X_m) is close to 1, it means that the sequence X_m is well-matched to its own cluster, and dissimilar to the other neighboring clusters. The dissimilarity between two sequences is calculated as the Hamming distance.

Analysis of the functional relationships in each group

To evaluate the biological meaning of our results, it is necessary to analyze the functional relationships among the miRNAs in the same group. We predicted the target genes of the miRNAs within each group using microCosm Targets release version 5. We chose for further analyses the target genes that were bound by more than half the miRNAs within each group.

We investigated the biological process and molecular function categories of Gene Ontology (GO), and the entries for the KEGG pathways enriched with the target genes in each group, to verify the functional relationships among the human miRNAs assigned to the same group. These analyses were conducted using the DAVID Bioinformatics Resources (Huang, Sherman & Lempicki, 2009).

Next, we looked for the shared functions of the conserved miRNAs and their targets across species. We extracted information for homologues between human and mouse genes from the Ensembl database (Flicek et al., 2013), and then determined whether the miRNA members in the same group shared the conserved targets and their biological roles.

miRNA functional family construction

We ran our algorithm independently multiple times with random starts on the same datasets, and with a varying number of clusters (5–100). For each number of clusters, we repeated the algorithm 10 times, because of the possible existence of many local maxima for the mixture model, and we selected the one that maximized the likelihood value. With the highest likelihood value, the chosen number of clusters was 81. All the grouping results are shown in Table S1.

We computed two score measures, Score_S (silhouette score) and Score_M (match score), to validate our grouping results (see Methods section). We calculated the average values for each group, and then compared them with the random and hierarchical clustering results. To confirm that our approach identified relevant miRNA families, we also compared the results with the family information in miRBase (Griffiths-Jones et al., 2008; Kozomara & Griffiths-Jones, 2011). In miRBase, the miRNA families have already been defined, and several of them are known to be sequentially and functionally conserved. We compared our family construction results with the miRBase family information, random groupings and hierarchical clustering results using the average Score_S and Score_M values. Table 1 shows that the Score_S and Score_M values for our results are high. The average Score_M for our results is similar to that for the miRBase families, and the Score_S is much better in our results. The match score, Score_M simply evaluates how well the sequences are represented by the PWM in their own group, and does not explain how well the similar data are collected together by dividing the incompatible instances into several groups, because the score does not measure the differences with samples in other groups. Therefore, the scores for our results indicate that our approach can be used as a way to construct a sequence-based family, by assigning similar miRNA sequences to an identical group well, while assigning dissimilar sequences to different sets.

Table 2 shows that the well-known miRNA family members defined in miRBase are grouped together in our experiment. For example, all the let-7 miRNA sequences, members of one of the most well-known families, are collected into the same group, cluster 10. The let-7 miRNA family members are known to be highly conserved across species in both their sequences and functions, and the members play roles in tumor suppression, and cell differentiation, proliferation and development (Roush & Slack, 2008). The mir-181 family is grouped into cluster 44 in our experiment. Its members are considered to be oncogenic miRNAs that down-regulate the Tcl1 overexpressed in mature B-cell lymphomas, and Hox protein, a repressor of differentiation processes in mammals (Pekarsky et al., 2006; Naguibneva et al., 2006). The members of the miRNA families that show corresponding results, including let-7, mir-15, mir-181, mir-196, and so on, share similar sequences, and the sequence conservation in our groupings can also be detected with WebLogo (Fig. S1) (Crooks et al., 2004). Our results confirm that our approach effectively groups previously known miRNA family sequences together.

Although many of our groupings are similar to previously constructed miRNA families as shown in Table 2, several results differ. For example, we cannot find sufficiently similar sequence patterns in the sequences of miRBase MIPF0000018 family (mir-154) (Fig. 1A, Table 3). In our experiment, the miRNAs in MIPF0000018 were allocated to several different groups, including clusters 3, 37, and 52, and the similarities between the sequences in our groups were shown more clearly. For further confirmation, we calculated the silhouette measure (Score_S) and match score (Score_M) for MIPF0000018 (Table 3). The Score_S and Score_M for the miRBase family were markedly lower than the results for our groups.

In more detail, the miR-381 sequences are included in the MIPF0000018 family with the mir-154 miRNAs. However, in our analysis, the mir-381 sequences are grouped with the mir-466-3p sequences in cluster 52. In this group, the average Score_S is 0.767 and the average Score_M is 0.990, which means that the sequences in the group are highly conserved. These scores show that the mir-381 sequences are sufficiently similar to those mir-466-3p, as expected. The sequence similarities are clearly shown using WebLogo in Fig. 1B. The first five bases of the seed sequences are identical, AUACA in the sequences. From this analysis, it might be supposed that mir-381 was inherited from the same ancestral miRNA gene as mir-466-3p, but not the same as that of mir-154. The inference is reasonable because the mature form of the mir-381 miRNAs is also a 3′-donor sequence in their secondary structure, like mir-466-3p. The property that the miRNA members within one miRBase family are split into several groups in our results is also found in other miRNAs including mir-188 and mir-506 (Table S1). Conversely, MIPF0000013 (mir-25) and MIPF0000069 (mir-32) in miRBase are merged into one group, cluster 62, in our results. Actually, the seed sequences of mir-32 are completely identical to those of mir-25, AUUGCA. Our grouping of cluster 62 is strongly supported by the observation that mir-32 and mir-25 have similar roles. For instance, these miRNAs can lead to cancer by inhibiting apoptosis, because mir-32 suppresses the expression of Bim protein, a pro-apoptotic factor, as mir-25 does in cancer cells (Petrocca et al., 2008; Ambs et al., 2008). Therefore, we established that our method might identify functionally related families.

Functional analysis for the constructed family

To verify the biological functions of our miRNA groups more comprehensively, we conducted the Gene Ontology (GO) and KEGG pathway enrichment analysis of the human target genes of the miRNAs within each group. The functions of miRNAs are strongly related to the biological roles of their target genes, because the miRNAs recognize the target mRNAs and inhibit their expression. We used the target information produced in microCosm Target version 5. The enrichment results for the target genes in each group are assessed in Tables S2 and S3. We show the most significant part of the results of the GO enrichment analysis using biological process terms in Table 4, and for the KEGG pathways in Table 5. In most of the groups, the target genes are significantly involved in several GO categories. Similarly, we checked the biological relationships among the miRNAs in our groups by the KEGG pathway enrichment analysis, and found frequently and statistically enriched pathways based on the target genes of the members in each group. These results for the GO annotation and KEGG pathway enrichment analyses are consistent with the biological functions of the miRNAs previously reported in the literature. For example, the target genes of cluster 63 were related to the cell cycle in our enrichment analyses. In fact, it is already known that the mir-16 gene family in cluster 63 regulates cell-cycle progression (Linsley et al., 2007; Xia et al., 2009).

Finally, we checked the conserved target information for the human and mouse miRNAs categorized into the same group. Figure 2A shows an example of the conserved target analyses in the two different species. The human and mouse miRNAs within cluster 63 share orthologous target genes. These analyses help to clarify the functions of the miRNAs in each group because the target relationships are conserved across most species. Using our miRNA grouping results, we can speculate about their functions based on our knowledge of other miRNAs, with already-known functions, in the same group. As an example, it has previously been shown that the expression of Bcl2 protein is negatively regulated by miR-15 and miR-16, which are assigned to the same group, cluster 63 (Cimmino et al., 2005). Because the Bcl2 protein inhibits cell apoptosis, its overexpression leads to leukemias or other cancers (Cory & Adams, 2002). In our experiment, there were many other miRNAs within the cluster 63, such as mir-322, mir-424, mir-497, and so on, grouped with the miR-15 family members, whose functions are not yet clearly known. From our results, we can infer that these miRNAs may induce cancers by a similar mechanism to that of miR-15 and miR-16 (Fig. 2B). In practice, microCosm predicts that mir-503 binds to the Bcl2 gene transcript. Therefore, our grouping analyses, together with the conservation information, provide clues to the biological effects of functionally unknown miRNAs. Furthermore, the example demonstrates that, as well as identifying the miRNA functions, our approach can help to discover miRNA-mRNA modules in the complex gene regulatory networks and to understand the combinatorial effects of miRNAs in cellular processes.