Calling SNPs with samtools: how to specify how many supporting observations are required?

Question

[Lina F]

Hi all,

I am interested in calling SNPs for some alignment data that I have. I am looking at solutions with samtools mpileup and freebayes.

I noticed that freebayes has an option to specify how many observations have to support a variant before it is called (the -C option).

Is there something similar for the samtools mpileup option?

I've been using samtools as follows:

samtools mpileup -uf reference.fasta Sample_sorted.bam -d 8000 | bcftools call -mv -Oz -o variants.vcf

However, I don't know how many observations supported the resulting variants.

An example SNP call from the above command is as follows:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample_sorted.bam
Ecoli_genome    5159    .   A   G   228 .   DP=303;VDB=0.658475;SGB=-0.693147;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;AC=2;AN=2;DP4=0,1,120,171;MQ=60   GT:PL   1/1:255,255,0
Ecoli_genome    5162    .   A   T   228 .   DP=310;VDB=0.633206;SGB=-0.693147;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,120,172;MQ=60 GT:PL   1/1:255,255,0
...

Any insight is appreciated!

Answer 1

[Brian Bushnell]

"DP4=0,1,120,171" tells you the number of reads at that position, in this format:

DP4=refPlus,refMinus,varPlus,varMinus

So, in that case (for the first variant) there was 1 read supporting the reference, and 120+171=291 supporting the variant.

Incidentally, if you're trying out variant callers, the BBMap package also has one callvariants.sh) that is extremely fast and easy to use. It has a flag "minreads" that can be set to specify the number of observations required to call a variant, as well as various other filters.

[Lina F]

Thanks for pointing me towards the DP4 field!