This is a test version of Biostars. For the public version, visit https://www.biostars.org.
average depth of each region in a BED file calculated from a BAM file

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!

bed bam coverage ngs
• 5,542 views • link
updated by Shruteek • Shruteek • written by Floydian_slip Floydian_slip

2 answers

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.

• 0 views • link
written by Ram Ram

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

• 0 views • link
written by Floydian_slip Floydian_slip

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
• 0 views • link
written by Ram Ram

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.

• 0 views • link
updated by Ram Ram • written by Floydian_slip Floydian_slip

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
• 0 views • link
written by Ram Ram

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!!!

• 0 views • link
written by Shruteek • Shruteek

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!

• 0 views • link
written by Floydian_slip Floydian_slip

Log in to answer this question.

More posts like this