Samtools multisample SNP caller

[amcrisan]

Hi -

I am running a comparison between running samtools SNP calling in the multi sample mode, vs using it on each file individually. I seem to be different getting per-sample genotypes when using one method compared to the other, I wanted to know if anyone had experienced similar difficulties before.

Both methods of calling SNPs use same command, except in one case I use the -b command and provide paths to the bam files I want, and in the other case I pass each file individually and aggregate them later. Just to be on the same page this is the command:

samtools mpileup -uf ref.fa aln1.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf

where in some cases, instead of aln1.bam I provide -b list, which contains a list of bam files.

When using the multi-sample way (-b command), the vcf file also contains the genotype information for each individual sample. The output format it comes in is GT:PLP:GQ, and the issue that I have is the GT seems to be wrong when I use the multi-sample call. For example, one of my samples is said to have a 1/1 (ALT,ALT) genotype at a specific position, BUT when I review in IGV this is clearly incorrect, furthermore when I run the sample on it's own (not multi-sample mode) this positions DOES NOT show up as containing a SNP. So there is something going on when I pass a list, that seems to mess up the GT (that's what I've concluded).

I think that it's a really weird error. I am trying ti diagnose what the problem could be and am wondering if anyone else has noticed it before.

[amcrisan]

After some looking into this further and getting some suggestions and help from others, this is what I think the difficulty I had was.

The reason why the multi-sample vs. single sample variant calling is different is because the multi-sample version uses information across all samples to call a genotype. There is an obvious advantage to this approach, which is if you have low coverage but some evidence for the genotype that is present in all the other samples then you are more likely to believe that the variant is real in low coverage data. This means that by aggregating information across samples, variant detection is more sensitive compared to a single sample approach.

I had suspected it was doing this but still couldn't get why it would call a variant when NO evidence of the non-reference allele was present. Then a piece of advice pointed me in I think the right direction: I should consider the prior. If you're familiar with Bayeisan methods, the posterior probability of a variant is the prior * the likelihood. The prior would be developed on a per-position basis and in the case of the multi-sample variant calls across all of the samples. If you had overwhelming evidence for BB genotypes at that position, the prior for BB would large and the influence of the likelihood would be negligible. The only time that you really could overcome the influence of the prior would be if at that position you had a very large depth, in which case the likelihood distribution would be so narrow that the likelihood for BB would be 0 or very close to it (if the genotype were actually AA), but in areas where the coverage is low, the likelihood distribution is quite wide, meaning that it is quite unsure of what genotype it could be and the prior would dominate and potentially call a false positive.

[hmk]

Hi,
I'm working on multi-sample variant calling using samtools. And I encountered a problem: what value should I set for the option 'm' in bcftools view.
Because at https://github.com/samtools/samtools/wiki/FAQ, it is said that samtools "Recent versions (tagged as 0.1.19 or newer) implement a new calling model (bcftools view -m) where multi-sample calling comes at no sensitivity cost. It also fixes some of the known issues with multiallelic sites. This is currently the recommended way of calling.". And I can not see a default value for 'm' when I simply type './samtools-0.1.19/bcftools/bcftools view'.
What value should I set for 'm', or just type '-m' with no argument?
Another small question is that if I could do multi-sample variant calling on individuals from different breeds but the same species using samtools?
Thanks!

[hmk]

Luckily, I solved the problem by googling. Here is the answer https://sourceforge.net/mailarchive/...sg_id=30716820, and a good value for 'm' to use is 0.99.