I have been using using the samtools pipeline for variant calling for a while now and have just moved over to using mpileup. I am using the suggested method:
- samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf
- bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf
I am aware at the vcfutils step that '-d' specifies the minimum read depth and '-D' specifies the maximum read depth. However, what is not clear to me is, for example, if -D = 100 and the coverage of a SNP is 200, what happens? Do 100 random reads get used, or more worryingly is that SNP not reported?