It's pretty easy with BEDOPS bedmap --count and bam2bed:

$ bedmap --echo --count --delim '\t' regions.bed <(bam2bed < reads.bam) > answer.bed

The file answer.bed will contain regions from the sorted reference file regions.bed, along with the counts of reads from reads.bam that overlap (or map to) those regions in the last column.

One base of overlap is considered when doing counts. If you want something more stringent, use --fraction-map or --fraction-ref or other options, depending on what you need. See bedmap --help or the online docs for more information on overlap criteria.

If you want to do more than --count reads, there are other operators, like --bases, --bases-uniq, --echo-overlap-size, and others, which can be useful for calculating overlap, coverage, and other statistics derived from mapped reads. See the online documentation for a full discussion.

The file regions.bed should be sorted. If you don't know the sort order of your regions of interest, just sort them first with sort-bed:

$ sort-bed regions.unsorted.bed > regions.bed

Then run bedmap as described. Easy-peasy, lemon-squeezy.

Hello,

I think featureCounts is your friend.

fin swimmer

Hello,

just for fun here is another way using gnu parallel, samtools view and awk.

awk -v OFS="\t" '{print $1, $2+1,$3,$4}' genes.bed|parallel --colsep '\t' -k "samtools view your.bam {1}:{2}-{3}|wc -l|awk -v OFS=\"\t\" '{print \"{1}\",int({2})-1,\"{3}\",\"{4}\",\$1}'"

<del>It has just a minor flaw: The coordinates in your bam file are 0-based. samtools view uses 0-based. So your region will start 1 base early. But this might be ok in your case.</del>

fin swimmer