This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BCFtools not calling variant, despite adequate coverage

I have four BAM files that I know to have a G->A mutation at chr20:63,439,704

An example in IGV is below: enter image description here

However, when I run bcftools mpileup and call:

bcftools mpileup -Ou -f hg38.fa Q201_5bamheader.sam | bcftools call -m > Q201_5_allsites.vcf

it doesn't call the variant. The output of bcftools mpileup is:

chr20 63439704 . G A,<*> 0 . DP=24;I16=4,9,4,6,504,20028,370,13690,0,0,0,0,270,5980,250,6250;QS=0.565217,0.434783,0;VDB=0.144519;SGB=-0.670168;RPBZ=-1.96622;MQBZ=0;MQSBZ=0;BQBZ=-0.877058;SCBZ=0;MQ0F=0.958333 PL 0,6,5,39,35,7

And for bcftools call, the output is:

chr20 63439704 . G . 29.5119 . DP=24;VDB=0.144519;SGB=-0.670168;RPBZ=-1.96622;MQBZ=0;MQSBZ=0;BQBZ=-0.877058;SCBZ=0;MQ0F=0.958333;AN=2;DP4=4,9,4,6;MQ=0 GT 0/0

I already tried setting the prior to 1, but it still didn't call it. Can anyone help me troubleshoot?

bcftools variant calling
• 1,916 views • link
updated by GenoMax GenoMax • written by mpiersonsmela • mpiersonsmela

1 answer

unless I'm wrong all the reads in this IGV screenshot look "blank" which is usually the color for the reads of low quality, duplicate, etc...

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

what is the output of samtools mpileup for this position ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

The output is:

chr20   63439704    G   23  ,,,,,aa,.A.A,,a,a.aA.aA FFFFFFFFkFFFFFFFFFFFFFF
• 0 views • link
updated by GenoMax GenoMax • written by mpiersonsmela • mpiersonsmela

The reads are not duplicates. I think the color might be due to also matching the chr20 alt contig but I'm not sure.

An example line from the sam file is:

A00808:347:H5GKLDSXY:4:1322:13774:30859 163 chr20   63439615    0   150M    =   63439796    331 AGGAGACACCGATGAGGGTGAAGGTTGCCGCAAGGAGCCTGCCGTTCCAGGTCTGGGGGTACTTGTCCCCGTAGCCAATGGTGGTCAGCGTGATCTGTGGGACCGCAGGCTCTAGTCACACGAAGGGCCTGCTCACACCCCTGAGGGCAG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFFFFFFFF:FFFFF:F:FFF:FFFFFFFF:FFFFFFFFFFFFFFFFFFFF,FFFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFFFF  NM:i:0  MD:Z:150    MC:Z:150M   MQ:i:0  AS:i:150    XS:i:150    XA:Z:chr20_KI270871v1_alt,+24898,150M,0;
• 0 views • link
updated by GenoMax GenoMax • written by mpiersonsmela • mpiersonsmela

your read is discarded because mapping quality=0

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Strange, I wonder why. Could it have anything to do with the "chr20_KI270871v1_alt"?

• 0 views • link
written by mpiersonsmela • mpiersonsmela

yes, this seems to be the problem. I will try realigning to a version of hg38 without alt contigs.

• 0 views • link
written by mpiersonsmela • mpiersonsmela

Log in to answer this question.

More posts like this