I want to count the number of reads that fall in some predefined regions of a genome. Normally, this can be easily done with bedtools multicov. I have sorted and indexed bam files of my alignments. I created a BED file with my regions of interest (excel, then saved as tab-delimited, and added .bed extension).

Here's what BED file look like:

Chr1    500000  505001  
Chr1    1000000 1001001 
Chr1    2000000 2000251 
Chr1    2500000 2500101

I run bedtools with this command:

bedtools multicov -bams ./*.bam -bed bedfile.bed

I have 5 bam files in the same directory, as well as 5 corresponding .bai index files. I get read counts for all 5 bam files but only for the last feature of the bed file, whatever the last line is. If I add more features it will still give me the counts for the last one only. Any idea what's wrong?

Thanks, TP