search

samtools / bcftools

This is the official development repository for BCFtools. See installation instructions and other documentation here http://samtools.github.io/bcftools/howtos/install.html
http://samtools.github.io/bcftools/
Other
633 stars 241 forks source link

bcftools mpileup + call seemingly confident in allele only present in a few reads? #2185

Open emos8710 opened 1 month ago

emos8710 commented 1 month ago

Hi!

I have an issue that is puzzling me. I have two salmonella samples that I'm comparing and expect them to be highly similar if not identical. I have mapped them with bowtie2 (three different versions as I was investigating if it was the mapping that was the issue) and created a mpileup and called variants with bcftools v 1.14 and 1.17 (no difference in the result)

bcftools call -m -v -O v --ploidy 1

I'm getting some variant calls that I don't understand. Here is the best example. There are three copies of each sample as I was checking if there was any difference between bowtie2 versions.

NC_003197.2 1349066 .   C   T   720.935 .   DP=815;VDB=0.00653457;SGB=-21.0274;RPBZ=-12.8027;MQBZ=-8.78292;MQSBZ=4.75308;BQBZ=6.99271;SCBZ=0;FS=0;MQ0F=0.00368098;AC=3;AN=6;DP4=4,76,70,101;MQ=37   GT:PL   0:0,112 0:0,112 0:0,118 1:255,0 1:255,0 1:255,0

I have also tried running the mpileup with a MQ threshold of 20, which barely makes a difference. For sample 1 there are 23 reads supporting REF and 126 reads supporting ALT. For sample 2 there are 8 reads supporting REF and 112 reads supporting ALT (plus 2 reads supporting a second ALT). However bcftools seems pretty confident in the REF call for sample 1 and very confident in the ALT call for sample 2. Why is this? I would expect more uncertainty, so I could filter this position :)

Lilu-guo commented 1 month ago

Yes, I also found that the specified MQ threshold didn't seem to work. For example, when using -q 20, it still uses 30 as the threshold; and when using -q 40 (greater than 30), it uses 40 normally.

Maybe we can discuss bcftools further. It will be more efficient to use WeChat or other instant tools.

pd3 commented 1 month ago

I am unable to comment without seeing the data, ideally a small bam slice + corresponding fasta reference chunk.

EDIT: This script can be used to create a small test case https://github.com/pd3/mpileup-tests/blob/main/misc/create-bam-test

emos8710 commented 1 month ago

Thanks for your reply, I will get back with some data!