pileup output different using Maq and Samtools commands

[mr_boourns]

Dear all,

We are trying to detect allelic specific expression with RNA-seq data by mapping against a masked genome (with the reference allele substituted for N's at known SNP positions). We have followed the approach of the Degner/Pritchard paper (2009) where they perform the mapping step using Maq. Our approach has been to use the pileup command in Maq to derive the allele counts at each position. We thought we were getting good results. However, we recently realized that the counts from the Maq produced pileup file do not agree with the counts from a visual inspection of the BAM file produced from the mapping. (This BAM file was produced by converting the .map file from Maq to a .bam file using Samtools). However, when we use the Samtools pileup command on this BAM file, the counts in the pileup file do agree with a visual inspection of the counts when viewing the BAM file in a genome browser.

For instance, for the Maq produced pileup file we get 24 T's and 0 C's at a particular base, whereas in the Samtools produced pileup file we get
14C's and 35T's.

Does anyone have any idea of what is happening here? Should we trust the Maq or the SAM tools results?

Thanks so much for the help,

John

[colindaven]

I would use Samtools as it is under active development and Maq, to my knowledge, is not any more.

However for SNP calling try the mpileup pipeline described on the Samtools webpage instead of using pileup. It's important to generate the I16 info fields for variant calling in our experience.