This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Bowtie2 Reads that aligned concordantly exactly 1 time

Hi,

I have mapped some Illumina reads onto my reference genome using Bowtie2. I only specified --nounal and left all other parameters at default.

111247311 reads; of these:
  111247311 (100.00%) were paired; of these:
    19249548 (17.30%) aligned concordantly 0 times
    81185571 (72.98%) aligned concordantly exactly 1 time
    10812192 (9.72%) aligned concordantly >1 times

I wish to extract the 81185571 reads that aligned concordantly 1 time.

I can extract concordant reads with samtools view -f 0x2, but I wonder how Bowtie2 defines reads that map concordantly exactly 1 time?

Can I extract these using the absence of XS:i or is it based on a MAPQ score threshold?

Thanks!

bowtie2 samtools
• 7,767 views • link
updated by Ram Ram • written by AW AW

2 answers

How To Extract Reads-Pairs Aligned Concordantly Exactly 1 Time?

• 0 views • link
updated by Ram Ram • written by GouthamAtla GouthamAtla

Thanks, that's a very helpful link!

My question was slightly different though. I want to know how Bowtie2 defines reads that aligned concordantly exactly 1 time?

Does it use the absence of XS:i or is it based on a MAPQ score threshold?

Thanks

• 0 views • link
updated by Ram Ram • written by AW AW

I vaguely recall that it looks for alignments where AS:i == XS:i, which will have a MAPQ of 0 or 1. Note that an alignment can also have a MAPQ of 1 if this is not the case, so you can't just samtools view -cq 2 blah.bam and get the same number.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

Thanks for your answer.

I have tried the following to get the 81185571 read pairs that aligned concordantly 1 time.

1. Extract concordant pairs where both reads map only once (i.e. unireads). This is based on absence of XS:i: flag.

Number of pairs remaining (absence of XS) = 74507750

2. Extract concordant pairs where AS:i == XS:i for both reads.

Number of pairs remaining (AS != XS) = 88621585

3. Extract concordant pairs where AS:i > XS:i for both reads.

Number of pairs remaining (AS > XS) = 83088150

Still cannot get the same number reported?

Any other suggestions would be greatly appreciated!

• 0 views • link
updated by Ram Ram • written by AW AW

You'll have to go through the source code then. There's no documentation on this.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

Yes, I found that the aln_sink.cpp file contains the instructions for printing number of concordant reads etc and that repetitive mapping appears to be specified by the pairMax value, but I have no idea where or how the pairMax value is set.

Thanks for your help!

• 0 views • link
updated by Ram Ram • written by AW AW

Welcome to the joys of figuring out undocumented parts of packages. At some point you have to ask yourself if getting the answer is really worth the hassle.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this