Samtools / Bcftools Missing Random Obvious Single Base Indels
0
3
Entering edit mode
10.5 years ago
Torst ▴ 980

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
samtools bcftools indel variant-calling • 3.5k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

@Torst, have you tried any other callers? Just curious if they do call this...

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 2698 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6