Question: Bowtie2 Reads that aligned concordantly exactly 1 time
1
gravatar for AW
4.4 years ago by
AW350
United Kingdom
AW350 wrote:

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!

 

 

 

 

 

samtools bowtie2 • 3.4k views
ADD COMMENTlink modified 10 months ago by Biostar ♦♦ 20 • written 4.4 years ago by AW350

Thanks, thats 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

 

 

 

ADD REPLYlink written 4.4 years ago by AW350
0
gravatar for geek_y
4.4 years ago by
geek_y9.9k
Barcelona
geek_y9.9k wrote:

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

ADD COMMENTlink written 4.4 years ago by geek_y9.9k
0
gravatar for Devon Ryan
4.4 years ago by
Devon Ryan92k
Freiburg, Germany
Devon Ryan92k wrote:

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.

ADD COMMENTlink written 4.4 years ago by Devon Ryan92k
0
gravatar for AW
4.4 years ago by
AW350
United Kingdom
AW350 wrote:

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!

 

 

ADD COMMENTlink written 4.4 years ago by AW350

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

ADD REPLYlink written 4.4 years ago by Devon Ryan92k
0
gravatar for AW
4.4 years ago by
AW350
United Kingdom
AW350 wrote:

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!

ADD COMMENTlink modified 4.4 years ago • written 4.4 years ago by AW350

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.

ADD REPLYlink written 4.4 years ago by Devon Ryan92k
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: 531 users visited in the last hour