Why is an indel clearly visible in IGV, but not reported?

[id0]

I have a BAM file where there are several indels that look real when viewed in IGV, but for some reason they are not being called by any of the variant callers that I tried. Even using samtools mpileup with bcftools does not produce them. The coverage is over 1000x and the indels are fairly frequent, so this is not based on just a few poor quality reads.

I know they are there. I can see them. All I want is the frequency. Why is this so difficult?

[chadn737]

1000x plus is extremely high coverage.

Are these highly repetitive regions, or do you have this coverage for the entire genome? Extremely read depth (like 1000) is typical of repetitive regions and leads to a lot of false positives. So variants with super read depths are often filtered out.

It would help if you posted exactly what you commands/arguments you used in calling variants. For instance, if you followed the examples given in the Samtools manual, you probably filtered out everything with a read depth of over 100x.

[id0]

This is targeted resequencing, so high coverage is expected. However, this is indeed a repetitive region.

The command I used for mpileup (adjusting the depth cutoff):
samtools mpileup -uf hg19.fa -d 100000 -l interval.bed sample.bam | bcftools view -vcg - > sample.vcf

Here are the expanded and the collapsed IGV screenshots of one deletion:

[chadn737]

Is this really "one deletion"? This is obviously a highly repetitive sequence its easy to see that many of these may be misaligned. I don't think you can be so confident looking at these alignments in IGV and say that there is an indel and that the problem is with the variant callers.

[id0]

I meant a deletion in one specific location. Yes, it does seem like there are multiple deletions there. However, the reads reach the end of the repetitive region on both ends. Doesn't that add more confidence?

Regardless, if they are showing up in the alignment as different, why can't they be quantified?

[xied75]

Is this really "one deletion"? This is obviously a highly repetitive sequence its easy to see that many of these may be misaligned. I don't think you can be so confident looking at these alignments in IGV and say that there is an indel and that the problem is with the variant callers.

Does a local re-alignment help situation like this?

[id0]

Does a local re-alignment help situation like this?

I tried with and without realigning. It does not make a significant difference.

[chadn737]

Can you show a picture of this same region in expanded view, but slightly larger so that we can see more of the sequence to the right?

Here's the thing. I could take that very bottom read in the expanded view introduce a gap in the same spot all the gaps start and move it over 6 bps and it would match perfectly and now you would have a 6 bp deletion. Its not coincidence that these gaps appear once you hit the start of that repetitive region.

[id0]

Can you show a picture of this same region in expanded view, but slightly larger so that we can see more of the sequence to the right?

Here's the thing. I could take that very bottom read in the expanded view introduce a gap in the same spot all the gaps start and move it over 6 bps and it would match perfectly and now you would have a 6 bp deletion. Its not coincidence that these gaps appear once you hit the start of that repetitive region.

Here it is...

[abi]

I am having same trouble with Samtools/mpileup

Did you get answer to this ? Can anyone comment ?

The BAM files for alignment map were sorted and indexed.

Then, here is what I did for pileup for processing for a small region of interest:

samtools mpileup -r 0:4,000,679-5,000,894 -uf e.coli.fa sorted.bam > sorted.region.mpileup

Then to obtain SNPs and InDels from mpileup file, here is what I did:

bcftools view -bvcg contigalignsorted.region.mpileup > sorted.region.bcf

bcftools view sorted.region.bcf > sorted.region.vcf


However, when i see my alignment of BAM file on IGV for this region, I notice that the final VCF file obtained through commands above, has a lot of SNPs and InDels missing! I can show a snapshot of a picture of SNP/InDel missed if required.


What went wrong in the calls ?

[abi]

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: