Disclaimer: I may use coverage and 'mean read depth' interchangeably in this post. I'm refering to the average, per-base read depth.
I'm running and compaing some mean coverage estimates for some specific bed regions on my bam files using bedtools; however, I'm having trouble finding the correct way to do so while taking into consideration the positions that were not covered by any read.
Earlier, I wanted to obtain the mean coverage on my bam files for the whole genome; and that seemed to go rather smoothly using bedtools genomecov ; it has both '-d' and '-dz' parameters to estimate the coverage either with 1-based or 0-based coordinates respectively (the first takes into account all bases, whereas the second only non-zero positions).
Later I saw a lot of topics saying that, for specific regions, the better tool would be bedtools coverage, and that seemed to work rather well to obtain 1-based coordinates read depth (using the -d parameter), but I can't find a way to obtain an estimate that doesn't discard zero (non-covered) positions... none of the parameters from the manual seem to do it (at least not in a straightforward way it seems), which seems odd seen as bedtools genomecov does have it.
Do you know of any other way to get that using bedtools?
Second disclaimer: I am aware that other software is able to do it. I've already successfully tried samtools depth. However, my advisor asked me to compare those results with what you get with bedtools specifically... thank you in advance.
