How to get multi mappers from bowtie 1

Question

[Amaranta_Remedios]

So I know there is a lot of info out there but I can't really work out who is "right". I am running bowtie 1 using small RNA-Seq data and I don't have a problem getting only the unique mappers, bowtie can do that by adding the -m 1 option. However I still can't work out how to extract the multi mappers, when using another option like -best, or -m 5 for instance. Any suggestion will be appreciated.

Thanks!

Answer 1

[Carlo Yague]

Multimappers have a mapping quality of 0, so you can pipe the output of bowtie into a samtools view command that will remove all reads with mapping quality above 0.

samtools view -b -U -q 1 mapped.bam > multimappers.bam

(-q 1 selects for reads with mapping quality >=1 ; -U inverts that selection)

[Amaranta_Remedios]

I will try that, thanks for the answer !

[xiaoleiusc]

Hi Carlo, Thanks for your answer, but all my short RNA reads from bowtie 1 does not have mapping quality score (all of them have 255). In addition, there is no NH value in the output sam file.

[Carlo Yague]

If bowtie1 didn't assign mapping quality, then the easiest/best way might be to remap the reads using a more modern aligner, then filter as suggested above. If for some reason you do not want to do that, another solution would be to extract alignments whose read name occur more than once in the .bam file (because if it occurs more than once, then it is a multimapper). Here is how I would do it (code untested).

# First, lets extract read names that are found at least twice in the bam file
samtools view mapped.bam | cut -f 1 | sort | uniq -d > read_name_multi.txt

# Then filter the bam file based on those read names using picardtools (you could also use grep, but it is slower)
java -jar FilterSamReads.jar INPUT= mapped.bam FILTER=includeReadList READ_LIST_FILE=read_name_multi.txt OUTPUT=multimappers.bam