Calculating total no. of multimapped reads using samtools for RNA-Seq data analysis

[nareshmvr]

Dear all ,

please excuse me if my questionaire is wrong.

I am confused in estimating proper command for calculating number of multimapped reads using samtools for my paired end RNA-Seqdata, when i tried searching forums , i found two
threads for calculating total no. of multimapped reads

http://seqanswers.com/forums/showthread.php?t=16427

http://seqanswers.com/forums/showthread.php?t=11967

I tried using both , the results were very different , please suggest me proper command to calculate multimapped reads using samtools

[BnaT]

I'd try filtering with FLAG 0x100 or 256. Sometimes I have seen inconsistencies depending on the mapper being used (e.g GEM mapper).

[nareshmvr]

I have used Tophat for mapping my genomes .

[blakeoft]

If I recall correctly, I wasn't able to actually find this in the TopHat documentation, but take a look at this discussion. According to the top answer, the mapping quality encodes the number of mappings.

255 = unique mapping
3 = maps to 2 locations in the target
1 = maps to 3-4 locations
0 = maps to 5 or more locations (up to the number defined in "--prefilter-multihits")

So I guess you could do something like

samtools view file.bam | awk '$5 < 255 {count++} END {print count}'

Edit: this might be inefficient, I'm really not sure. But I think it'll get the job done.