(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
(...)
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
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.
will do, thanks.
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
Log in to answer this question.
my first guess is strand bias. it appears that the alternate allele is consistently on one strand and not the other.
Nice observation. How could I tell mpileup to ignore this parameter ?
Got an answer on seqanswers: http://seqanswers.com/forums/showthread.php?p=45786#post45786. I'll try it tomorrow.
BTW did you check quality of changed nucleotides ? if low they are might not be reported.