Why Does Mpileup Skip My Mutation ?

Question

[Pierre Lindenbaum]

(cross posted on seqanswers: http://seqanswers.com/forums/showthread.php?t=12542 )

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.

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

Thanks !

Pierre

EDIT:

...and here is the output of pileup (not mpileup):

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

[Aaronquinlan]

my first guess is strand bias. it appears that the alternate allele is consistently on one strand and not the other.

[Pierre Lindenbaum]

Nice observation. How could I tell mpileup to ignore this parameter ?

[Pierre Lindenbaum]

Got an answer on seqanswers: http://seqanswers.com/forums/showthread.php?p=45786#post45786. I'll try it tomorrow.

[Rm]

BTW did you check quality of changed nucleotides ? if low they are might not be reported.

Answer 1

[Pierre Lindenbaum]

Here is the solution , suggested by swbarnes2 on seqanswers:

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.

$ gunzip -c snps.vcf.gz | grep chr1 | grep "120458492"
chr1    120458492    .    G    A    99    .    DP=26;AF1=0.5;CI95=0.5,0.5;DP4=10,2,12,2;MQ=60;PV4=1,0.37,1,1    PL:GT:GQ    255,0,255:0/1:99

[lh3]

Try "mpileup -E"

[Pablo Marin-Garcia]

to see the difference between -B and -E (the recommended option by lh3) see lh3 answer here: http://biostar.stackexchange.com/questions/10038/mpileup-variant-call

[Aaronquinlan]

Ah yes, I've been burnt by this as well. BAQ certainly reduces false positives, but I have also found hundreds of true positives that it misses.

[Aaronquinlan]

will do, thanks.