Why does mpileup skip my mutation ?

[lindenb]

(cross posted on biostar: http://biostar.stackexchange.com/questions/9961 )

Hi all,

I know, by capillary sequencing, that one of my samples contains a mutation at position:120458492.

Some reads were aligned with bwa and I can clearly see my mutation using samtools tview.

       120458491 120458501
CCTGCTCTGGGGAGCTATGCCAGGATGGGTGCC
........R........................
.....C..A..................C.....
........C. ......................
............  ...................
........A.....    ...............
,,,,,,,,a,,,,,,,,,,,,,,,,  ......
,,,,,,,,a,,,,,,,,,,,,,,,,,,,    .
........A........................
........A........................
.............................C...
........A........................
.................................
........A........................
.................................
.................................
........A........................
......................C..........
........A........................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
........A........................
.................................
 ................................
  ...............................
   .....A........................
      ,,,,,,,,,gg,,,,,,,,,,,t,,,,
       .A........................
        .........................

however, we I call the mutations using samtools mpileup

${SAMTOOLS}/samtools mpileup -uf  ${HG19} sample.bam |\
  ${SAMTOOLS}/bcftools/bcftools view  -bvcg - > snps.bcf
${SAMTOOLS}/bcftools/bcftools view  snps.bcf | gzip --best > snps.vcf.gz

But I can't see the mutation in the VCF.

here is the output of pileup (not mpileup):

(...)
chr1    120458492   G   26  C.AaaAA.A.A..A.A,A...A,AA^].    JddfhQacfhhgggahehhhhaB[hg
(...)

How can I know why mpileup (or bcftools ?) skipped this mutation ?

Thanks !

Pierre

[swbarnes2]

Try redoing mpileup with -Buf. That's worked for me. Like you, I observed that my sanger-verified mutation looked fine in pileup, but when I looked at it in the mpileup run without -B, the Baq calculations destroyed the quality scores of that locus, so it wouldn't call a SNP. Running it with -B should make your mpileup look like your pileup at that locus, and then the SNP will have high enough quality to be counted.

[lindenb]

That worked ! many thanks !

[lh3]

try -E
...