Question: coverageBed: what am I doing wrong?
0
gravatar for jon.brate
3.1 years ago by
jon.brate200
Norway
jon.brate200 wrote:

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

coverage bam bedtools • 1.9k views
ADD COMMENTlink modified 3.1 years ago by Aaronquinlan10k • written 3.1 years ago by jon.brate200
1

You probably want the -split option...

ADD REPLYlink written 3.1 years ago by Devon Ryan89k

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 REPLYlink written 3.1 years ago by jon.brate200

Why don't you use ggbio to plot coverage?

ADD REPLYlink written 3.1 years ago by Irsan6.8k

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 REPLYlink written 3.1 years ago by jon.brate200
2
gravatar for Devon Ryan
3.1 years ago by
Devon Ryan89k
Freiburg, Germany
Devon Ryan89k wrote:
bedtools coverage -d -split -abam SRR_1.420.bam -b features.bed > SRR_420.cov

This produces the expected results.

ADD COMMENTlink written 3.1 years ago by Devon Ryan89k

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 REPLYlink modified 3.1 years ago • written 3.1 years ago by jon.brate200

It looks like I have 2.23.0 installed locally.

ADD REPLYlink written 3.1 years ago by Devon Ryan89k

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

ADD REPLYlink written 3.1 years ago by Devon Ryan89k

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

ADD REPLYlink written 3.1 years ago by jon.brate200
1

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 REPLYlink written 3.1 years ago by Devon Ryan89k

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

ADD REPLYlink written 3.1 years ago by Devon Ryan89k

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

ADD REPLYlink written 3.1 years ago by jon.brate200
2
gravatar for Aaronquinlan
3.1 years ago by
Aaronquinlan10k
United States
Aaronquinlan10k wrote:

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 COMMENTlink modified 3.1 years ago by Devon Ryan89k • written 3.1 years ago by Aaronquinlan10k
0
gravatar for harold.smith.tarheel
3.1 years ago by
United States
harold.smith.tarheel4.3k wrote:

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 COMMENTlink modified 3.1 years ago • written 3.1 years ago by harold.smith.tarheel4.3k

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 REPLYlink written 3.1 years ago by jon.brate200
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1090 users visited in the last hour