How to put header back to .bam file?

[eilosei]

I've filtered my .bam file for uniquely mapped reads with this command:

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

Then I wish to convert the .unique.bam file 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.sorted.unique.bam


It seems with the grep command I left out the header of the bam file. Can anyone tell me how to fix this problem? Thank you!

[swbarnes2]

throw a -h into "samtools view", but I won't vouch for how that will work with the grep command embedded in there.

[krobison]

throw a -h into "samtools view", but I won't vouch for how that will work with the grep command embedded in there.

It most certainly won't behave; the grep will strip out the header. A bit challenging (at least for me) to write a grep expression to save both, but with perl/awk/etc it is straightforward. The appropriate perl to replace the grep would be

perl -n -e 'print $_ if (/^\@/ || /XT:A:U/) '

There is also the reheader command in samtools

samtools view -H file.bam > header.sam  # extract header only
samtools reheader header.sam file.unique.bam

It might even be possible to pipe the first command to the second, but I haven't tested that

[volks]

samtools view -H file.bam > header.sam  # extract header only
samtools reheader header.sam file.unique.bam

It might even be possible to pipe the first command to the second, but I haven't tested that

try:

samtools view -H file.bam > header.sam
samtools view file.bam | grep "XT:A:U" | cat header.sam - | samtools view -Sb - > file.unique.bam

[eilosei]

try:

samtools view -H file.bam > header.sam
samtools view file.bam | grep "XT:A:U" | cat header.sam - | samtools view -Sb - > file.unique.bam

This code works! Thank you!

[eilosei]

Thank you for the suggestion. But the reheader command came out gibberish. What's the problem?


There is also the reheader command in samtools

samtools view -H file.bam > header.sam  # extract header only
samtools reheader header.sam file.unique.bam

It might even be possible to pipe the first command to the second, but I haven't tested that[/QUOTE]

[sophiespo]

Thank you for the suggestion. But the reheader command came out gibberish. What's the problem?


There is also the reheader command in samtools

samtools view -H file.bam > header.sam  # extract header only
samtools reheader header.sam file.unique.bam

It might even be possible to pipe the first command to the second, but I haven't tested that

[/QUOTE]


Sorry, this is an old question and it might be solved by now but for anyone who is searching, samtools will print to STDOUT unless you direct it somewhere.

samtools view -H file.bam > header.sam  # extract header only
samtools reheader header.sam file.unique.bam > file.new.bam