Computational Identification of piRNAs Using Features Based on RNA Sequence, Structure, Thermodynamic and Physicochemical Properties

1. INTRODUCTION

Since its discovery, RNA interference (RNAi) and its effector non-coding RNAs (ncRNAs) namely, microRNAs (miRNAs), small interfering RNAs (siRNAs), and PIWI-interacting RNA (piRNAs) have revolutionized our understanding of the mechanisms regulating gene expression [1, 2]. The common key event in all the RNAi regulatory mechanisms is the binding of small ncRNAs to Argonaute (AGO), forming a ribonucleoprotein complex termed as RNA-induced silencing complex (RISC) [3]. The AGO family of proteins is phylogenetically divided into two sub-families: somatic AGO [4], and germ line specific PIWI (P-element induced wimpy testis) clade [2, 5]. Both miRNAs and siRNAs act through the AGO family and constitute the RISC complexes known as the miRISC and siRISC, whereas the small ncRNAs interacting with PIWI are known as the piRNAs [6].

piRNA is a recently discovered class of ncRNAs, which are in the length range of ~24-32 nucleotides [1, 7, 8]. Initially, piRNAs were described as repeat-associated siRNAs (rasiRNAs) because of their origin from the repetitive elements such as transposable sequences of the genome [9]. However, later it was identified that they acted via PIWI-protein [10]. In addition to having a role in the suppression of genomic transposons, various roles of piRNAs have been recently reported like regulation of 3’ UTR of protein-coding genes via RNAi [11], transgenerational epigenetic inheritance to convey a memory of past transposon activity [12], and RNA-induced epigenetic silencing [13]. Furthermore, piRNA sequences are comparatively diverse than any other class of cellular ncRNAs and they constitute the most prevalent class of ncRNAs [7].

The overall mechanism of piRNA biogenesis is substantially different than that of the other ncRNAs. It includes siRNA, miRNA and long ncRNAs (lncRNAs) [13]. piRNAs are generated from long, single-stranded RNAs, which are transcribed from the genomic loci termed as piRNA clusters [14-19]. Earlier studies suggested multiple mechanisms of piRNA biogenesis in different cell types, tissues and organisms [10]. However, recent reports proposed a single and unified model that explained the mechanism of piRNA biogenesis across a range of evolutionarily diverse organisms [20-23]. These studies indicated that piRNA biogenesis could be divided into two parts: 1) the piRNA-dependent amplification loop, known as the ‘Ping-Pong cycle’, and 2) the piRNA-independent generation of phased trailing piRNAs. The former pathway begins with a maternally inherited 1U-biased piRNA, known as the initiator piRNA [24]. The PIWI-bound initiator piRNA cleaves the complementary single-stranded long transcript sequence into a pre-pre-piRNA with a terminal 5’ monophosphate [23]. The PIWI-bound pre-pre-piRNA undergoes subsequent RNA cleavage from its 5’ end to produce responder pre-piRNA. The responder pre-piRNA further undergoes 3’ end-processing to generate mature responder piRNA with 10A-bias. Subsequently, the mature responder piRNA enters the piRNA biogenesis cycle, acting as an initiator piRNA, and in turn, produces a new responder piRNA with 1U-bias, which is identical to the original initiator piRNA. The pathway operates as an amplification loop and hence, it is termed the ‘Ping-Pong’ cycle [7]. Thus, this arm of the unified model is dependent on a maternal initiator piRNA and also on the availability of long single-stranded piRNA precursor sequences. In contrast, phased piRNA generation is a piRNA-independent process, which operates through interaction with the mitochondrial endonuclease phospholipase D family member 6 (PDL6), an endonuclease present on the outer mitochondrial membrane. PDL6 cleaves the remaining 3’ end of the pre-pre-piRNA in a repeated manner, producing an array of tail-to-head pre-piRNAs, which are known as phased trailing piRNAs [21]. Thus, the current model proposed that the initiator piRNAs produce responder piRNAs and trailing piRNAs, which were known as the secondary and primary piRNAs, respectively, in the old model [22]. Hence, owing to the extensive role of piRNAs in regulating various biological processes, the conserved features of their biogenesis, and their sequence diversity, genome-wide identification of novel piRNAs would be of great importance to help understand piRNA-guided gene regulatory mechanisms across different cell types and organisms.

In piRNA identification, the general alignment-centered algorithms such as basic local alignment search tool (BLAST) [25] and Multiple Expectation Maximization for Motif Elicitation (MEME) [26] are well known for their robustness, but not for accuracy due to divergence of piRNAs across different species [27]. Recently, next-generation sequencing (NGS)-based methods have emerged as a powerful platform to identify piRNAs in a high-throughput mode [28]. However, in addition to piRNA, sequencing-generated data may also harbor reads from several other small ncRNAs such as miRNAs, endogenous siRNAs (endo siRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), smaller fragments of mRNAs and long ncRNAs in the piRNA length range. Therefore, the length range alone cannot be a characteristic parameter to screen for piRNAs in small RNA-based sequencing platforms. Hence, more accurate and robust algorithms are required to distinguish the real piRNA sequences from the pseudo piRNAs (non-piRNAs) of similar lengths by employing various discriminative features.

Recent studies employed multiple approaches to identify piRNAs with high accuracy in the genome of various organisms (Table S1 ^{(242.4KB, pdf)} ). Betel et al. proposed an algorithm based on the position-specific scoring matrix (PSSM) of 10 nucleotides upstream and downstream of piRNA sequences identified from the pachytene stage of mouse spermatogenesis, and predicted mouse-derived conserved piRNAs with ~60-70% precision [29]. Nevertheless, developed on a homogeneous dataset and trained on static position-specific nucleotide usage, the method shows position-specific features for only one species and hence, has limitations in piRNA identification in other species. Zhang et al. developed an algorithm using k-mer nucleotide frequency to capture the sequence-based characteristics in piRNA sequences [27]. Since the algorithm adopted dynamic nucleotide frequency calculation over the static PSSMs, it was able to identify piRNAs in the organisms beyond the training set species. However, there was a need to include additional features to improve the accuracy of piRNA identification apart from the sequence-based features. Other studies integrated additional features than sequence alone like structure and thermodynamic energy, and developed different piRNA prediction algorithms. The algorithm ‘PIANO’ was developed to predict piRNAs in Drosophila melanogaster using structural features and transposon interactions [30, 31], adopting local contiguous sequence-structure triplet-elements [27]. To predict mouse piRNAs using sequence motifs, ‘Pibomd’ was developed [32]. The above tools were trained on homogeneous datasets generated from single species (Table S1 ^{(242.4KB, pdf)} ). ‘Accurate piRNA prediction’ [33], ‘GA-WE’ (a genetic algorithm-based weighted ensemble method) [34] and ‘2L-piRNA’ (a two-layer ensemble classifier) [35] utilized sequence, structure and k-mer spectrum profile-based features to achieve better performance. Although these tools demonstrated a high accuracy level in predicting piRNAs, there was a need to develop a more robust algorithm capable of analyzing heterogeneous datasets from multiple organisms with high accuracy.

In this study, a comprehensive and robust machine learning-based algorithm is developed to predict piRNA sequences relying on the hybrid features such as spectrum profile or k-mer features (k=1 to 5), secondary structure, thermodynamic energy and physicochemical properties of RNA dinucleotides extracted from the piRNA sequences of eight species: Homo sapiens, Mus musculus, D. melanogaster, Caenorhabditis elegans, Danio rerio, Gallus gallus domesticus, Xenopus tropicalus, and Bombyx mori (Table S1 ^{(242.4KB, pdf)} ). During the 10-fold cross-validation (10-CV), our method reached an overall accuracy of >98% and a sensitivity of >98% in the above species. When compared with the existing algorithms, our hybrid predictive model piRNAPred demonstrated higher accuracy level. In conclusion, piRNAPred is the most updated piRNA identification method developed by integrating features of sequence, structure, thermodynamic energy and physicochemical properties extracted from a heterogeneous dataset across phylogenetically diverse organisms.

2. MATERIALS AND METHODS

3. Features utilized

3.1. Spectrum Profile or k-mer Features

One of the important contributing features of predictive model development is the k-mer string [33, 35, 38]. In machine learning, a k-mer string is defined as particular k-mer tuple (1, 2, 3, 4, 5 or even higher) of nucleotide or amino acid sequences that can be used to identify some representative motifs within DNA or proteins. The overall idea of using different k-mer strings is the identification of differential nucleotide usage between the authentic piRNA and pseudo piRNA sequences. However, when the dataset sequences are of different length, spectrum profile or k-mer nucleotide composition (k-MNC) is used [34, 39]. k-MNC is defined as a total number of a particular nucleotide divided by the length of the sequence. e.g., a piRNA sequence with five guanine residues and of length n, where n=19 nucleotides, the percentage composition of guanine residues will be 26.31%. Hence, there are 4 mononucleotide composition-based features (A, T/U, G and C).

Open in a new tab

(k, any amino acid or nucleotide)

Similarly, there will be 16-dinucleotide composition-based features (AA, AU, AC, AG, UA, UU, UC, UG, CA, CU, CC, CG, GA, GU, GC, GG). From k=1 to 5, they are designated as mono-, di-, tri-, tetra- and penta-nucleotide compositions. Hybrid predictive models (k=1, 2/1, 2, 3 etc.) involve higher-order combinations like hybrid of mono- and di-nucleotide composition is (md), mono-, di- and tri-nucleotide is mdt, mono-, di-, tri- and tetra-nucleotide is mdtt and mono-, di-, tri-, tetra- and penta-nucleotide is mdttp. Thus, k-MNC can be utilized to identify the motifs or regions of interest, which can distinguish the positive sequences from the negative ones in terms of differential nucleotide usage. Here, we employed a hybrid model of k=1-5, resulting into a total of 1364 strings of feature space.

4. Algorithm development

We used SVM^light software package (http://svmlight.joachims.org/) to train the 1516 features of D^1684p+1665n and to develop the predictive model in 10-CV mode. SVM^light is an implementation of a supervised MLT called support vector machines (SVM), which was proposed by Vapnik [44]. It works on the principle of pattern-recognition in the given dataset during training, and predicts the test-set based on these patterns. We used radial basis function (RBF) kernel having kernel width parameter (gamma, g) and regularization parameter (c) to build the SVM-based models in the range of 0.000001 to 5 and 0.00001 to 2000, respectively. Apart from SVM, we utilized other MLTs (Random Forest (RF)), Bagging, classification via regression, J48 pruned tree, naive Bayes and Instance Based Learner (Ibk)) using Weka package (Table 2).

5. Feature Selection to achieve dimensionality reduction

An optimal classifier depends upon features used for training. Therefore, to screen out the maximum relevant sub-space out of 1516 features, we employed attribute selected classifier in the Waikato Environment for Knowledge Analysis (Weka) package [45], and retrieved reduced feature space and provide 198 maximum relevance features (Table S4 ^{(242.4KB, pdf)} ). Further, we checked the performance of these 198 features in different MLTs (Table S5 ^{(242.4KB, pdf)} ).

6. Statistical assessment of algorithm

The performance of predictive models on the test set was assessed with the help of the statistical equation Mathews correlation coefficient (MCC). It is generally applied to calculate correlation between the actual and predicted values along with other statistical measures namely, percentage sensitivity (Sn), specificity (Sp), accuracy (Ac), true positive (TP), true negative (TN), false positive (FP), and false negative (FN) [35, 46]. These equations are provided below:

7. RESULTS

8. DISCUSSION

Cell growth, homeostasis and phenotypic expression of cellular functions have underlying complex regulations operating both at the gene expression and epigenetic levels [8, 10]. Gene expression regulation is majorly governed by small ncRNA-based surveillance and silencing in the biological systems [7, 47]. Among different players of RNAi, piRNAs are the most recently discovered and highly diverse class of ncRNAs [2]. They were first reported to govern the gonadal cell development by regulating the expression of transposons in the mouse germ cells [48-50], and loss of their expression leads to sterility [32-34]. Recent findings also suggest their roles in post-transcriptional gene regulation [11, 51], transgenerational epigenetic inheritance [12], regulation of mRNA decay [52], and localization in germplasm part of the embryo during embryogenesis [53]. Although many piRNAs have been recently discovered in some species, their identification has largely remained elusive in most of the organisms due to a lack of robust tools capable of identifying piRNA across phylogenetically diverse organisms with high accuracy.

The identification of piRNAs is majorly carried out by the NGS-based methods in which, many sequences from other ncRNAs fall within the range of piRNA sequence-length. Hence, to differentiate the piRNAs from the false-positive sequences, a more accurate and robust algorithm was warranted [33]. Recently, various methods were published which employed different features and used piRNAs

from different species. Integration of the k-mer spectrum profile in the algorithm was found to better distinguish the real piRNAs than the pseudo ones, and these k-mer nucleotide strings escalated to hybrid models by combining different k-mer nucleotide combinations [33]. Zhang et al., proposed a k-mer approach to identify piRNAs from five organisms i.e. rat, mouse, human, fruit fly and nematode [27]. Further, all the ncRNAs have a unique secondary structure, which can be computed by SSTE [31]. This approach was utilized to classify piRNAs in D. melanogaster [30]. However, this tool had a homogeneous training dataset, and this limited the discovery of novel piRNAs in the other non-related organisms [54]. ‘Pibomd’ extended the use of sequence-based motifs in the prediction of M. musculus piRNAs [32]. However, a comprehensive dataset was needed to capture all the possible sequence-based motifs and predict piRNAs from various evolutionary diverse organisms. Other algorithms like GA-WE [34] and 2L-piRNA [35] used pseudo dinucleotide composition (pseKNC), PSSM, mismatch profiles, etc., and achieved high accuracy of the classifier. However, the existing piRNA prediction algorithms did not combine various other features like secondary structure, thermodynamic energy and physicochemical properties of RNA, and did not utilize piRNA sequences from phylogenetically diverse organisms, which we considered crucial for the prediction of piRNAs with high accuracy.

In the present study, we collected 1684 positive piRNA sequences from eight different species (Table S1 ^{(242.4KB, pdf)} ). The length range of 1684 piRNAs was 24-33 nucleotides, and they were statistically distributed in 21 followed by 24, 29 and 31 nucleotides. To compare piRNAs and non-piRNAs, we investigated the composition-specific properties in detail and calculated the composition up to 5^th order of nucleotides that resulted in 1364 vector space. The accuracy of the classifier increased from k=1 to 5, suggesting that a combination of higher-order nucleotide usage (5-mer in our case) would better distinguish piRNAs from the non-piRNA sequences (Table 1). Apart from training the classifier on the differential nucleotide usage by incorporating k-MNC, we also trained our classifier on the relative position-specific occupancy of nucleotides at 1^st and 10^th position in the primary and secondary piRNAs, respectively, which resulted into a feature space of 16. However, their accuracy (57.12%) suggested that these features alone were not enough to distinguish the piRNAs from the non-piRNA sequences. Further, the high accuracy of the triplet element features suggested that sequence and secondary structure together could be a potent and discriminative property to segregate piRNAs from the pseudo sequences of similar length. However, there was moderate accuracy in the classifier trained on the properties based on Gibbs free energy and physicochemical values of dinucleotides. Importantly, the hybrid model piRNAPred which incorporated the k-mer spectrum profile, SSTE, thermodynamic properties of continuous dinucleotides and physicochemical features, demonstrated the highest accuracy among all the models with an MCC of 0.97.

CONCLUSION AND FUTURE IMPLICATIONS

In conclusion, the piRNAPred demonstrated the highest accuracy in predicting piRNAs as compared to the existing piRNA prediction tools. We hope it would be helpful in expanding our current understanding of piRNA biology by predicting novel piRNAs in different organisms. In the future, we will update more piRNAs from the other species, and develop a user-friendly interface. The scripts, predictive models, datasets and other supplementary material are provided on https://github.com/IshaMonga/piRNAPred.

SUPPLEMENTARY MATERIAL

Ethics Approval and Consent to Participate

Not applicable.

Human and Animal Rights

No animals/humans were used for studies that are the basis of this research.

Consent for Publication

Not applicable.

Availability of Data and Materials

All materials including computational workflow, predictive models, scripts and datasets utilized in this work are provided on GitHub repository at: https://github.com/IshaMonga/piRNAPred.

Funding

We acknowledge the Indian Institute of Science Education and Research, Mohali (IISER Mohali) for providing the financial support for this work.

CONFLICT OF INTEREST

The authors declare no conflict of interest, financial or otherwise.