Problem to use Bedtools after filtering uniquely mapped reads with samtools

[eilosei]

Hi there,

I've aligned my ChIP-seq data with BWA and filtered for uniquely mapped reads with this command:

samtools view file.bam | grep "XT:A:U" > file.unique.bam

Then I tried to convert the file.unique.bam to .bed file with Bedtools but failed. The command I used is:

bamToBed -i file.unique.bam > file.bed

And I got error message like:

BgzfStream ERROR: read block failed - invalid block header
BamHeader ERROR: could not read magic number
BamReader ERROR: Could not load header data for file.unique.bam

Does anyone have this problem before? How can I solve it?

Thanks!

[arvid]

There are two reasons why it doesn't work:

1. "file.unique.bam" is not in the BAM format (look at "samtools view" for correct usage).
2. You stripped the SAM headers with your grep command.

You might need find a better way of filtering for uniquely mapped reads which preserves the SAM headers, or you add the SAM headers back to the header-stripped SAM file you created with your first command afterwards (look at "samtools view").

Then, convert it to BAM with SAMtools and feed it into BEDTools.

[jbrwn]

Originally Posted by eilosei

Hi there,

I've aligned my ChIP-seq data with BWA and filtered for uniquely mapped reads with this command:

samtools view file.bam | grep "XT:A:U" > file.unique.bam

Then I tried to convert the file.unique.bam to .bed file with Bedtools but failed. The command I used is:

bamToBed -i file.unique.bam > file.bed

And I got error message like:

BgzfStream ERROR: read block failed - invalid block header
BamHeader ERROR: could not read magic number
BamReader ERROR: Could not load header data for file.unique.bam

Does anyone have this problem before? How can I solve it?

Thanks!

grab header

samtools view -H .bam > new.sam

stick on the unique reads

samtools view .bam | grep ... >> new.sam

convert to bam

samtools view -Sb -o unique.bam new.sam

also check out samtools reheader.