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
You probably want the
-split
option...Thanks, but I doesn't seem to change anything. In the
features.bed
I only specify the start and stop coordinate for the counting:Should this matter?
Why don't you use ggbio to plot coverage?
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.