Rnaseq : Overall Statistics

Question

[Nicolas Rosewick]

Hi,

Very simple question : Is there a tool to count the number of reads that are aligning in an annotated region in the genome. I've a gff file and a bam file (from tophat - paired-end reads 76x2). I thought maybe using bedtools..

Thanks,

N.

Answer 1

[Istvan Albert]

Apparently the -L flag of samtools that can also do this:

 -L FILE  output alignments overlapping the input BED FILE [null]

So you could do a:

$ ~/bin/samtools view -c results.bam chrI:1-100000
180    
$ cat query.bed 
chrI    0    100000
$ ~/bin/samtools view -c -L query.bed results.bam 
180

Answer 2

[Ido Tamir]

bedtools coverage -counts -abam the.bam -b the.gff

Although the PE reads complicate this a little bit, I would say. If this is for counting for genes, I would suggest you look at something like HtSeq-count http://www-huber.embl.de/users/anders/HTSeq/doc/count.html

bedtools AFAIR does not take split reads in the 'abam' file into account.

[Madelaine Gogol]

If you use the -split argument, it will treat split bam entries as distinct intervals (if that's what you want to do).