How to get the rRNA ratio from a RNAseq dataset

Question

[niu2rseq]

Hello,

I want to know if there is any way using the bedtools and miRdeep2 output bed file to get the rRNA ratio in my miRNAseq fastq data. Thank you very much!

I have a gtf file, a genome.fa, a bed file from the miRdeep2. Thanks!

[Devon Ryan]

Why don't you just align to an rRNA fasta file and then count the reads? That'd be quick and easy.

[niu2rseq]

Thank you for this idea! Here is what I did:

I downloaded the sequence M_ncrna.fa from Ensembl.

And from this I pulled out all rRNA sequences here:

awk '/rRNA/{print $0;getline;print $0}' M_ncrna.fa > M_rRNA.fa

and then I build a bowtie2 index for this M_rRNA.fa

and then I did:

bowtie2 M_rRNA my.fastq -S aligned.sam

I don't really care about the sam file but from the output oe file I got:

1684608 reads; of these:
  1684608 (100.00%) were unpaired; of these:
    1680288 (99.74%) aligned 0 times
    48 (0.00%) aligned exactly 1 time
    4272 (0.25%) aligned >1 times
0.26% overall alignment rate

Does this make sense? I doubted it but I am not sure where the problem is. Any suggestions?

[Devon Ryan]

Without checking, my guess is that your awk command isn't sufficient to print the entire sequence. rRNAs can be quite large, so you'll miss a lot of sequence unless it's all on one line.

[niu2rseq]

Hi Devon,

Here the ncrna.fa is quite small, only 998K. It contains only 7228 lines and 3614 sequences. There are only 504 rRNA sequences there.

[Devon Ryan]

Fair enough then. That does seem to be too low an alignment number, though. Are these raw reads or is this after collapsing things?