SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   Bioinformatics (http://seqanswers.com/forums/forumdisplay.php?f=18)
-   -   Problem to use Bedtools after filtering uniquely mapped reads with samtools (http://seqanswers.com/forums/showthread.php?t=16229)

eilosei 12-14-2011 05:11 PM

Problem to use Bedtools after filtering uniquely mapped reads with samtools
 
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 12-15-2011 04:57 AM

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 12-21-2011 04:51 PM

Quote:

Originally Posted by eilosei (Post 59606)
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
Code:

samtools view -H .bam > new.sam
stick on the unique reads
Code:

samtools view .bam | grep ... >> new.sam
convert to bam
Code:

samtools view -Sb -o unique.bam new.sam
also check out samtools reheader.


All times are GMT -8. The time now is 07:39 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.