indel calling problems with mpileup/bcftools

[alexwong.carleton]

Hello,

I'm using samtools and bcftools to call snps from a targeted resequencing project. Nucleotide polymorphisms seem to be fine, but indels are proving to be a problem. I can see the indels in, say, samtools tview, with high coverage. However, my samtools/bcftools output find only a few reads at the relevant position.

My commands (following local realignment) are:

samtools mpileup -uf concatenated.fasta -L 1000000 -d 1000000 Ec_A10.realigned.bam >Ec_A10.mpileup
bcftools view -bvc Ec_A10.mpileup >Ec_A10.var.raw.bcf
bcftools view Ec_A10.var.raw.bcf >Ec_A10_name.vars.tsv

The relevant line from the output is:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Ec_A10.test
concatenated 2000 . CATCCTCGATGT C 185 . INDEL;DP=9;VDB=0.0074;AF1=0.5;AC1=1;DP4=2,2,1,4;MQ=60;FQ=159;PV4=0.52,0.00017,1,1 PL 223,0,194

So, what is strange is that the reported depth is 9 reads, but there are clearly many more - I've also attached a screen shot from tview. Any suggestions on how to clear this up?

Thanks!

[abi]

The results are wrong sometimes

More importantly , it even shows up wrong SNPs which otherwise is not a SNP as per the bam file loaded on IGV.


0 4009334 . A C 125 . DP=5;VDB=3.277706e-02;RPB=-1.291774e+00;AF1=0.5;AC1=1;DP4=1,1,2,1;MQ=59;FQ=74;PV4=1,0,1,0.38 GT:PL:GQ 0/1:155,0,101:99
0 4009336 . T A 17.1 . DP=5;VDB=6.080000e-02;RPB=1.291774e+00;AF1=0.5;AC1=1;DP4=2,1,1,1;MQ=59;FQ=20.1;PV4=1,1,0.14,1 GT:PL:GQ 0/1:47,0,154:50

As above 4009334 is a SNP but not as per IGV picture below:

[swbarnes2]

Hello,

I'm using samtools and bcftools to call snps from a targeted resequencing project. Nucleotide polymorphisms seem to be fine, but indels are proving to be a problem. I can see the indels in, say, samtools tview, with high coverage. However, my samtools/bcftools output find only a few reads at the relevant position.

My commands (following local realignment) are:

samtools mpileup -uf concatenated.fasta -L 1000000 -d 1000000 Ec_A10.realigned.bam >Ec_A10.mpileup
bcftools view -bvc Ec_A10.mpileup >Ec_A10.var.raw.bcf
bcftools view Ec_A10.var.raw.bcf >Ec_A10_name.vars.tsv

The relevant line from the output is:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Ec_A10.test
concatenated 2000 . CATCCTCGATGT C 185 . INDEL;DP=9;VDB=0.0074;AF1=0.5;AC1=1;DP4=2,2,1,4;MQ=60;FQ=159;PV4=0.52,0.00017,1,1 PL 223,0,194

So, what is strange is that the reported depth is 9 reads, but there are clearly many more - I've also attached a screen shot from tview. Any suggestions on how to clear this up?

Thanks!

I only count 5 reads covering base 2000 in your picture.

Why don't you try counting the fastq directly by grepping a little bit of sequence from each allele, or by adding the putative indel allele to your reference, and realigning?

[alexwong.carleton]

swbarnes2 - Sorry, I should have made clear - the image captures only a small proportion of the reads covering the region. The mutation I'm looking for is an 11-bp deletion.

I'll try out both of your suggestions - thanks! I should think that one or both should work, although ideally it would be nice to not have to have solutions tailored for specific mutations.