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?