Question: Why Does Mpileup Skip My Mutation ?
7
gravatar for Pierre Lindenbaum
9.2 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum130k wrote:

(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
(...)
ADD COMMENTlink modified 8.2 years ago • written 9.2 years ago by Pierre Lindenbaum130k

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

ADD REPLYlink written 9.2 years ago by Aaronquinlan11k

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

ADD REPLYlink written 9.2 years ago by Pierre Lindenbaum130k

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

ADD REPLYlink modified 12 months ago by RamRS30k • written 9.2 years ago by Pierre Lindenbaum130k

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

ADD REPLYlink written 9.2 years ago by Rm8.0k
5
gravatar for Pierre Lindenbaum
9.2 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum130k wrote:

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
ADD COMMENTlink written 9.2 years ago by Pierre Lindenbaum130k
3

Try "mpileup -E"

ADD REPLYlink written 9.2 years ago by lh332k
3

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

ADD REPLYlink written 8.8 years ago by Pablo Marin-Garcia1.8k
1

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.

ADD REPLYlink written 9.2 years ago by Aaronquinlan11k

will do, thanks.

ADD REPLYlink written 9.2 years ago by Aaronquinlan11k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2044 users visited in the last hour