Question: Differences in read numbers - Bowtie2 output
0
gravatar for geek_y
4.7 years ago by
geek_y9.4k
Barcelona/CRG/London/Imperial
geek_y9.4k wrote:

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 tim

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

 

samtools alignment bowtie2 • 2.2k views
ADD COMMENTlink modified 4.7 years ago by Devon Ryan89k • written 4.7 years ago by geek_y9.4k

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 thanAS:i.

 

 

ADD REPLYlink modified 4.7 years ago • written 4.7 years ago by Ashutosh Pandey11k

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

ADD REPLYlink modified 4.7 years ago • written 4.7 years ago by geek_y9.4k

That's nice, but completely irrelevant.

ADD REPLYlink written 4.7 years ago by Devon Ryan89k

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. 

ADD REPLYlink written 4.7 years ago by geek_y9.4k
0
gravatar for Devon Ryan
4.7 years ago by
Devon Ryan89k
Freiburg, Germany
Devon Ryan89k wrote:

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.

ADD COMMENTlink modified 4.7 years ago • written 4.7 years ago by Devon Ryan89k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1498 users visited in the last hour