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!



"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.