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

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 • 5.4k views
ADD COMMENT
1
Entering edit mode

You probably want the -split option...

ADD REPLY
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?

ADD REPLY
0
Entering edit mode

Why don't you use ggbio to plot coverage?

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

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

This produces the expected results.

ADD COMMENT
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
...
ADD REPLY
0
Entering edit mode

It looks like I have 2.23.0 installed locally.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
8.7 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])

enter image description here

ADD COMMENT
0
Entering edit mode
8.7 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.]

ADD COMMENT
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: enter image description here

ADD REPLY

Login before adding your answer.

Traffic: 1615 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6