samtools mpileup parameter settings

[xxy39]

Hi all,

I am using samtool mpileup to do SNP calling. I am experiencing some difficulties with the parameter settings.

For the single end data with a 200X coverage, samtools gave me the similar number of SNPs as SOAPsnp, GATK-UGT and Atlas-SNP2. However, for a pair-end data with 1000X coverage, samtools gave much less SNPs than the other callers.

For the pair end data, I tried the following code:
# samtools mpileup -ug -C50 -d1000 -f Feb11.2011.combined.25chr.fa 826T.fixed.cleaned.bam | bcftools view -bvcg - > 826T.samtool.mpileup.bcf
# bcftools view 826T.samtool.mpileup.bcf | vcfutils.pl varFilter -D1000 > 826T.samtool.mpileup.c50d1000.vcf

I tried the parameters "-C", "-d", and "-D" with several different numbers, but they all gave me the exactly same results. Based on the samtools instruction, the -d and -D are associated with coverage, and changing them around should make some difference. And I tried with older version (v.1.12) and new version (v.1.16) samtools. These parameter settings seem not working at all. Do anyone know what might go wrong? What should I do to get more SNPs called by samtools?

Any help appreciated!

[DZhang]

Hi,

You may want to examine the bam file using any genome browser. IGV from Broad is a good option, for example. Check the location(s) that other programs called a SNP but samtools did not. You may gain some insights from this practice. All these programs may behave strangely in different situations...

Douglas

https://www.contigexpress.com