Differences in read numbers - Bowtie2 output

Question

[GouthamAtla]

Hi All,

I have done bowtie2 alignment with the following command:

bowtie2 -x Index_BaseName trim_11CS.fq -S 11CS.sam -p 30

The output I got is:

2679302 reads; of these:
  2679302 (100.00%) were unpaired; of these:
    1760057 (65.69%) aligned 0 times
    846305 (31.59%) aligned exactly 1 time
    72940 (2.72%) aligned >1 times
34.31% overall alignment rate

But when I see the stats with samtools:

samtools view -S -q 42 11CS.sam | wc -l

or

samtools view -c -S -q 42 11CS.sam

The output is 69678

I am trying to extract 846305 (31.59%) aligned exactly 1 time

I am really confused with the numbers reported by bowtie2 and samtools. I am trying to find out uniquely mapped reads with highest mapping quality assigned by bowtie2 i.e 42

[Ashutosh Pandey]

Not all the reads that have been aligned uniquely will be assigned with the highest MAPQ. You may find reads with MAPQ of 30 and still part of the group of reads that aligned exactly 1 time. It is advisable to filter your bam file using MAPQ scores rather than selecting reads that are uniquely aligned. In case you still want to retrieve the uniquely aligned reads you can use XS:i: filed in the SAM output of bowtie. Can you see XS:i flag in your sam output?

Here is the description of that flag from Bowtie2 manual.

Alignment score for the best-scoring alignment found other than the alignment reported. Can be negative. Can be greater than 0 in --local mode (but not in --end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than AS:i.

[GouthamAtla]

I have XS:i:somedigit in my sam file that I got from bwa-mem alignment

[Devon Ryan]

That's nice, but completely irrelevant.

[GouthamAtla]

I mean I do not have XS tag in bowtie2 output but it's there in BWA alignment. Wih BWA alignment I am using mapQ 30 as my threshold as mentioned Devon.

Answer 1

[Devon Ryan]

There's no perfect correspondence between MAPQ and what bowtie2's output metrics consider as unique (see the reply from Ashutosh Pandey). In general, you should completely forget the incorrect concept of "uniquely aligned", since there's actually no such thing. Just pick a reasonable MAPQ threshold (10 is good for many purposes, 5 for others, etc.) and use those.