per bp depth intersected with gff file

Question

[Biomonika (Noolean)]

Hi all,

I have a bam file and gff file and need to report per base coverage for all intervals (features) from gff file. Instead of having 0 coverage for all bps outside those regions, I simply need them not to be reported. I am familiar with bedtools and examined few tools, such as bedtools intersect, bedtools genomecov and bedtools coverage, but I don't seem to find the exact commands (yet). Number of reads per feature won't help me, I truly need per bp report. Thanks a lot:)

[Devon Ryan]

Have you considered just using samtools? Namely: samtools depth -b foo.bed sample.bam? That won't double count overlapping paired end reads and should only print regions specified in the BED file. Remember to convert from 1-based to 0-based coordinates!

Answer 1

[Devon Ryan]

You should be able to use bedtools coverage piped to awk. The formatting of the output from coverage with the -d option may or may not be that useful for you, but you can convert it to whatever you want with awk.