Making correct RNA references

Often, I need to have a nice RNA reference for alignments or similar tasks. What sometimes happens is that the reference sequence you download from a database (for example hairpin.fa from miRBase, or tRNAs from GtRNAdb) contains the real RNA sequence – contains uracils and not thymines. By definition, this is a correct RNA sequence. But, if you want to do an alignment the same way you are used to do alignments to a genome reference you might get some unexpected results. They are often caused by Us in the sequnce. You will get some alignment from sequences that have fewer Us than number of mismatches but definitely the alignment results will be biased. My rule no. 1 when I start to work with any RNA-based reference is to make sure they have Ts and not Us and I also like to have them in one-line fasta format.

For example for miRBase miRNA hairpin reference:

$ wget ftp://mirbase.org/pub/mirbase/21/hairpin.fa.gz # Get the sequence
$ zcat hairpin.fa.gz | perl -pe '/^>/ ? print "\n" : chomp' | tail -n +2 | sed '/^[^>]/ y/uU/tT/' > hairpin.oneline.fasta # Convert to one-line fasta and replace all Us for Ts
$ grep -A1 "^>hsa" hairpin.oneline.fasta | sed '/^--$/d' > hsa_hairpin.fasta # Get only human hairpin sequences

Leave a comment