I am using samtools to call variants in a haploid genome (yes I know it is designed for diploid). It finds SNPs easily, and most of the indels, but there are a few indels it can NOT seem to find no matter what parameters I give it. I use bwa mem for the alignment of MiSeq PE 150 reads, and I get about 150x coverage, all good. Using tview etc, there is clearly a G deletion at position 600, but bcftools will not seem to call it, despite calling some nearby (~100bp before). All the mapping and base qualities seem to be fine.
This is the command I am using:
samtools mpileup -B -g -u -f ref.fa aln.bam | bcftools view -v -
It happens whether I use the -B
or not.
Any help appreciated.
Below is the INDEL position pileup +/- 2 bases:
Contig_Plasmid 598 C 185 ....,,,..,,,.,,,,,,..........,...,.,,,....,,.$,$,.,....,.,,...,,......,,,...,,,,,,,,,,,,,,....,...,,,.,.....,.,,.,....................,,,,,,,,.....,.,.,..,,,,,,
,,,,,,,.,,.....,,,,,,..,.., BHHHCBFHHGGGHGGGGGGHHGHHGDHHHGHHHGHHHHFHHHHHHAHHHHHHGGGHHFHHIGHFHHHHHHHHHHHHGHHHHHHHHHEHD5FH2HH?GHHHGIHHIHHHFHHGGHHHBHHHHGGHHHHHHHHHFHHHHGHHFGHHHHHHHHHHHGHGHHHHHHH?HHBH
FFFFFHHHHHHB@HCDH
Contig_Plasmid 599 A 183 ....,$,$,..,,,.,,,,,,..........,...,.,,,....,,,.,....,.,,...,,......,,,...,,,,,,,,,,,,,,....,...,,,.,.....,.,,.,....................,,,,,,,,.....,.,.,..,,,,,,,,
,,,,,.,,.....,,,,,,..,.., EGHHCBFHHGGFHFGGGGGHHHFHHFHHHGHHHFHHHHHHHHHHHHHHHHFHHHHHHHHFHGHHHFHIHHHHHHFHHHHGHHHHAHGDHH2HH1GHHHHHHHHHHHFHHHHHHHGHHHHHHHHHHHHHHHFHFHHHHH5HHGHHHHHHHHHHHHHHHHGII4HHFHFF
FFFHHHHHHB?HCDH
Contig_Plasmid 600 T 181 .$.$.$.$,.-1G.-1G,-1g,-1g,-1g.-1G,-1g,-1g,-1g,-1g,-1g,-1g.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G,-1g.-1G.-1G.-1G,-1g.-1G,-1g,-1g,-1g.-1G.-1G.-1G.-1G,-1g,-1g,-1
g.-1G,-1g.-1G.-1G.-1G.-1G,-1g.-1G,-1g,-1g.-1G.-1G.-1G,-1g,-1g.-1G.-1G.-1G.-1G.-1G.-1G,-1g,-1g,-1g.-1G.-1G.-1G,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g.-1G.-1G.-1G.-1G,-1g.-1G.-1G.-1G,-1
g,-1g,-1g.-1G,-1g.-1G.-1G.-1G.-1G.-1G,-1g.-1G,-1g,-1g.-1G,-1g.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G.-1G,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g.-1G.-1G.-1G.-1G.-1G,-1g.-1
G,-1g.-1G,-1g.-1G.-1G,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g,-1g.-1G,-1g,-1g.-1G.-1G.-1G.-1G.-1G,-1g,-1g,-1g,-1g,-1g,-1g.-1G.-1G,-1g.-1G.-1G,-1g 4FHEBHHFFFHFGFGGGHHHHHHFHHHGHHHFHHGHHHHHHHHHHHHHHGHHHFHHIFHGFHHHHHIHHHGHFHHHHHHHGH5H55HHDHGFGHHHFHHHHHHHFHHHHGHHFHHHHHHHHHFHGHHHHHFHHHIGGHHGHHHHHHHHHHHGHHEH2HH4HGHHFFFFFHHHHHHBAHDDH
Contig_Plasmid 601 G 179 ,*********************************************************************************************************************************************************************************^]t BHHFFFHFGGGGGHHGHHHHHHHGGHHHFHHHHHHHHHHHHHHHHHHHHHFHHHFHBHHHHGIIHHHHIFHHHHHGHHHHHFDHHDHHFGHHHHHHHHHGHDHHGFHHGBHHHHHGHHHHHHHHHGHFHHHIHHHHHHHHHHHHHHHHHHHHHHHH4HGHHGGGGGHHHGHHFFHCDHB
Contig_Plasmid 602 G 179 ,..,,,.,,,,,,..........,....,.,,,....,,,.,....,.,,...,,......,,,...,,,,,,,,,,,,,,....,...,,,.,.....,.,,.,....................,,,,,,,,.....,.,.,..,,,,,,,,,,,,,.,,.....,,,,,,..,..,, BHHFFFHFGGGGGHHGHHHHHHHGGHHHFHHHHHHHHHHHHHHHHHHHHHFHHHFHBHHHHGIIHHHHIFHHHHHGHHHHHFDHHDHHFGHHHHHHHHHGHDHHGFHHGBHHHHHGHHHHHHHHHGHFHHHIHHHHHHHHHHHHHHHHHHHHHHHH4HGHHGGGGGHHHGHHFFHCDHH
I hope this gets answered, we've occasionally noticed similar situations and could not explain it, left us scratching our heads that some other cutoff/threshold might have been triggered
@Torst, have you tried any other callers? Just curious if they do call this...
I also found this odd behaviour too, which I filed as an Issue on the samtools github site: https://github.com/samtools/samtools/issues/79