coverageBed: what am I doing wrong?
3
0
Entering edit mode
5.1 years ago
jon.brate ▴ 250

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

bedtools coverage bam • 2.8k views
1
Entering edit mode

You probably want the -split option...

0
Entering edit mode

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?

0
Entering edit mode

Why don't you use ggbio to plot coverage?

0
Entering edit mode

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.

2
Entering edit mode
5.1 years ago
bedtools coverage -d -split -abam SRR_1.420.bam -b features.bed > SRR_420.cov


This produces the expected results.

0
Entering edit mode

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

0
Entering edit mode

It looks like I have 2.23.0 installed locally.

0
Entering edit mode

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

0
Entering edit mode

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

1
Entering edit mode

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.

0
Entering edit mode

0
Entering edit mode

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

2
Entering edit mode
5.1 years ago

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


0
Entering edit mode
5.1 years ago

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

0
Entering edit mode

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: