[Samtools-help] samtools mpileup genotype likelihood calulation

[Günter J. <Gue...@ce...>]

Dear Samtools-Team,

I encountered that samtools mpileup followed by bcftools for variant 
calling produces different results when the bam file is sorted 
differently. Since the bam file has to be sorted by coordinate, the 
order of reads sharing the same start position may vary in different 
versions of the bam file, i.e. if reads in a coordinate sorted bam file 
are shuffled and then sorted by coordinate again, the order of reads 
mapping to identical start positions usually differs from the order of 
the same reads in the initial bam file.

This however, does change the quality of the variants that are called by 
bcftools. Is suspect that the genotype likelihoods change. I would have 
expected that changing the order of reads mapping to the same start 
position does not influence the variant calling, when using the exactly 
same samtools mpileup/bcftools commands. I use the current version 1.3., 
but also older versions show that problem.

Is there any explanation why this is the case? The problem occurs with 
every bam file I tested so far. It seems like samtools mpileup does not 
use all of the reads overlapping a variant position for likelihood 
estimation?! Is there any option in bcftools or samtools mpileup to 
prevent this from happening, i.e. to produce the same variant quality 
scores/genotype likelihood values regardless of the order of reads?

I would appreciate any comments on that.

Best regards,
Günter