I am having an issue with samtools that I wondered if you might know. I have been using samtools to call snps using a .bam file which I also created using samtools. I was only interested in snps and not indels, so I used the -I option in mpileup, in addition I was hoping to filter the snps using "vcfutils.pl varfilter -w 1000" in order to filter out an snps within 1kb of a gap.
The issue, however, is that I only seem to be able to use one option or the other. That is, if I use -I to eliminate indels, then I can not use "-w 1000" to filter them, and I end up with the same number of snps (in fact the same file) no matter what value is used for -w. However, when I do not use the -I option, then I end up with a much smaller number of snps when using "-w 1000."
I have provided the commands below:
"-I -w1000" --> 26,395 snps ( results in same file as immediately below)
samtools mpileup -ugf referenceGenome.fasta sortedIndexed.bam -r Gm02:1000000-15000000 -E -S -C 60 -D -d 10000 -I | bcftools view -Nbcvg - > tmp.bcf
bcftools view tmp.bcf | vcfutils.pl varFilter -d 8 -D 100 -w 1000 - > result.vcf
"-I -w100" --> 26,395 snps ( results in same file as immediately above)
samtools mpileup -ugf referenceGenome.fasta sortedIndexed.bam -r Gm02:1000000-15000000 -E -S -C 60 -D -d 10000 -I | bcftools view -Nbcvg - > tmp.bcf
bcftools view tmp.bcf | vcfutils.pl varFilter -d 8 -D 100 -w 100 - > result.vcf
"-w1000" --> 7,444 snps
samtools mpileup -ugf referenceGenome.fasta sortedIndexed.bam -r Gm02:1000000-15000000 -E -S -C 60 -D -d 10000 | bcftools view -Nbcvg - > tmp.bcf
bcftools view tmp.bcf | vcfutils.pl varFilter -d 8 -D 100 -w 1000 - > result.vcf
"-w100" --> 22,001 snps
samtools mpileup -ugf referenceGenome.fasta sortedIndexed.bam -r Gm02:1000000-15000000 -E -S -C 60 -D -d 10000 | bcftools view -Nbcvg - > tmp.bcf
bcftools view tmp.bcf | vcfutils.pl varFilter -d 8 -D 100 -w 100 - > result.vcf
Does anyone know if mpileup -I an varFilter -w are incompatible?
