A bioinformatic-assisted workflow for genome-wide identification of ncRNAs

Abstract

With the upcoming of affordable Next-Generation Sequencing technologies, the number of known non-protein coding RNAs increased drastically in recent years. Different types of non-coding RNAs (ncRNAs) emerged as key players in the regulation of gene expression on the RNA-RNA, RNA-DNA as well as RNA-protein level, ranging from involvement in chromatin remodeling and transcription regulation to post-transcriptional modifications. Prediction of ncRNAs involves the use of several bioinformatics tools and can be a daunting task for researchers. This led to the development of analysis pipelines such as UClncR and lncpipe. However, these pipelines are limited to datasets from human, mouse, zebrafish or fruit fly and are not able to analyze RNA sequencing data from other organisms. In this study, we developed the analysis pipeline Pinc (Pipeline for prediction of ncRNA) as an enhanced tool to predict ncRNAs based on sequencing data by removing transcripts that show protein-coding potential. Additionally, a feature for differential expression analysis of annotated genes as well as for identification of novel ncRNAs is implemented. Pinc uses Nextflow as a framework and is built with robust and well-established analysis tools. This will allow researchers to utilize sequencing data from every organism in order to reliably identify ncRNAs.

Figure 1.

Overview on the process of generation of training data for the test runs. (A) Coding sequences (CDS; blue boxes) of H. sapiens, A. thaliana and S. cerevisiae were taken directly from RefSeq. Sufficient data for ncRNAs are available for human and thale cress; whereas the dataset of S. cerevisiae is too small to train CPAT solely on. Therefore, all RefSeq entries of ncRNAs of the phylum ascomycota were used as non-coding training data for CPAT. CDS were split. Right part: a portion of the CDS was used to contaminate the set of ncRNAs (white box) with coding RNAs in a ratio of 5:1 (blue, hatched box). This was done to simulate not annotated coding transcripts, which might still be in the dataset after filtering out all annotated CDS of the genome. CPC2 was used to predict ncRNAs within this mixed RNA pool. 80% of the sequences predicted as non-coding were used as non-coding training set (white box) for CPAT and 20% as the non-coding test set. Left part: Remaining CDS were split again: 80% for the coding training set, 20% for the coding test set. (B) The training datasets of ncRNAs (white and hatched boxes) CDS (blue box) were combined in a ratio of 1:1. 10-fold stratified cross-validation (CV) was used to calculate the model-specific ‘optimal’ cut-off. In each iteration the weighted Youden's index was used to calculate the cut-off. The mean of cut-offs from all 10 iterations is used to predict ncRNAs based on their coding probability calculated by CPAT.

Figure 2.

Overview of the process of generation of training data used in Pinc. (A) Right part: The filtered transfrags (white box) were subjected to CPC2 prediction of ncRNAs. 80% of the sequences predicted as non-coding were used as non-coding training set for CPAT and 20% as the non-coding test set. Left part: Coding sequences (CDS; blue boxes) are taken from the provided genome annotation and were split: 80% for the coding training set, 20% for the coding test set. (B) The training set that consists of ncRNAs (white box) was combined with the CDS training set (blue box) in a ratio of 1:1. 10-fold stratified cross-validation (CV) is used to calculate the model-specific, ‘optimal’ cut-off. In each iteration the weighted Youden's index was used to calculate the cut-off. The mean of cut-offs from all 10 iterations is used to predict ncRNAs based on their coding probability calculated by CPAT.

Figure 3.

Graphical overview on Pinc. Raw Sequencing reads are filtered based on quality and length using fastp. Subsequently, HISAT2 aligns the reads against the reference genome. StringTie assembles aligned reads into transfrags. Transfrags of already annotated features are removed by filtering for putative novel ncRNAs based on gffcompare's transfrag classification code. Together with the protein-coding RNAs from the reference annotation an organism-specific model is trained using CPC2 and CPAT to assess the coding probability of all putative, novel, non-coding transfrags. As edgeR requires the total count of reads mapped to each transfrag for a differential expression analysis, HTSeq-count was used to count the reads.