I'm trying to use bedtools intersect and group by to get a count of reads by region they overlap. This should be straightforward, but the output I'm getting is just a number (presumably, the count of all reads that overlapped any region).

Process:

bamToBed mybam.bam > mybam.bed
sort -k1,1 -k2,2n -k3,3n mybam.bed > sort.mybam.bed
bedtools intersect -a sort.mybam.bed -b regions.bed -wa -wb | groupBy -g 1-3 -c 8 -o count > groups.bed

As far as I can tell from the documentation, that should do just what I need. I receive no errors when running this, but no columns print in the output. I understand that the input files must be sorted on the columns to be output, but I think I've done that.

Any thoughts?