0
$\begingroup$

I have a dataset of paired-end reads, which can overlap. So when i run bcftools mpileup this overlap gives a double count in the depths so we see DP = 2 even though it should be DP = 1.

The question is similar to

Double-counting coverage of overlapped read pairs

However this is more in relation to bcftools mpileup. If I use the same test data as in the link above

cat samtest.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

samtools mpileup samtest.sam | head -52 | tail -9
1   44  N   1   T   <
1   45  N   1   A   <
1   46  N   1   T   6
1   47  N   1   A   6
1   48  N   1   A   6
1   49  N   1   T   W
1   50  N   1   T$  6
1   51  N   1   a   <
1   52  N   1   a   <

samtools mpileup -x samtest.sam | head -52 | tail -9
1   44  N   1   T   <
1   45  N   1   A   <
1   46  N   2   T^Sa    <<
1   47  N   2   At  <<
1   48  N   2   At  <<
1   49  N   2   Tt  <<
1   50  N   2   T$a <<
1   51  N   1   a   <
1   52  N   1   a   <

So here we see that samtools handles the double count issue, and if we use the -x / --ignore-overlaps we see the double count as we would expect.

The issue i am having is how do i get bcftools mpileup to NOT count this overlapping region (pos 46 - 50) twice. If i am just runnning bcftools mpileup in the simplest form similar to the first samtools mpileup command i get the following:

bcftools mpileup -f ref2.fa samtest.bam | head -74 | tail -9

1   44  .   T   <*> 0   .   DP=1;I16=1,0,0,0,27,729,0,0,50,2500,0,0,6,36,0,0;QS=1,0;MQ0F=0  PL  0,3,27
1   45  .   A   <*> 0   .   DP=1;I16=1,0,0,0,27,729,0,0,50,2500,0,0,5,25,0,0;QS=1,0;MQ0F=0  PL  0,3,27
1   46  .   T   <*> 0   .   DP=2;I16=1,0,0,0,21,441,0,0,50,2500,0,0,4,16,0,0;QS=1,0;MQ0F=0  PL  0,3,21
1   47  .   A   <*> 0   .   DP=2;I16=1,0,0,0,21,441,0,0,50,2500,0,0,3,9,0,0;QS=1,0;MQ0F=0   PL  0,3,21
1   48  .   A   <*> 0   .   DP=2;I16=1,0,0,0,21,441,0,0,50,2500,0,0,2,4,0,0;QS=1,0;MQ0F=0   PL  0,3,21
1   49  .   T   <*> 0   .   DP=2;I16=1,0,0,0,27,729,0,0,50,2500,0,0,1,1,0,0;QS=1,0;MQ0F=0   PL  0,3,27
1   50  .   T   <*> 0   .   DP=2;I16=1,0,0,0,21,441,0,0,50,2500,0,0,0,0,0,0;QS=1,0;MQ0F=0   PL  0,3,21
1   51  .   A   <*> 0   .   DP=1;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0  PL  0,0,0
1   52  .   T   <*> 0   .   DP=1;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;MQ0F=0  PL  0,0,0

NOTE: the ref2.fa is simply the sequence of the read repeated 3 times, just to create a basic test sample set.

So here I again have the DP=2 and i get the same result if i use the -x or -Q 1 as suggested for samtools mpileup in the question posted earlier in the link above.

Do anyone know how to get bcftools mpileup to just count DP=1 in pos 46 - 50 similar to the samtools pileup results?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.