not an expert in bowtie2, rather BWA, but, by definifion the multiple mapping read is the one, which has more than 1 location in the genome, where the alignment scores are identical. Easiest way to identify such reads is getting the ones that have the mapping quality field in SAM file set to 0. That is getting the multiple mapping reads according to their definition

Thanks for sharing your opinion. MAPQ score of 255 engenders uniquely aligned reads; not to mention with the highest mapping score. Based on some literature, and by logic too, isn't it obvious that any number between 0 and 255 will represent multiple aligned reads. So, is it a subjective matter to chose quality score (depiction as per the 5th field in the SAM file) to ascertain multiple reads.

I believe the following will simply do the trick:

> samtools view -b -q your_alignment_file.bam > filtered.bam

Additionally, could you also comment on the discrepancy (if any) with the alignment summary/ statistics from Bowtie2.

MAPQ score can be different depending on the aligener you use and it can be interpreted in different ways. For example, I've been using GMAP/GSNAP as the sole aligner for all my NGS analysis and the MAPQ scores it gives are

40: Maps uniquely to one location on the genome

3: Maps to two locations on the genome

2: Maps to three locations

1: maps to 5 - 10 locations

0: maps to more than 10 locations

There are scores that are in between and those mainly happen if there mismatches in the read e.g. MAPQ 35, 30, 25 ..etc

Filtering reads reads on the basis of MAPQ score depends on the experiment. If it is an RNA-seq experiment then in most cases you should not be very stringent with the filtering. There two ways to filter out low MAPQ reads

1. using samtools

$ samtools view -b -q 30 input.bam > output.bam

2. Using bedtools, only if you want to get transcript read counts in BED format

$ bedtools multicov -bams input.bam -q 30 -bed genome file in gff,gtf, or bed format > output.bed

best of luck

Thanks Waleed. So broadly, higher MAPQ score ensures precise, consolidated unique alignment.

Exactly but filtering low MAPQ reads depends on the experiment:

1. RNA-seq: Should not be very stringent with filtering because you already allow a number (5 - 10) of mismatches to map exon-exon junctions.

2. For experiment like ChIP-seq, RIP-seq, CLIP-seq ...etc where identifying the precisie binding sites of the proteins is absolutely critical then a stringent filtering is required here. With this type of experiments I filter out reads with MAPQ < 25 (this is based on GSNAP MAPQ score).

I used Bowtie2 only a few times and I believe it has the same MAPQ scoring as TopHat, for a reference here is the interpretation of the bowtie MAPQ scores:

0 = maps to 5 or more locations

1 = maps to 3-4 locations

3 = maps to 2 locations

255 = unique mapping

Dear,

i used BWA mem alignment, now i have reads with MAPQ 0 that i would like to extract into a new sam or bam file. Do you know if it is possible to do that with Samtools or another tool?

Untill now i tried:

samtools view -bSq 1 input.sam > output.bam

which filter out all reads with mapping quality lower then 1, but i would like to extract them to a new sam/bam file and work on them.

Thank you for your help!