Understanding bedtools coverage for datasets with in-dels.

Question

[Shakunthala Natarajan]

Hello. I have long reads (ONT) and I have mapped them to a reference genome. I have the sorted bam file in hand. As a next step, I wanted to calculate the read coverage depth per position using bedtools with the command:

bedtools coverage -d -split -ibam sorted_bam_file > output_file

I want to further understand how bedtools deals with insertions and deletions, since my reads are long reads (ONT).

Can someone please help me with how to proceed further to understand how in-dels are dealt with during coverage calculation with bedtools as mentioned above?

Thank you!

[GenoMax]

You may want to try pandepth https://github.com/HuiyangYu/PanDepth for this. You probably don't need per-base coverage since the files would become gigantic.

mosdepth can generate per-base coverages and has similar similar functionality: https://github.com/brentp/mosdepth

Both these programs would likely be much faster than bedtools.

[Shakunthala Natarajan]

Hi thank you for bringing this to my notice! I would try this. But presently, I am using bedtools as a part of a script that I am running. To understand the output of the script's working, I need to understand how coverage is being calculated by bedtools when in-dels are present. Essentially, I want to back track from my results, get to the details of how coverage calculation is done by bedtools in the presence of in-dels. Is there some way in which I could do this?

[Pierre Lindenbaum]

or you can use samtools depth in.bam

[Shakunthala Natarajan]

Hi thank you for suggesting samtools. But I would like to understand the coverage calculation by bedtools in the presence of in-dels as mentioned above.

Answer 1

[Brian Bushnell]

You can calculate positional coverage, including indels, using BBTools, like this:

pileup.sh in=mapped.bam basecov=basecov.txt

It has various options that you can see by running the shellscript with no arguments, but the most salient is "delcoverage=t" which is the default. That means bases covered by a deletion event will have their coverage incremented. You can set it to false, and then deletion events will not have their coverage incremented. There's no similar option for insertions, because those bases don't exist in the reference, so obviously it makes no sense to classify them as covered based on an insertion.

This does not require a sorted bam file, it works fine with an unsorted sam file too.

[Shakunthala Natarajan]

Hello Brian! Thank you for suggesting BBTools. I will try this!