Count intervals of multiply mapped reads overlapping the genome

[sunkissbean]

Hi everyone,

I have some SAM/BAM files containing the alignments of small RNA-seq reads to mm9 that are created by BWA. I'm interested in calculating where they are mapped to in the genome. I noticed that there are a lot of reads that are mapped to multiple loci (multi-mappers) in genome. Therefore, I first separated the unique-mappers from multi-mappers and counted intervals of unique-mappers overlapping different regions of the genome using the bedtool. Now here comes the issue: as multi-mappers are mapped non-uniquely to various regions in the genome, if I simply use the same method as the uniquely-mappers, I will overestimate the number of counts, how can I count the normalized intervals of multi-mappers overlapping different regions of the genome? Is there any tool in the galaxy or R that I can use? Thank you for the help in advance!

[dpryan]

The thing with small RNAs is that many of them have multiple copies in the genome. For that reason, it's often easier to simply align to the transcriptome and then either remove multimappers or use RSEM/eXpress/etc. to get estimated per-feature counts.

[sunkissbean]

That's a very good suggestion! I tried using RSEM to quantitate the reads from SAM/BAM files that were generated by BWA. Somehow, I always got an error message reading "tool error" (btw, I used RSEM from the tool shed in Galaxy for this purpose). I feel that the error came from when I prepared the reference as the reference file that I generated was very small (<1 MB). Specifically, I used the RSEM prepare reference function to generate the reference from RefSeq DNA fasta file. Then, I used the RSEM calculate expression to quantitate the expression level of small RNAs by using a SAM/BAM file that was generated by BWA and the reference file. I am really confused what went wrong now.