samtools
failed to detect one mutation
--Hello,
i have run bcftools mpileup with: bcftools mpileup --threads 16 -C 50 -E -Ou -f genome.fa mybam.bam | bcftools call --threads 16 -cv -Oz -o mybam.vcf.gz
but it failed on a single mutation position, see here IGV panel with the mutation:
the corresponding sam file on this position:
samtools view mybam.bam chr21:43104255-43104428 | grep "43104346" M02792:126:000000000-K9RHL:1:1114:12346:24983 163 chr21 43104253 0 147M = 43104346 240 TTAACTGTCTTTGAAAAGAACATGAAGTTTTTATAATTTACATGAAAAAAAGGCAAACAAACCTGGCTAAACGTCGGTTTATTGTGCAACCGAGAGCACCTGTCTCCATGACGACATGCTCCAATTTTGAAATAAAATGAACAGTTG BEDHBFFDGFEEDEGGHHFIDECKFJGEFGFEBEBJEFEBDEEKGJJJJHGHDIEJJEEJJDEIKIFGBJJDBEHBHDFFBEEKFGFFJDEBFHFHGEFHHKEHHHHEELFDBFFEDLFHIHEJEFFFKFJKFBJJJDKGKFEGEFD XA:Z:chr21,+6495931,147M,0; ZA:Z:CGCT ZB:Z:CCGT BC:Z:5 MD:Z:147 RG:Z:SAD190004_S5 NM:i:0 AS:i:147 XS:i:147 QX:Z:CCB ABC RX:Z:CGC-CCG xc:i:1 zd:Z: xm:i:1 M02792:126:000000000-K9RHL:1:1114:12346:24983 83 chr21 43104346 0 147M = 43104253 -240 GAGCACCTGTCTCCATGACGACATGCTCCAATTTTGAAATAAAATGAACAGTTGACTCTGTAAGGGAAAATGAGAGCTGATTATTTTGCTGGGAAGATATCAAACACATGGAATATGTCAGCAGCATGACATACACTATCAAATTAC HHFLEHGEFFFGIKFEGDEHEKEEFGGIKGEIJJEHEFEBFFGEEHFFKHEJEHEGGGEDBGHHHHFGGEFHHGIFHFHEJCEJJJEFHFIHHFIHEBEHKGFFLEKEEHHGEBEEEHLHGLHFLEEHEDDCDCDHCDEHDECHCCD XA:Z:chr21,-6496024,147M,0; ZA:Z:CGCT ZB:Z:CCGT BC:Z:5 MD:Z:147 RG:Z:SAD190004_S5 NM:i:0 AS:i:147 XS:i:147 QX:Z:CCB ABC RX:Z:CGC-CCG xc:i:1 zd:Z: xm:i:1
i don't know which parameter to use to recover the mutation.
any ideas ?
If you're able to compile things, it would be interesting to know if https://github.com/samtools/bcftools/pull/1679 cures this. I believe it may under-count AD in rare circumstances though, which I think is why it never got merged. I don't know the specifics on it sadly, nor how frequent a problem it was (compared to the rate of missing mutations for example).
There is also bcftools mpileup--indels-2.0, but be warned at the moment it will greatly increase the number of false positive calls because it doesn't yet have appropriate filtering enabled (WIP).
i think i found the problem, i have a XA tags in sam file and it's a secondary alignment reported by bwa-mem, does 'bcftools call' discard secondary alignment in default behaviour ?
That's correct, bcftools mpileup skips secondary alignments by default. This can be controlled with the option
--skip-any-set STR|INT Skip reads with any of the bits set [UNMAP,SECONDARY,QCFAIL,DUP]
