Missing SNPs in samtools output eventhough it display VDB values

[k2bhide]

Samtools Version (0.1.18) error: I called SNPs for 17 samples from same organism, however some of the samples failed to list any SNPs on some of the chromosomes even though I expected some SNPs. Hence I looked in IGV browser and pileup information in order to verify whether there are any alternate alleles at the base position. It was clear that these positions had all mapped reads with alternate bases as compared to reference bases.

I removed –v option from bcftools commands to get depth information for all the base positions irrespective of being called variant sites.

Here is the command I used
samtools mpileup -6 -Buf ../gff_and_Refgenome/ref.fasta in.bam | bcftools view -bcg -> in.var.raw.bcf
bcftools view in.var.raw.bcf > vcfutils.pl varFilter -D50 > in.var.fil.vcf
which gave following result for one of probable SNP positions, failed to be listed as SNP (AtoT) in final result file.

Chr2 617743 . A . 28.2 . DP=95;VDB=0.0535;;AC1=2;FQ=-30 PL 0

I would really appreciate if anyone would help me understand why such SNPs were ignored in final file eventhough it display VDB values which are related to SNP calling in samtools.

[mbayer]

Hi k2bhide,

here are a couple of things to try:

- try the mpileup without the -B option, in case this is filtering out SNPs
- try without the filtering script
- set the minimum base quality for mpileup to 0, in case the variants you are observing happen to be read errors
- set the SNP probability to low for bcftools (e.g. "-p 0.9999")

All of the above will make the command output a lot of dross/false positives but it is useful for the purpose of understanding what's going on.

You then need to tighten things up again to get your actual results. My approach is to make the initial SNP calling very relaxed and then filter appropriately, but if a SNP isn't there after the initial calling then it's obviously gone for good and you can't get it back.

I typically use something like this:

samtools mpileup -Q 0 -D -g -u -f <fastaFile> <bamFile> | bcftools view -vcgA -p 0.9999 - > SNPs.vcf

cheers

Micha