Hi All,
I have been searching a lot and haven't really found a good conclusion to my problem.
All I am trying to do, is to count the number of multi-mapped reads across different mapping algorithms (Bowtie1, Bowtie2, BWA MEM, STAR, Tophat, GMAP).
I know that HTSeq count required the flag NH:I: to indicate if an alignment is unique or not. However, only Star seems to output that correctly.
Now I have found bamstat.py
from RSeQC (http://rseqc.sourceforge.net/) and it does it automatically. However, they only look on the alignment quality and I wonder whether this is the best way to do it (the correct way).
Also, I have been trying to forbid every mapper to output multi-mapped reads and only STAR actually seems to be doing this according to bamstat.py.
These are the options I chose for each mapper to forbid multi-mapped reads in the SAM (x I set to 1).
BWA:
mem: outputs randomly one of multi-mapped alignments or all (-a) - no way to put limitations
Bowtie1 and 2 :
-k prevents to output more than x alignments
GMAP:
-n X -Q (don't print alignment if more than X are found)
Star:
outFilterMultimapNmax x(print alignments if not more than X are found)
Can you give me some help on how to do it correctly?
Best,
Tomi
What exactly is your problem? You want to count multi-mapped reads correctly (as your header says), or avoid multi-mapping reads (as your text says)?
I am trying to do both.
I want to avoid outputting multi-mapped reads and subsequently confirm that there really are none.