average depth of each region in a BED file calculated from a BAM file

Question

[Floydian_slip]

Hi, I have a list of regions in a BED file and I want to know the average depth of each region as calculated from a BAM file, ideally as the fourth column in the output file (like an extra column in the BED file). Is there some combination of Samtools/Bedtools and awk that I can use to get this? Thanks!

Answer 1

[Ram]

I'd recommend taking a look at mosdepth. It has an option where you can give it a BED file and it will calculate depth per interval. Plus, it is a LOT faster than samtools. It doesn't count every single read though (unlike samtools), but that shouldn't matter as it only excludes reads that you shouldn't be counting anyway. Give the manual a good read.

[Floydian_slip]

Thanks for the reply, RamRS! Running into installation issues with mosdepth. Any way either samtools or bedtools can do this for me?

[Ram]

Yeah samtools should be able to do this, but trust me - mosdepth is far quicker. What problems are you running into? If you have access to conda, you should just

conda create -n mosdepth_env -c bioconda mosdepth
conda activate mosdepth_env

[Floydian_slip]

Thanks, RamRS. I was able to install it using docker but ran into error:

$docker run -v ~/Documents/ quay.io/biocontainers/mosdepth:0.2.4--he527e40_0 mosdepth  -n --fast-mode -b ~/Documents/redesign_sorted.bed ~/Documents/sample ~/Documents/ecs-ready.bam
[E::hts_open_format] Failed to open file ~/Documents/ecs-ready.bam
SIGSEGV: Illegal storage access. (Attempt to read from nil?)

I also installed the latest htslib from http://www.htslib.org/download/

Any idea how I can fix this? Thanks so much for your help.

[Ram]

Sorry, I don't know Docker, so I can't say if that error is Docker based or mosdepth based. I recognize the HTSLIB error though. What is the output of

samtools view -H ~/Documents/ecs-ready.bam | head

[Shruteek]

Thank you for mentioning mosdepth. I had never heard of it before. I needed to extract the average per-base coverage of every aligned fragment (of which I had ~25,000,000 per sample) in 10 sequencing run samples. Samtools depth just calculated the coverage of all bases that were covered by any region in the BED file, instead of calculating coverage on each region separately. Samtools bedcov did what I asked, but took hours to run on a single sample. I installed mosdepth and I have no idea how, but mosdepth ran in about 5 minutes per sample and gave me exactly what I needed. I appreciate it!!!

Answer 2

[Floydian_slip]

I ended up using "samtools bedcov" and then writing a quick script to divide the total base pairs (last column) by the length of the region. I matched that up with Samtools depth and it matches. Thanks for your help, RamRS!