Dear All,

I am interested in counting the number of mapped reads in equal size intervals on genome (say 200 bp), my question is whether featurecounts any option for doing this. According to what I have found so far it sounds like it is more appropriate for counting reads mapping to exons, genes, and promoters. In order to use it for this purpose I define a set of features like following:

library("Rsubread")
library("BSgenome.Hsapiens.UCSC.hg19")
lengths = seqlengths(Hsapiens)
chrlengths = lengths[1:24]
intervalnum <- ceiling(chrlengths/200)

start <- rapply(as.list(intervalnum),function(x) seq(1,x*200,200))
end <- rapply(as.list(chrlengths),function(x) c(as.list(seq(200,x,200)),x))
chr <- rep(names(intervalnum),intervalnum)
ann <- data.frame( 
GeneID=as.character(1:sum(intervalnum)), 
Chr=chr, 
Start=start, 
End=end, 
Strand = rep("+",sum(intervalnum)),
stringsAsFactors=FALSE) 

bams <- dir(".","_sorted.bam$", full=TRUE)
fc <- featureCounts(files=bams, nthreads=5, readExtension3=150, annot.ext=ann)

Is there a better and faster way? (this way, it takes too long for each bam file)

An alternative might be to use windowCounts in the csaw package. You can tinker around with the width, spacing and ext parameters to get it to do what you want; in your case, it seems that values of 200 bp would be appropriate for the first two parameters, and 150 bp for ext. I don't think it would be appreciably faster than featureCounts though, especially since it isn't parallelized; but is speed really a problem? Run it once, leave it for 10-20 minutes depending on your computer's performance/number of BAM files/number of reads, and save the count matrix (from featureCounts) or the SummarizedExperiment object (from windowCounts) to file for future use. It's not like you have to run it every time you open up your data set.

P.S. featureCounts will count whatever features you give it, so it's completely valid to give it fixed intervals. However, if you do so, note that reads spanning the interval boundaries will not be assigned to either interval by default; you'll have to set allowMultiOverlaps=TRUE for such reads to be counted, which is especially important for small intervals. This isn't a problem for gene-based features as genes are usually sparsely distributed (in mammalian genomes, at least).