bedtools coverage double-counting

Question

[igor]

I am using bedtools coverage to calculate coverage at various regions, but I also want to know how many reads are overlapping any of the regions. Sometimes the regions are close, so the reads span more than one and those reads get counted twice. Is there any way to prevent that or maybe an alternate method that can do that?

[dariober]

Have a look also at htseq-count, it has options to decide what to do with ambiguous overlaps. And also featureCounts (about which I have no experience though).

Answer 1

[harold.smith.tarheel]

You could use bedtools intersect with the -c flag to report the number of overlaps for each read, then change any values >1 to 1 and sum (or filter out the reads with zero overlap and count the number of unfiltered reads).

[dariober]

bedtools intersect -c seems a good strategy to first get reads overlapping one and only one feature (i.e. count = 1). Then you need covergeBed to actually count the (now unambiguous) reads in each feature.