How to quantify reads at a location?

[charlesli]

I have two bam files (H3K56ac ChIP and input) that I want to use to score regions (promoters) that I have defined in a bed file. I've been trying to use seqMINER to do this, but i'm pretty sure there is an easier method? Suggestions for a good tool to do this would be much appreciated. Thanks.

[GenoMax]

Genomecov from Bedtools.

[shi]

Dear Charlesli,

You may try the featureCounts function included in the bioconductor package Rsubread - http://bioconductor.org/packages/rel.../Rsubread.html

Cheers,
Wei

[AlexReynolds]

BEDOPS bedmap is useful for this task. This application takes elements that overlap other elements in UCSC BED files and applies operations to them (counting them, listing their IDs, calculating statistics from score column data, etc.).

For this task, we can use bedmap --count to show the number of reads that overlap some regions of interest (promoters, distal elements, etc.).

First, prepare your regions-of-interest with BEDOPS sort-bed:

$ sort-bed unsortedRegions.bed > sortedRegions.bed

This preparation step only has to be done once for any input to BEDOPS tools. This step allows BEDOPS tools to operate faster and with a lower memory profile than alternatives. Sorting only has to be done once, as BEDOPS tools export data in sorted order.

Next, let's assume your reads are in BAM format. We can use BEDOPS bam2bed to convert this data to sorted BED, and we then pipe it to bedmap to apply the --count operation:

$ bam2bed < reads.bam \
| bedmap --echo --count sortedRegions.bed - \
> answer.bed

The file answer.bed contains results in the following format:

[ region-1 ] | [ count of reads over region-1 ]
[ region-2 ] | [ count of reads over region-2 ]
...
[ region-N ] | [ count of reads over region-N ]

If your regions are promoters, for example, then this result tells you how many reads are contained within each promoter's genomic coordinates, which you specify in sortedRegions.bed.

The default overlap criterion is one base. In other words, one base of overlap between the read and region is sufficient for inclusion. You can adjust this overlap parameter with other options, if you need more stringency.

Note that this result is a sorted BED file, as well, and can be piped to other BEDOPS tools (or other utilities which process BED data).

[charlesli]

bedops looks great. thanks for the detailed explanation. thanks to others for suggestions as well.