SAMTOOLS mpileup only INDELs, no SNPs?

[jkozubek]

Hi All,

I used BWA to create BAM files. Then I used samtools mpileup to create a VCF file. However, there were only indels (no SNPs) in my VCF file.

Right now, I am using samtools fillmd to add MD tags to my BAM files.

/home/directory/bin/samtools fillmd -Aeb /home/directory/bam/${file1}.bam /home/directory/Mouse/mm10.fa > /home/directory/bam/${file1}.recode.bam

With a cursory look at my output, I now see the SNPs in the sequence reads (amid lots of =======). However, fillmd seems like it is taking a Very long time. Does anyone know whether the MD tags are essential to SNP calling? In other words, should this fix the issue? Other ideas?

By the way, my vcf line was:
/home/directory/bin/samtools mpileup -B -m 5 -q 5 -S -D -u -f /home/directory/Mouse/mm10.fa -b /home/directory/bam/bam.list | /home/directory/bin/bcftools view -bvcg - > /home/directory/bam/mm10.bcf
/home/directory/bin/bcftools view /home/directory/bam/mm10.bcf | /home/directory/bin/vcfutils.pl varFilter > /home/directory/bam/mm10.vcf


Jim

[jkozubek]

I believe I did not specify -I parameter in bwa to specify Illumina 1.3+ ASCII

[jkozubek]

This was not the answer.
In fact, upon further investigation, I discovered the sequence data we received has very low base quality. I set samtools mpileup to -Q 0 and now get lots of SNPs.

[jkozubek]

Sequencer is coding in Sanger/Illumina 1.9 (ASCii + 33) scale. I noticed BWA recoded the values by subtracting 33. So it adjusted all the base calls. I checked and the base calls are -33 in the BWA output. However, samtools still considers the code (ASCii + 33) by default, so it too subtracts 33. Therefore, the quality code was readjusted twice.