coverageBed: what am I doing wrong?

Question

[jon.brate]

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:

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):

Any Idea what is happening? Link to the input files

[Devon Ryan]

You probably want the -split option...

[jon.brate]

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?

[Irsan]

Why don't you use ggbio to plot coverage?

[jon.brate]

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.

Answer 1

[Devon Ryan]

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

This produces the expected results.

[jon.brate]

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

[Devon Ryan]

It looks like I have 2.23.0 installed locally.

[Devon Ryan]

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

[jon.brate]

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

[Devon Ryan]

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.

[Devon Ryan]

It turns out that there's already a bug report about this.

[jon.brate]

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

Answer 2

[Aaronquinlan]

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])

Answer 3

[harold.smith.tarheel]

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.]

[jon.brate]

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: