How to work out coverage (100X) of WES using exome bed file?

Question

[amy__]

Hello,

I know that there are lots of coverage questions on here but quite a few are 10 years + and really confusing.

I have WES fastq files, I have mapped these to hg38 and have their bam files. I also have the bed files for the Exome panel used in the sequencing.

I wish to work out the coverage (like 100X). I know that bedtools has a way but I am finding it hard to work out.

Do I need to get a list of the gene locations even though I have the exon region bed files?

Thanks! Amy

Answer 1

[ATpoint]

Adapting from Tools To Calculate Average Coverage For A Bam File?:

samtools view -u -L exome.bed your.bam \
| samtools depth -a \
|  awk '{sum+=$3} END { print "Average = ",sum/NR}'

The first one subsets for alignment covering the exome (use the bed file with the targets of the capture), and the rest then calculated the average coverage.

[Pierre Lindenbaum]

samtools depth can use a bedfile

samtools depth -b exome.bed -a in.bam |  awk '{sum+=$3} END { print "Average = ",sum/NR}'

[amy__]

Thank you so much, I will try this!

[amy__]

Hi ATpoint, do you know what the variable $3 is? Thanks! Amy

[ATpoint]

What do you mean? The third column in the depth output (that is what $3) means is what we want to sum up. Does that make sense?

[amy__]

Hi, yes you are right, I'm sorry! Slurm was messing with the command! Thanks for all your help

Answer 2

[Pierre Lindenbaum]

https://github.com/brentp/mosdepth

mosdepth  -by exome.bed  --no-per-base output in.bam

look into the summary Total region to get the mean depth