samtools mpileup giving QUAL of 0
1
0
Entering edit mode
6.6 years ago
blur ▴ 220

Hello all,

An innocuous mpileup gives QUAL of 0 for all bases. 0 shouldn't even be a legal QUAL value (as it is the result of an exponent, according to VCF format documentation). This is the command we're running:

samtools mpileup --VCF -r ChrA_C_glabrata_CBS138:1-100,000 --output-tags DP,DP4 -u -f CBS138.fasta 12.sorted.bam > 12.mpileup.txt

A sample output line:

ChrA_C_glabrata_CBS138    337    .    G    A,C,T    0    .    DP=31;I16=4,11,7,8,530,19116,501,17779,6,6,2,2,298,6778,286,6564;QS=0.5,0.4,0.0666667,0.0333333;VDB=0.523807;SGB=-0.688148;RPB=0.768403;MQB=0.461075;MQSB=0.999074;BQB=0.952765;MQ0F=0.741935    PL:DP:DP4    0,7,3,30,26,5,36,31,9,5:30:4,11,7,8

Thanks!

 

samtools mpileup vcf • 1.9k views
ADD COMMENT
0
Entering edit mode

Another example, that might explain the diploid assignment:

Total count: 95
A      : 0
C      : 67  (71%,     33+,   34- )
G      : 0
T      : 28  (29%,     12+,   16- )
N      : 0
---------------

ADD REPLY
0
Entering edit mode

As I said previously, samtools assumes you have a diploid organism. Sites like the one above suggest that either you don't have a haploid organism, you have multimappers, you have incorrect alignments, or you have notable genetic heterogeneity.

ADD REPLY
0
Entering edit mode
6.6 years ago

0 is certainly legal, since 10^0 = 1, meaning that there's a 100% chance that the call is wrong. This is unsurprising given that all four possible base changes are apparently present at similar depths and samtools assumes a diploid genome.

ADD COMMENT
0
Entering edit mode

my genome is haploid, and all the SNPs I got look like the line above...

Also, my DP doesn't match the read depth - I opened these files in IGV beforehand, and for example, if IGV has 94 reads, the DP is marked for the same position as 31....

mpileup problem?

ADD REPLY
0
Entering edit mode

There's some filtering that happens before the DP metric is made. This is unlikely to be an mpileup problem per se but rather either a usage problem or simply due to having a bunch of multimappers (otherwise, the depths of the alternate sequences makes no sense). Are you trying to use this for some sort of pooled input?

ADD REPLY
0
Entering edit mode

no... this is what worries me...

ADD REPLY
0
Entering edit mode

What do things look like in IGV in this area?

ADD REPLY

Login before adding your answer.

Traffic: 2818 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6