How to get deletion count from bam file in IGV?

[nottsbio]

I am trying to get a count of all bases, insertions and deletions from a bam file viewed in IGV genome viewer for each individual nucleotide. I can get the counts for the bases as a wig file (example shown below):

track type=wiggle_0
#Columns: Pos, Combined Strands A, Combined Strands C, Combined Strands G, Combined Strands T, Combined Strands N
variableStep chrom=E17F span=1
1 0 241 0 0 0
2 337 0 0 0 0
3 414 0 1 1 0
4 6 412 8 5 0
5 4 411 30 5 0
6 475 6 2 0 0
7 9 1 417 69 0

However, I also really want to get the count of deletions for each position. It is possible to get this count if you go through each position individually and write down the numbers (as shown in the image). But I was thinking surely there must be a way to do this in IGV? I have tried igvtools but it still doesn't give me a deletion count. Any guidance or other suggestions will be greatly appreciated!

Thanks
-T

[dpryan]

I'd be surprised if IGV included a way to do that, it would almost never be used. You could either parse the output of "samtools mpileup" or just directly use a variant caller, as applicable.

[nottsbio]

I'd be surprised if IGV included a way to do that, it would almost never be used. You could either parse the output of "samtools mpileup" or just directly use a variant caller, as applicable.

Thanks for the suggestions! I'm not very experienced in the bioinformatics, would you be able to clarify what you mean by these? I'm not sure how I would parse the output of mpileup. This is what I am trying with it but I am unsure what options to use to get the output I want.

$ samtools mpileup -f reference.fasta -Q 13 input.sorted.bam -o output.txt

Also directly using a variant caller, does that mean calling variants from my bam file - what sort of output can I get from that?

Sorry for all the questions and thank you for your suggestions, I really appreciate it!!

[dpryan]

Maybe we should go about this in a different way. What is your biological goal? That is, what is the biological question that you're trying to answer by doing this?

[nottsbio]

I have a mutant sample with a deletion and I want to see if the deletion count in this region is significantly higher compared to the wild type (using a new genome sequencer - so I want to see if it can detect deletions). But I want to get the base (ACGT, insertion, deletion) counts for every position so I can tell if the detected deletion is significant or if the sequencing device just has a high error rate. Hope that makes more sense!

[dpryan]

Right, so you want to call variants between the mutant and wild-type samples, paying attention to only a small region. Just use a variant caller and ensure that it finds the deletion in the mutant but not wild-type sample. I'd recommend GATK's haplotype caller, since it enables joint variant calling. If you really wanted, you could compare the genotyping probabilities between the samples (it's one of the fields in the VCF file produced). This assumes that the deletion isn't very large, in which case other tools targeted specifically toward finding large deletions might give better results.

[nottsbio]

Thank you for your help, I will try this!