bcftools doesn't call suspected indel

[rwmills]

Hello-

I am looking for CRISPR-mediated indels using samtools/bcftools 1.2, but bcftools call doesn't report an apparent indel.

If I use this command:
samtools mpileup -guf ref.fasta sortedReads.bam | bcftools call -c - > out.vcf

and inspect the output vcf, I find in the expected indel location that only ref exists:

Human	3078	.	A	.	283.236	.	DP=574;VDB=0.00187095;SGB=-0.590765;RPB=0.623577;MQB=0.639386;BQB=0.0765187;MQ0F=0;AF1=0;AC1=0;DP4=531,0,5,0;MQ=40;FQ=-281.989;PV4=1,0.0104987,1,1	GT:PL	0/0:0
Human	3079	.	A	.	283.236	.	DP=578;VDB=0.958818;SGB=-0.616816;RPB=0.337168;MQB=0.757481;BQB=0.00621451;MQ0F=0;AF1=0;AC1=0;DP4=524,0,6,0;MQ=40;FQ=-281.989;PV4=1,4.05418e-05,1,2.10095e-14	GT:PL	0/0:0

However, if I backup a step and use
samtools mpileup -f ref.fasta sortedReads.bam > out.pileup

and inspect the output pileup file I find the below, which I believe indicates that site 3078 is reference, but 3079 has a 3bp indel. If so, what options should I pass to the caller to report these indels when I setup my pipeline (I don't want to manually inspect every pileup!)?

Human	3078	A	537	*............................................................................................G......................+3ATT..................................................................................+3ATT................................................................................+3TTT..+3TTTCGGG............................................................................................................+3ATT...............+3ATT..................................................................................................................................	HHHHE?H;GGAHHE>ECF/FHC/GHHHHHGHAHHHHHHF/FHHH/H/HHHEHEHHHHHGHHHGCHHG?HHHHHH?HHHHHHGGHHHHGFHHHEEGH0FEAHHGG>HHHHHHHAGE/?GHHHHHHF/<67HHHHHHHGEHHHG/HHHHHHH/>HCHHHGHHHGH>HHHHGFHHAGGHHHHEH//H9H9HFGHCHH/GH:BHHFG/CEGG/?/:F>EFFFEHH/FHHGFGG/GH/HFHHHHH/HGEHHHHHGHHHHHHHHGGHH//G?CFEFFHHH9<85;=/GG3H?E3HBH?HGFHFH/GHG?HGHGHHFCHHH/HHECHFEHHHHEHHGHHFA/GEFFAH/EHHH/HHHEEGGFAGH/EHF/HHHHHGHEGGHH?HEH//GHHG?HFFHH:GHFGGHFH/EHHGH:HGF/FHHHHHHHGHGEHGCHFHHHEHH/HHHHHHHHHHHGHHHHHHHHHHHHHHHHHHHHHHHHFHHHHHHHHHHHHHAHGHHHHHHHHHFHHHHHHHHHHHHHHHHHHHGHHGHHHHGHGHHGHC03B4
Human	3079	A	531	*.+3TTG.+3TTG.+3TTG..+3TTG.+3TCG...+3TTG.+3TTG..+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG...+3TTG.+3TTG.+3TCG.+3TTG...+3TTG.+3TTG....+3TTG.+3TTG...+3TTG.+3TTG.+3TTG.+3TTG.+3TTG...+3TTG.+3TTG..+3TCG.+3TTG.+3TTG...+3TTG.....+3TTG.+3TTG....+3TTG..+3TTG.+3TTG...+3TTG.+3TTG.......+3TTG.+3TTG..+3TTG..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG...+3TTG.+3TTG..+3TTG.+3TTG.+3TTG..+3TTG..+3TTG.+3TTG...+3TTG...+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG...+3TTG...+3TTG..+3TTG.+3TTG.+3ATG.+3ATG..............+3CTG.+3CTG.............................-5CACAT....+3GTG...+3TAG.+3TAG..+3TCG.+3TCG..+3TCG...+3TGG..+3TGG....+3TGG.+2TG....+3TTC.+3TTC..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTTGG+3TTGT+3TTG.+3TTG...+3TTG..+3TGG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG...+3TTG.+3TTG....+3TTG..+3TTG.+3TTG.....+3TTG..+3TTG.+3TTG...+3TTG.+3TTG..+3TTG....+3TTG.+3TTG.....+3TTG.+3TTG..+3TTG...+3TTG..+3TTG.+3TTG...+3TTG.+3TTG.+3TTG.+3TAG.+3TTG.......+3TTG.+3TTG.+3TTG..+3TTG.+3TTG...+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG...+3TTG.....+3TTG.+3TTG.+3TTG.+3TTG..+3TGG.+3TTG..+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG..+3TTG.+3TTG..+3TCG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG..+3TTG.+3TTG..+3TTG..+3TTG...+3TTG.+3TTG.+3TTG.+3TCG.+3TTG.+3TTG.+3TTG..+3TTG...+3TCG.+3TTG.+3TTG.+3TTG.+3TTG....+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.-5CACAT.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG..+3TTG.+3TTG.+3TTG.+3TTG.+3TTG.+3TTG..+3TTGCC.+3TTGC	HH;;/H/FH9/E//E/AFEHF;6HHHHG/HGHHHHHHGBGHHH?;GHEEEHHHGHHHHHHHH;HEHGEHHHHH/HHHHHHG/HHFF//CHG///G?EE9HHHG;HHHHHHHEH/F/HHH=HG/2G723HHHHEHH/EHHHA>H>EEFHH/EH>HHHHHHHCH/GHHHEEHH2///;/H46/=H6H6/D/;FCCB;52/G//5;5/9=//75=99///>D/HHEGHH>HE/H>HHHHH6HEAHHEHH6HHH6HHHHE/=6C/G8/D9//F;;901/GG3G8?GH;H9HH9GEHF=HH/HHHFHG=FHHH/HHF/H5/HGFHHHHHHH:AC/FFFEH/?HHH?HH;6HHHEEH/64HH//HHFHHEH/HHGFF;HEC4;;H6=HEEHHHHFEE;GH>FH;GH;;//CAHHH;HH/HEA;G?C;HG;HHD/;:E;H;;;;;;;;;;;=;;;;;;;;;;;<;;;;;;H=<;;;;;;H;;;;;2;;2;;;<;;=;;>;;;<;;;;;;;;;;;;;;;;;;;;;;H<;<;;<G;10:1

[Richard Haigh]

Hi rwmills

I am having similar problems with bcftools call - I am using old test data which I know contains 60 snps and 5 indels and Varscan has no problem identifying indels from the data but bcftools will only see the snps whatever modifiers I add to the command. Did you ever resolve this?

Richard

[rwmills]

Hi rwmills

I am having similar problems with bcftools call - I am using old test data which I know contains 60 snps and 5 indels and Varscan has no problem identifying indels from the data but bcftools will only see the snps whatever modifiers I add to the command. Did you ever resolve this?

Richard

Hi Richard-

I was never able to get bcftools to make this call even with substantial data grooming before passing to the caller. In my case, it is likely due to the algorithm being confounded by extreme read depth (>8000), as VarScan's overview mentions. VarScan easily calls this insertion as well as a deletion I hadn't noticed before.

[Richard Haigh]

Hi Richard-

I was never able to get bcftools to make this call even with substantial data grooming before passing to the caller. In my case, it is likely due to the algorithm being confounded by extreme read depth (>8000), as VarScan's overview mentions. VarScan easily calls this insertion as well as a deletion I hadn't noticed before.

Thanks for the reply. I have also tried with data with various read depths and all no good - I'm starting to think that the samtools 1.2 version of bcftools call is just broken for indels. In the end I gave up and have managed to bodge my pipeline using the newer version of Varscan (3.2.9) which gives all the right indels but now has a more acceptable vcf output. A simple bit of editing with awk and the files went into snpEff with no problem.