Prediction of piRNAs using transposon interaction and a support vector machine

Training and testing sets

Two datasets were built for D. melanogaster: one contained real piRNAs and the other contained pseudo piRNAs. We downloaded 987 piRNAs from the NCBI GenBank database (GI: 157361675–157362817) [24] and 12,903 piRNAs from the NCBI Gene Expression Omnibus with the accession number GSE9138 [14]. By using short sequence alignment software, SeqMap [25], highly similar sequences were removed. After removing redundancy, 13,848 non-redundant piRNAs were kept. We downloaded 261,500 Drosophila transposons from the UCSC Genome Browser (Apr. 2006 dm3) [26]. We aligned 13,848 piRNAs to the transposon sequences using SeqMap with a maximum of three mismatches allowed. Among 13,848 non-redundant piRNAs, 9,758 (70.4%) could be aligned successfully, suggesting that they can target transposons.

Since DNA sequences are not random sequences, there are some differences between coding and non-coding RNAs. Because piRNAs are non-coding RNAs, we used non-coding RNAs as a negative control to generate our pseudo piRNA dataset. We downloaded 102,655 Drosophila ncRNA sequences from the NONCODE v3.0 database [27]. First, we removed all piRNAs from this dataset. We then randomly selected one ncRNA sequence and randomly cut out a short sequence of 20–30 nt as one candidate sequence. By this double-randomization process, we were able to obtain about 200,000 candidate pseudo piRNAs. Next, we mapped all these candidate sequences to the transposons with a maximum of three mismatches, and those sequences that did not map to the transposons were removed from the candidate sequence dataset. Accordingly, we produced 38,919 non-redundant candidate pseudo piRNAs. We then randomly selected some candidate pseudo piRNA sequences to simulate the length distribution of real piRNAs. Finally, we obtained 9,240 sequences that formed the pseudo piRNA dataset as the negative dataset for SVM classification.

Cross-species test set

We applied the SVM classifier trained with Drosophila piRNAs to human, mouse and rat data. In total, 32,152 human, 75,814 mouse and 66,758 rat piRNAs were downloaded from the NONCODE v3.0 database [27]. Transposons of the three species were downloaded from the UCSC Genome Browser [26], including 8,537,572 human, 7,320,714 mouse and 6,380,192 rat transposons.

Structure-sequence triplet elements

The main function of piRNAs is to silence transposons. To target transposons, piRNAs need to bind with their target sequences. In piRNApredictor [18], 1,364 k-mer strings (k = 1, 2, 3, 4, 5) were used to describe piRNA sequences. Although this k-mer approach is a good way to characterize and extract sequence content features from piRNAs, it is purely a mathematical method that might lack biological insight and significance. In our program, we analyzed piRNA-transposon interaction information using RNAplex [28]. When a piRNA binds with a transposon, there are two statuses for each nucleotide of the piRNA, paired or unpaired (see Figure 1). The paired nucleotides of piRNAs are indicated by opening brackets "(" whereas the paired nucleotides of transposons are indicated by closing brackets ")". The unpaired nucleotides in both piRNAs and transposons are indicated by dots “.”. For any three adjacent nucleotides of a piRNA, there are 8 (2³) possible structure compositions: "(((", "((.", "(..", "(.(", ".((", ".(.", "..(" and "…". Combining the middle nucleotide of each three adjacent nucleotides, we can form 32 (4 × 8) different triplet elements that contain both structural information of transposon-piRNA alignment/pairing and piRNA sequence information, which we call structure-sequence triplet elements (see Figure 1). For instance, the triplet element “(((U” indicates that three contiguous piRNA nucleotides are aligned with a transposon and the middle nucleotide is U.

Support vector machine

Support vector machines (SVMs) have been widely applied in the classification of biological signals. For a given dataset, x _i ∈ R _n (i = 1,…N) with corresponding labels y _i (y _i  = +1 or −1, representing real and pseudo piRNAs respectively in this work), SVM gives a decision function ${\displaystyle \int }\left(x\right)=sgn\left({\displaystyle {\sum }_{i=1}^N}{y}_i{\alpha }_{i}K\left(x,x_i\right)+b\right)$, where α _i represents the coefficients to be learnt and K is the kernel function. The LibSVM3.12 package (http://www.csie.ntu.edu.tw/~cjlin/libsvm/) [29] was used to perform the analysis. For optimizing the SVM classifier, the penalty parameter C and the RBF kernel parameter γ were adjusted using the grid search strategy in LibSVM.

Prediction system assessment

SVM classification

In one of these tests, Piano correctly recognized 935 out of 975 real Drosophila piRNAs, and detected 874 out of 924 pseudo piRNAs as negative cases (Additional file 1: Table S2). We calculated the average value of ten tests. Piano gives a sensitivity of 95.89 ± 0.50%, specificity of 94.61 ± 0.81%, accuracy of 95.27 ± 0.34%, and precision of 94.95 ± 0.71% (Figure 2). The high performance of Piano indicated that real and pseudo piRNAs are quite different in terms of structure-sequence triplet elements. The triplet elements combine both structural information of piRNA-transposon alignment/pairing and sequence information of the middle nucleotide of three contiguous piRNA nucleotides. Such a structure-sequence triplet element was previously used to classify real and pseudo miRNAs [30], suggesting that this structure-sequence feature might be common for small ncRNAs. Although pseudo piRNAs are also antisense to transposons due to their alignment (see Methods), they can be effectively distinguished from real piRNAs by the triplet elements, demonstrating that piRNA-transposon interaction information is an intrinsic characteristic of piRNAs.

We calculated the average frequencies of the 32 structure-sequence triplet elements in the real piRNAs and pseudo piRNAs. Our data analysis indicated that "(((G" and "(((C" appear at higher frequencies in real piRNAs than in pseudo piRNAs. The group of two-paired nucleotides and one unpaired (e.g., "((.A") appears more often in pseudo piRNAs than in real piRNAs (Figure 3). We calculated the F-value to estimate the discriminative power of the different triplet elements [31],[32].

$F\left(x_j\right)=\left|\frac{{\mu }_j^{+}-{\mu }_j^{-}}{{\sigma }_j^{+}+{\sigma }_j^{-}}\right|$

For each feature x _j , j = 1, …, N, we calculated the mean ${\mu }_j^{+}$ (${\mu }_j^{-}$) and standard deviation ${\sigma }_j^{+}$ (${\sigma }_j^{-}$) using positive or negative examples, respectively. The results demonstrated that “…G”, “(.(G”, “..(C”, “..(G”, and “(..C” are the top five discriminative elements. Four of them contain continuously unpaired nucleotides, suggesting that binding stability between piRNA-transposon interactions is the key information in classifying real and pseudo piRNAs (Additional file 2: Table S1).

Application of Piano to other species

The high accuracy in predicting mammalian piRNAs achieved by the SVM classifier trained with Drosophila piRNAs suggests that the structure-sequence triplet element represents a conserved feature for piRNAs.

Comparison with other methods

As shown in Table 3, using the same datasets for the three species, Piano has prediction specificity of ~40%, which is much higher than that of piRNApredictor (~10%). Figure 4 shows the ROC curves (AUC) of piRNApredictor and Piano. AUC is a global performance measure because it integrates overall threshold values [33]. Clearly, Piano achieves better performance than piRNApredictor in identifying piRNAs.

Prediction of Chilo suppressalispiRNAs

Rice striped stem borer (SSB) is an important rice pest that causes huge yield loss. To date, no piRNAs have been reported in SSB. We applied our program to predict piRNAs from small RNA-Seq data; 2,170,655 short sequences in total. From this data, 82,639 piRNAs were predicted. The whole prediction procedure takes ~7 hours on an Ubuntu server (Sugon X8DT6, 2 CPU processors, each has 12 threads, 48 G memory). An interesting discovery is that insect piRNAs might have a different length distribution than mammalian piRNAs. The mammalian piRNAs have a length peak at 29–30 nt, whereas that in Drosophila is 24–26 nt and that in SSB is 27–28 nt (Figure 5). These findings are consistent with previous results [14].

piRNA target sequences

The main function of piRNAs is to target and silence transposons. In this study, we analyzed piRNAs and their target sequences in human, rat, mouse, fruit fly and rice stem borer. We calculated the percentage of piRNAs targeting different categories of transposons. Our data analysis indicated that the majority of human piRNAs (95.0%) target SINE transposons. In mouse, 67.5% of piRNAs target SINE and 24.9% target LINE transposons. In rat, 65.6% of piRNAs target SINE and 29.0% target LINE transposons. In Drosophila, 66.8% of piRNAs target LINE and 26.4% target LTR transposons. In SSB, 42.4% of piRNAs target LINE and 44.0% target SINE transposons (Figure 6). These results indicate that piRNAs may have somewhat different mechanisms of action in different species [34],[35].