Count multi-maped reads independent of mapper

Question

[Tomislav Ilicic]

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

[David Langenberger]

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)?

[Tomislav Ilicic]

I am trying to do both.

I want to avoid outputting multi-mapped reads and subsequently confirm that there really are none.