arq5x / bedtools

A powerful toolset for genome arithmetic.
http://code.google.com/p/bedtools/
GNU General Public License v2.0
140 stars 85 forks source link

coverageBed output incorrect with -abam option, correct with -abam stdin #55

Closed nkrumm closed 11 years ago

nkrumm commented 11 years ago

I cant seem to get correct output from coverageBed when providing the bam file directly:

bedfile.bed (51,011 lines)

...
Example1Contig  912832681   912833713   MyFeature
...

using samtools view on bamfile.bam shows no reads in this interval:

$ samtools view bamfile.bam Example1Contig:912832681-912833713
$ 

Piping this directly into coverageBed reports zero:

$ samtools view bamfile.bam Example1Contig:912832681-912833713 | 
coverageBed -abam stdin -b bedfile.bed | grep MyFeature
Example1Contig  912832681   912833713   MyFeature   0   0   1032    0.0000000

As well as piping the entire bam file into coverageBed:

$ samtools view bamfile.bam | coverageBed -abam stdin -b bedfile.bed | grep MyFeature
Example1Contig  912832681   912833713   MyFeature   0   0   1032    0.0000000

BUT, using the -abam <file> option, a totally different count is registered:

$ coverageBed -abam bamfile.bam -b bedfile.bed | grep MyFeature
Example1Contig  912832681   912833713   MyFeature   1355    1032    1032    1.0000000

Any ideas as to what is going on (or what I missed here!) are appreciated. Thanks, Nik

arq5x commented 11 years ago

Hi Nik, it looks like you are piping SAM to bedtools, not BAM. Try using -b or -u (preferably):

$ samtools view    -u    bamfile.bam Example1Contig:912832681-912833713 | 
coverageBed -abam stdin -b bedfile.bed | grep MyFeature