samtools coverage for WES using bed file

[amyhouseman]

Hello,

I was wondering if it's possible to add both bam and bed files as input for samtools coverage?

Thanks,
Amy

[whitwham]

I am not sure I understand your request. Could you expand on what you want and, more importantly, why you want it?

Understanding your problem would help us know how to fix the issue.

[amyhouseman]

Hello,

I have whole exome sequences, I have aligned them to hg38. I want to work out the % coverage across the reference genome at the points that the bed file states were sequenced. I may be wrong but samtools coverage seems to provide these results in your example samtools coverage -r chr1:1M-12M input.bam

However, rather than just a part of the chromosome I would like to do this using the regions of the exome capture stated in the bed file.

I may be misunderstanding the use of samtools coverage, so apologies if that is the case. I have used samtools depth to work out the average depth but would like some sort of coverage across the parts of the whole exome regions to go alongside.

Sorry if I am wrong, and if this doesn't make sense. I am still learning. Thank you for your reply! Amy

[whitwham]

Now I understand. You are right, samtools coverage only accepts one region where samtools depth can also accept bed files.

It should be possible to add a bed file option to samtools coverage.

[jkbonfield]

There's also samtools bedcov which is explicitly designed for that scenario, but the output format is a little different. Does that suffice?

[amyhouseman]

That's great, I'll have a go at bedcov, thank you!

[amyhouseman]

Hello, bedcov worked well - is there a way to get some similar stats that samtools coverage outputs like:

Number of reads:
Covered bases:
Percent covered:
Mean coverage:
Mean baseQ:
Mean mapQ:

from the bedcov output?
and also, what the columns of the bedcov file mean.

Thanks,
Amy

[nicolo-tellini]

Hello,

I think that what @amyhouseman proposed could be a really great improvement over running the command window by window.

Thanks for your availability
Nicolò

[amyhouseman]

I actually ended up using qualimap for this rather than bedcov, it gave me more stats and output results that were useful

[codebuzer]

Hello Everyone,
I am trying get the gene wise coverage for that I use both samtool(command line) and pysam (python wrapper), everything is good like when I align the results of both the approach covered bases and covered percent are same but a problem occur when covered bases are more than this (end - start + 1 ), and because of that out covered percentage increased even more than 100% (while using pysam) but for the same region samtool using command line gives different like less than 100% or 100% but not more than 100%. I wanted to know how samtool handle these type of situation when at the time of alignment number of covered bases are more than the number of region itself. I am new to this so it might be possible that my question could be wrong.
I will be thankful if anyone here will explain the algorithmic approach of samtool to calculate the coverage especially when covered bases are more than point of the interested region.

[nicolo-tellini]

Hi,

is there any advance on that?

thanks

[whitwham]

It is on my to-do but I have not got around to doing anything about it.