How to retrive the number of reads aligned to various genomic features?

Question

[AlexAbdulkaderKheirallah]

Hello All,

I've got a ChIPseq dataset and want to generate a pie chart outlining proportion of various genomic features with which a factor associates. I came across various applications that can do it automatically such as ChIPseek or ChIPpeakAnno. However these programs take BED file with outlined peak locations and compares it with genomic annotation while what I want to do is to use a BED file containing information about genome annotation and query ChIPseq BAM file to retrieve a number of reads that align to those locations. I suppose I should be able to do it with samtools but not sure how to do it. Could you help? Last but not least ChIPseek or ChIPpeakAnno do not give me information on features I'm interested in such as CGIs or enhancer regions.

Thanks!

[AlexAbdulkaderKheirallah]

I was thinking a bit about James Ashmore's approach and I wonder whether in principle to you advise doing that or it is rather risky to use all the available reads for genome occupancy and whether genome occupancy should be assessed based on peak calls?

Answer 1

[James Ashmore]

This can be done using the bedtools coverage function. Given a BED file containing the regions you want to examine, and a BAM file of aligned reads, use the following command to output a tab delimited file which lists each region followed by the number of reads aligned to that region:

bedtools coverage -counts -a regions.bed -b sample.bam > reads.in.regions.bed

[AlexAbdulkaderKheirallah]

Thanks James