I have very noisy nanopore reads and am trying to call SNPs/indels.
I'm runinig into some trouble when truing to use the
samtools mpileup | bcftools call combination.
It seems to incorrectly be calling very long indels where it looks like there is no support. For example:
samtools mpileup -f ref.fa sample.bam -r Chromosome:198940-198940 produces:
Chromosome 198940 C 37 .,.-1A.-1A...-1A.-1A.-1A.,-1a.-1A.,-1a,-1a.,.-1A,,.-1A,,+1a,+1a,.-1A,,-1a.,-1a,-1a,,-1a,-1a.-1A,-1a. 06?5/<:8>/4<?3691465<3354.08:23:11@3C
This looks fine... however when I try to call the genotype with bcftools I get a very long indel.
samtools mpileup -ugf ref.fa sample.bam -r Chromosome:198940-198940 -t AD | bcftools call -mv produces:
Chromosome 198940 . CAAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA C 225 . INDEL;IDV=32;IMF=0.5;DP=62;VDB=0.145385;SGB=-0.662043;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,7,2;MQ=60 GT:PL:AD 1/1:255,27,0:0,9
htsbox pileup -f ref.fa sample.bam | grep 198940 it shows there is one read with this indels together with many other reads with smaller indels.
Chromosome 198940 C C,C-75AAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA,C-4AAAA,C-2AA,C-1A,T,C+1A,C+1T,C+3GGT 0/4:25,1,1,2,24,5,2,1,1
Is there any reason why it would be grouping these reads together into the big indel?
Any idea what is going wrong here?