Rnaseq : Overall Statistics
3
0
Entering edit mode
12.2 years ago

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.

rna-seq • 2.1k views
ADD COMMENT
2
Entering edit mode
12.2 years ago

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
ADD COMMENT
2
Entering edit mode
12.2 years ago
Ido Tamir 5.2k
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.

ADD COMMENT
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 1916 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6