Prediction of piRNAs using transposon interaction and a support vector machine

Abstract

Background

Piwi-interacting RNAs (piRNAs) are a class of small non-coding RNA primarily expressed in germ cells that can silence transposons at the post-transcriptional level. Accurate prediction of piRNAs remains a significant challenge.

Results

We developed a program for piRNA annotation (Piano) using piRNA-transposon interaction information. We downloaded 13,848 Drosophila piRNAs and 261,500 Drosophila transposons. The piRNAs were aligned to transposons with a maximum of three mismatches. Then, piRNA-transposon interactions were predicted by RNAplex. Triplet elements combining structure and sequence information were extracted from piRNA-transposon matching/pairing duplexes. A support vector machine (SVM) was used on these triplet elements to classify real and pseudo piRNAs, achieving 95.3 ± 0.33% accuracy and 96.0 ± 0.5% sensitivity. The SVM classifier can be used to correctly predict human, mouse and rat piRNAs, with overall accuracy of 90.6%. We used Piano to predict piRNAs for the rice stem borer, Chilo suppressalis, an important rice insect pest that causes huge yield loss. As a result, 82,639 piRNAs were predicted in C. suppressalis.

Conclusions

Piano demonstrates excellent piRNA prediction performance by using both structure and sequence features of transposon-piRNAs interactions. Piano is freely available to the academic community at http://ento.njau.edu.cn/Piano.html.

Background

Non-coding RNAs (ncRNAs) are important RNA molecules. Although they do not encode proteins, their roles in gene regulation are crucial [1],[2]. There are many types of long ncRNAs whose functions remain largely unknown [3]. Short ncRNAs, such as microRNAs (miRNAs) and piwi-interacting RNAs (piRNAs), are important post-transcriptional regulators [4]. piRNAs are produced from un-characterized precursors in both male and female germline cells. The discovery of piRNAs was a highly important break-through as they are involved in germ cell formation, germline stem cell maintenance, spermatogenesis and oogenesis [5]-[8].

The biogenesis of piRNAs is quite different from that of miRNAs. Although details of their biogenesis are currently unclear, several models have been proposed. In germline cells, piRNAs can be produced by the primary processing pathway and by a feed-forward loop, called the “ping-pong” pathway, which uses primary piRNAs to direct cleavage of complementary transposon sense transcripts [9]. These mature sense piRNAs will target complementary antisense piRNA precursors to create mature antisense piRNAs that can continue sense piRNA generation. piRNAs lack apparent structural motif and sequence conservation across different species, making their prediction a difficult task. piRNAs are generally understood to participate in transposon silencing during embryo development [10]. The majority of piRNAs are antisense to transposons. In the genome, piRNAs tend to occur in clusters and to be located in intergenic regions [5]. However, piRNAs are also found in somatic cells [11], and studying piRNA functionality is still a challenging task because of the wide variation of piRNA sequences.

piRNAs have been reported in human [12], mouse [6], rat [13], zebra fish [7], and fruit fly [14]. A typical experimental procedure to obtain piRNA data relies on immunoprecipitation of small RNAs bound to the protein PIWI and deep sequencing. However, with this method, it is still hard to identify piRNAs expressed at low levels or with restricted spatiotemporal expression. Therefore, computational prediction can provide an alternative approach to identify potential piRNAs. Unfortunately, homology sequence searching methods such as BLAST [15] or motif searching methods such as MEME [16] are not suitable for detecting piRNAs because sequence conservation is very low and no conserved structural motif has been detected in piRNAs.

The first de novo algorithm to identify piRNAs was a position-specific usage method that classifies piRNA sites along the genome using piRNAs starting with a uridine at their 5′ ends. A vector of 21 × 4 components was constructed containing 10 nucleotides upstream and 10 downstream of the starting U (i.e., +10 to −10, where U has the position of 0). The precision of this algorithm was only 61-72%, indicating that this tool is helpful for piRNA classification but still needs improvement [17]. Zhang et al. developed a k-mer based algorithm, named piRNApredictor, to predict piRNAs. piRNA and non-piRNA sequences from five model species were used as the training set. piRNApredictor has a high precision of >90% and a sensitivity of >60% [18]. piRNApredictor was integrated with mirTools 2.0 to predict piRNAs from small RNA-Seq data [19]. Moreover, iMir can be used to find piRNAs [20], but it mainly focuses on miRNAs. There is another program called "multiclass relevance units machine" that shows an excellent performance on piRNA classification [21]. However, it focuses on algorithm development and its software is not publicly available. proTRAC [22] and piClust [23] were developed to display known piRNA clusters, but they cannot be used to find new piRNAs.

Here, we present a new program, piRNA annotation (Piano), to predict piRNAs using piRNA-transposon interaction information. A support vector machine (SVM) was used to classify real piRNAs and pseudo piRNAs. Our analysis of Drosophila melanogaster data shows that Piano performs well in piRNA prediction, with over 90% prediction sensitivity, specificity and accuracy. The SVM classifier trained with Drosophila piRNA data can also accurately identify piRNAs of other species such as Homo sapiens, Mus musculus and Rattus norvegicus. Using small RNA-Seq data, Piano was successfully used to predict piRNAs for an important rice pest, the rice striped stem borer, Chilo suppressalis.

Methods

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.

Figure 1

Principle of triplet elements.

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

Prediction accuracy (ACC), specificity (Sp), precision (Pre) and sensitivity (Se) are widely used to evaluate the algorithm performance. The equations for these parameters are given below, with the following abbreviations: false positive (FP), true positive (TP), false negative (FN) and true negative (TN).

$$Se=\frac{TP}{TP+FN}\times 100\%$$

$$Sp=\frac{TN}{TN+FP}\times 100\%$$

$$Pre=\frac{TP}{TP+FP}\times 100\%$$

$$\mathrm{A}\mathrm{C}\mathrm{C}=\frac{TP+TN}{TP+TN+FP+FN}\times 100\%$$

Results and discussion

SVM classification

We used a SVM to classify real and pseudo piRNAs using 32-dimensional vectors of structure-sequence triplet elements. The training dataset was randomly divided into ten equally sized partitions. Each partition had the same ratio of positive samples to negative samples. Seven partitions were merged together as the training dataset. Two of the other partitions were merged together to validate the classifier for model selection. The tenth partition was used as the testing dataset. We used 10-fold crossing validation to improve the reliability. The training procedure was repeated ten times with different combinations of training set (seven partitions), validation set (two partitions) and testing set (one partition) (Table 1). We called our program that uses a SVM classifier with structure-sequence triplet elements to predict piRNAs, the piRNA annotation platform, abbreviated as Piano.

Table 1 Datasets for 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.

Figure 2

10-cross-validation results.

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|$$

Figure 3

The distribution of triplet elements in two datasets (pseudo piRNA vs. real piRNA).

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

To test the robustness of the program, we used the SVM classifier trained using the aforementioned Drosophila piRNA dataset to predict human, mouse and rat piRNAs. After aligning 32,152 human, 75,814 mouse and 66,758 rat piRNAs to relevant transposon sequences, 7,140 human, 14,495 mouse and 14,195 rat piRNAs were alignable and used in our cross-species application. The SVM classifier correctly recognized 6,690 out of 7,140 human (93.7%), 12,915 out of 14,495 mouse (89.1%) and 12,737 out of 14,195 rat piRNAs (89.7%). This gives an overall accuracy of 90.9% for the three cross-species datasets (Table 2).

Table 2 Cross-species validation results

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

Piano was compared with piRNApredictor, which was developed by Zhang et al. (2011). We used the same datasets to test the performance of these two methods. For each species, the testing data were composed of real piRNAs and pseudo piRNAs, all of which were mapped to the relevant transposon sequences (mismatch < =3). When predicting mouse piRNAs, compared with the algorithm proposed by Betel et al. (2007), piRNApredictor had high precision, 95.53%, and the sensitivity was 72.47% with the default parameter (t = 2). This means that piRNApredictor is good at recognizing positive but not negative samples. When comparing Piano and piRNApredictor with our datasets, Piano achieves higher sensitivity, specificity and accuracy than piRNApredictor (Table 3).

Table 3 Comparison between results from Piano and piRNApredictor

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.

Figure 4

The ROC curves of the two algorithms.

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].

Figure 5

The length distribution of piRNAs in five species ( D. melanogaster, C. suppressalis , R. norvegicus , M. musculus , and H. sapiens ).

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].

Figure 6

Percentages of piRNAs paired with different kinds of transposon in five species ( D. melanogaster , C. suppressalis , R. norvegicus , M. musculus , and H. sapiens ).

Conclusions

In this study, we developed a novel program for piRNA annotation called Piano. The program uses piRNA-transposon alignment/pairing and piRNA nucleotide content information (i.e., structure-sequence triplet elements) and achieves a high sensitivity, specificity and accuracy of over 90%. To the best of our knowledge, this is the best prediction performance achieved in comparison with other tools, such as piRNApredictor. Piano can be used not only for large-scale piRNA prediction from small RNA sequencing data but also for genome-wide annotation of piRNAs.