This is a test version of Biostars. For the public version, visit https://www.biostars.org.
coverageBed: what am I doing wrong?

I want to plot the coverage across a region of a chromosome from a .bam file. I thought this was pretty straight forward, but for some reason I am not getting the correct coverage using bedtools coverage.

Here is what the .bam-file looks like in IGV: IGV screen shot

Notice that there are almost no reads mapping to the introns. And when I calculate the coverage using: coverageBed -d -b SRR_1.420.bam -a features.bed > SRR_420.cov
I get this (NB the x-axis coordinates are wrong):
coverageBed

Any Idea what is happening? Link to the input files

bedtools coverage bam
• 7,027 views • link
updated by Aaronquinlan Aaronquinlan • written by jon.brate jon.brate

You probably want the -split option...

• 1 views • link
written by Devon Ryan Devon Ryan

Thanks, but I doesn't seem to change anything. In the features.bed I only specify the start and stop coordinate for the counting:

Supercontig_1.420   29000   41000   AGO_miRNA2

Should this matter?

• 1 views • link
written by jon.brate jon.brate

Why don't you use ggbio to plot coverage?

• 1 views • link
written by Irsan Irsan

I didn't know of this package. I will check it out, thanks! I used Gviz first, but found it a little too complicated for my needs. Also it is difficult to use to plot stranded coverage.

• 1 views • link
written by jon.brate jon.brate

3 answers 

bedtools coverage -d -split -abam SRR_1.420.bam -b features.bed > SRR_420.cov

This produces the expected results.

• 1 views • link
written by Devon Ryan Devon Ryan

Sorry, but I am a little confused. Which version of bedtools are you using? I am using version 2.24 and apparently coverageBed changed: "As of version 2.24.0, the coverage tool has changed such that the coverage is computed for the A file, not the B file."

Therefore I used: coverageBed -d -b SRR_1.420.bam -a features.bed -split > SRR_420.cov which gives the weird graph. I tried your command, it gives the following output which I don't fully understand (and the output is over 100MB):

Supercontig_1.420 29925 30020 SRR350207.1113251/1 50  + 29925 30020 0,0,0   1 95,  0,   1   1
Supercontig_1.420 29925 30020 SRR350207.1113251/1 50  + 29925 30020 0,0,0   1 95,  0,   2   1
Supercontig_1.420 29925 30020 SRR350207.1113251/1 50  + 29925 30020 0,0,0   1 95,  0,   3   1
Supercontig_1.420 29925 30020 SRR350207.1113251/1 50  + 29925 30020 0,0,0   1 95,  0,   4   1
...
• 1 views • link
written by jon.brate jon.brate

It looks like I have 2.23.0 installed locally.

• 1 views • link
written by Devon Ryan Devon Ryan

The -split option seems to not apply in more recent versions.

• 1 views • link
written by Devon Ryan Devon Ryan

That is weird. I had also version 2.20 installed and your suggestion worked fine!

• 1 views • link
written by jon.brate jon.brate

Yeah, I just tried in 2.25.0 and that doesn't work correctly unless one instead does bedtools coverage -a features.bed -b <(bedtools bamtobed -split -i SRR_1.420.bam) -d > SRR_420.cov. I'm guessing that something changed in 2.24.0 that broke some things.

• 1 views • link
written by Devon Ryan Devon Ryan

Good to know. I will keep my eyes open for a bug fix. Thanks for the help

• 1 views • link
written by jon.brate jon.brate

Thanks very much for reporting this. I have pushed a fix to the bedtools repository. It will be a part of the next release, but you can always update your code now by cloning from the repository.

Below is an updated plot based upon the following command:

$ bedtools coverage -a features.bed -b SRR_1.420.bam -d -split > features.bed.split.cov

#R
> s <- read.table('features.bed.split.cov')
> plot(s[,6])

enter image description here

• 1 views • link
updated by Devon Ryan Devon Ryan • written by Aaronquinlan Aaronquinlan

Why do you think that the coverage plot is wrong? It looks like the IGV peaks at 34 kb and 39 kb, albeit at base-pair resolution (which you specified with the -d flag).

[Edit: Sorry, I totally missed the issue regarding introns.]

• 1 views • link
written by harold.smith.tarheel harold.smith.tarheel

To me it doesn't look like the bedtools coverage is correct on the introns. I don't show the introns and exons here, but they correspond perfectly to the coverage peaks in the IGV screen shot. Here I fixed the x-axis on the coverage plot from bedtools: enter image description here

• 1 views • link
written by jon.brate jon.brate

Log in to answer this question.

More posts like this