3.1 The characteristics of the identified interactions in CLASH
We identified 12,170 miRNA-mRNA and 58,043 isomiR-mRNA interactions from the six CLASH samples (Supplementary Table S2). We observed each of the eight types of isomiRs, with a 3’ addition isomiR hsa-miR-4268 occurring the most frequently in 3,565 interactions, while 96 isomiRs occurring only 10 times (Supplementary Figure S1). Moreover, 98 isomiRs were found to participate in at least 100 interactions. Consistent with the previous studies [3], there were more 3’ isomiRs than other types in the CLASH data set (Supplementary Table S2). The number of isomiRs with nucleotide addition were more than other types. The number of SNP and MNP isomiRs were similar in the dataset. Despite the varied frequency of different types, there were at least 9 isomiRs from each of the eight types (Supplementary Table S3). Note that, these statistics were based on all the documented miRNAs in miRbase database. We also cross-referenced the miRNAs with curated miRNA database miRGeneDB [39] and found that 66 of the documented 268 miRNAs in miRGeneDB were in the CLASH dataset. Of them the most frequent isomiR was a 3’ addition isomiR of hsa-miR-615-3p that occurred in 824 interactions. The number of 3’ isomiRs was still higher (218) than 5’ isomiRs (57) and polymorphic isomiRs (10).
Since we considered isomiRs covered by at least 10 reads, all de novo identified isomiRs were supported by their recurrent occurrence. The mean and median number of reads supporting these isomiRs were 53 and 19, respectively. There were 200 miRNAs with at least one identified isomiR. The number of identified isomiRs for a miRNA varied from 1 to 98. On average, a miRNA had around 6 isomiRs of different types. Despite the existence of different types of isomiRs, it was evident that one miRNA preferred specific types of isomiRs. In other words, for a given miRNA, a specific type of isomiRs occurred much more frequently. In fact, for all 67 miRNAs with at least 100 isomiR-mRNA interactions, at least one type of isomiRs occurred in significantly higher than expected frequencies (corrected Binomial test p-value < 0.01).
The identified isomiRs were likely to be condition-specific. We compared the identified isomiRs in the CLASH interactions with the isomiRs in the CLEAR-CLIP interactions. The CLASH samples were from a healthy kidney cell line, while the CLEAR-CLIP samples were from a carcinomic liver cell line. Among the 1,226 isomiRs and exact miRNAs identified in the CLASH data, 1,203 (98.12%) were not identified in the CLEAR-CLIP data. If we considered the 1,095 isomiRs and exact miRNAs supported by at least 50 reads in the CLASH data, 1,078 (98.45%) were still not identified in the CLEAR-CLIP data. It was thus not the isomiR quality that made the difference of isomiRs in different experiments. In other words, isomiRs and their interacting sites are likely to be condition-specific.
We investigated the difference between the unique target sites of the 5’ and 3’isomiRs and those of the exact miRNAs in the CLASH dataset. Among the 5,742 CLASH target sites of exact miRNAs, 355 were common with the 5’ isomiRs’ while 2,353 were common with the 3’ isomiRs’. The lower number of common targets of exact miRNAs with the 5’ than the 3’ isomiRs corroborates that the 5’ isomiRs have altered target specificity. Interestingly, the 3’ isomiRs targeted much more unique mRNAs (5,021) than other types, indicating that the 3’ miRNA regions may also be important for miRNA/isomiR targeting.
3.2 DMISO had good performance on miRNA and isomiR target prediction
We evaluated DMISO by 10-fold cross-validation on the training data (Supplementary Table S4). It showed over 99% area under the receiver operating characteristic curve (AUROC) and area under the precision-recall curve (AUPR), and more than 93% F1 scores, precision, and recall.
We also evaluated DMISO on three independent datasets: the 20% left-out CLASH test data, the CLEAR-CLIP data, and the miRTarBase data (Table 1). The performance of DMISO on the left-out CLASH test data was similar to that of the above cross-validation on the CLASH training data. That is, the AUROC and AUPR were over 99%, and the F1 score, precision and recall were more than 93%. When DMISO was tested on the CLEAR-CLIP data, the performance was slightly lower (94% AUROC, 99% AUPR, 94% F1, 98% precision and 90% recall). Because we did not have negative pairs in the miRTarBase dataset, we could only evaluate the recall of DMISO on this independent dataset (Table 2). DMISO had a recall of 92%, almost as good as the recall on the CLEAR-CLIP and 20% left-out CLASH test data.
Table 1
Performance comparison on CLASH 20% test data and CLEAR-CLIP data.
Data set
|
Tools
|
Pos
|
Neg
|
AUROC
|
AUPR
|
F1
|
Precision
|
Recall
|
Specificity
|
CLASH test
|
DMISO
|
14053
|
8319
|
0.9945
|
0.9969
|
0.9581
|
0.9286
|
0.9895
|
0.8714
|
miRanda
|
14053
|
8319
|
0.5873
|
0.6921
|
0.3032
|
0.9851
|
0.1792
|
0.9954
|
RNA22
|
14053
|
8319
|
0.5002
|
0.6283
|
0.0023
|
0.7273
|
0.0011
|
0.9993
|
TargetScan
|
14053
|
8319
|
0.5568
|
0.6686
|
0.2181
|
0.9563
|
0.1231
|
0.9905
|
miRAW
|
14053
|
8319
|
0.6119
|
0.6889
|
0.6552
|
0.7305
|
0.5940
|
0.6298
|
miTAR
|
14053
|
8319
|
0.5970
|
0.6832
|
0.5300
|
0.7639
|
0.4057
|
0.7882
|
CLEAR-CLIP
|
DMISO
|
14684
|
1323
|
0.9427
|
0.9945
|
0.9398
|
0.9844
|
0.8990
|
0.8420
|
miRanda
|
14684
|
1323
|
0.5230
|
0.9211
|
0.1135
|
0.9790
|
0.0603
|
0.9856
|
RNA22
|
14684
|
1323
|
0.4993
|
0.9173
|
0.0004
|
0.6000
|
0.0002
|
0.9985
|
TargetScan
|
14684
|
1323
|
0.6001
|
0.9337
|
0.3670
|
0.9901
|
0.2252
|
0.9751
|
miRAW
|
14684
|
1323
|
0.5704
|
0.9283
|
0.7030
|
0.9367
|
0.5626
|
0.5782
|
miTAR
|
14684
|
1323
|
0.6904
|
0.9482
|
0.6596
|
0.9793
|
0.4972
|
0.8836
|
Table 2
Performance comparison on miRTarBase data.
Tools
|
Recall (Overall)
|
Recall (Functional)
|
Recall (Non-functional)
|
DMISO
|
0.9164
|
0.9168
|
0.8938
|
miRanda
|
0.7045
|
0.7054
|
0.3540
|
RNA22
|
0.0199
|
0.0195
|
0.0221
|
TargetScan
|
0.7632
|
0.7643
|
0.3053
|
miRAW
|
0.6734
|
0.6696
|
0.5531
|
miTAR
|
0.0090
|
0.0090
|
0.0044
|
The above analysis was on all interactions in the three test datasets. We further studied how well DMISO predicted the interactions involving different isomiR types instead of exact or wild-type miRNAs (Supplementary Table S5). DMISO showed consistently good performances on different types of isomiR-mRNA interactions. For instance, DMISO had an AUROC of 94%, an AUPR of 99%, a F1 score of 93%, a precision of 98%, and a recall of 89% on the CLEAR-CLIP data for the 5’ isomiR-mRNA interactions. Overall, although the performance was close for different types of isomiR-mRNA interactions, it was slightly better on 3’ isomiR-mRNA interactions on both CLASH and CLEAR-CLIP datasets.
We also compared DMISO with five other tools: miRanda, RNA22, TargetScan, miRAW, and miTAR. The input to every tool was the pairs of positive or negative miRNA/isomiR-mRNA target sites. DMISO showed a superior performance to all five tools in terms of every metric considered (Tables 1 and 2). For instance, DMISO had an AUROC of 94% and an AUPR of 99% on the CLEAR-CLIP data, compared with the highest AUROC of 69% and the highest AUPR of 95% by the other five tools. The five existing tools had a smaller precision and a much smaller recall than DMISO, which may be because DMISO was the first tool that took the isomiR-mRNA interactions into account to train the models. It also highlighted the importance of considering such interactions for future miRNA target predictions.
Among the other five tools, the deep learning-based miRAW and miTAR had slightly higher AUROC and AUPR than those of the three classical tools (TargetScan, miRanda and RNA22) and much higher recall scores. This indicates that miRAW was able to capture the information in the non-seed regions, unlike the three tools that focused on the seed regions. This implied the importance of the non-seed regions for miRNA target site identification [19, 20, 40]. On the miRTarBase dataset, DMISO had a recall that was at least 10% larger than the recall scores of other tools. In contrast, miRAW and miTar had a worse recall than miRanda and TargetScan, suggesting that the deep learning models used in the two tools might not consider certain well-studied features of miRNA target interaction captured by the DMISO model (Table 2).
3.3 Features analyzed
Despite the high accuracy of deep learning models to solve a problem, these models are infamous for being unable to interpret what they learn from the data. To address this problem, we used two feature interpretation methods, convolutional kernel analysis and input perturbation [22], as an attempt to unwind the learning process of the DMISO model.
Since DMISO has convolutional layers, we applied the convolutional kernel analysis method to extract the higher-level features learned by the model. In this process, the 10 kernel matrices, each of size 4 x 8, of the two convolutional layers for the miRNA (isomiR) and target branches were analyzed to find the weights of the A, T, C and G learned by the layers. In the beginning of the training process, the weight values of the kernels were randomly initialized in the forward pass and updated through backpropagation of the loss. After 500 epochs of training, 3 of the 10 kernel matrices in the miRNA (isomiR) branch and 7 of the 10 kernel matrices in the mRNA branch picked up certain weight values that might be commensurate with the protein binding motifs in the miRNA and mRNA sequences. For the remaining kernels in the respective branches, the weight values were too similar (25%) to be assigned to one of the 4 nucleotides. To find out the binding protein motifs, we compared the 4 x 8 kernel matrices of both the miRNA and target branches with the JASPAR motif database on vertebrates [41]. The top motifs that each kernel matched with high significance were the protein binding motifs GATA1::TAL1, ZFP42, RARA::RXRG, RARA::RXRA, ESR2, ZFP42, ZBTB26, etc., all of which are CCHH and CCCH type Zinc finger proteins. Zinc finger proteins are well known as RNA binding proteins, which are essential to bind with ribonucleoproteins in a RNA-induced silencing complex [42]. This analysis shows that DMISO was able to recognize the binding profile of RNA binding motifs in miRNA/isomiR and mRNA sequences through the kernels of the convolutional layers, which has not been considered in the existing miRNA target prediction tools.
The other feature interpretation method that we used was the input perturbation technique. In this case, for every miRNA and target sequence in a dataset, a 4-nucleotide long mask consisting of “N” was applied on the sequence starting at every position of the miRNA and the target sequences. The changed value of DMISO’s prediction probability after applying the mask was recorded for the region in the corresponding sequences. The mask was then slid across the miRNA and the target sequences by one nucleotide every time. The average changes in the prediction probability after scanning all the miRNA and target sequences within a dataset gave us the regions that were most significant to DMISO within all the miRNA and target sequences.
The input perturbation method was applied to the 20% CLASH test data and CLEAR-CLIP datasets (Fig. 2). For each position, the mean (middle blue) and variation (gray area) of the sensitivity of DMISO to the changes at that position for all miRNA/isomiR and mRNA sequences in a dataset was recorded. In both datasets, positions 1–9 of the miRNAs/isomiRs had the most variation, which confirmed the importance of the seed region. However, not all positions in the seed were of the same importance. For instance, the first position was associated with the lowest variation while positions 3–5 were the highest in the seed for both datasets. The sensitivity variation of DMISO was also similar for both datasets from position 1 until position 18, and then dramatically decreased to the end. This might be because of the difference in miRNA/isomiR lengths (the shortest miRNA/isomiR was 17 nt long). It also suggested the importance of almost all miRNA positions instead of only the seed for target binding. When the target mRNA sequences were changed, DMISO reacted more on the modifications in the 5’ regions of the target sequences, which further implied the importance of the 3’ regions of miRNAs/isomiRs that bound the 5’ region of the target mRNA sequences. In addition to the high sensitivity variation in the 5’ region, the positions around 50 (3’ region) of the target sequences showed a spike in sensitivity variation in both datasets. Since the 3’ region of the target corresponded to the seed region of the miRNA/isomiR, this confirmed the importance of the high-quality match in this seed region as well.
We also clustered miRNAs based on the binding sensitivity scores at every position of miRNAs. Our intuition was, if the model showed similar reaction patterns to the changes to two miRNAs, the two miRNAs might have similar binding patterns or features (Supplementary Figure S2). We found that most of these miRNA clusters had many common targets, and a cluster of miRNAs often work in the same pathways (Supplementary Table S6).