samtools sam to bam header

[mattbawn]

Hi

I am trying to sort my bam file using samtools 0.1.19-44428cd
with the following command:

samtools sort test_sequence.bam test_sorted

gives:
[bam_header_read] invalid BAM binary header (this is not a BAM file).
Segmentation fault (core dumped)

The bam file was generated in samtools with :
samtools view -bS test_sequence.sam > test_sequence.bam

and the header in the sam file (generated with bwa) is

@SQ SN:chr10 LN:135534747
@SQ SN:chr11 LN:135006516
@SQ SN:chr11_gl000202_random LN:40103
@SQ SN:chr12 LN:133851895
@SQ SN:chr13 LN:115169878
@SQ SN:chr14 LN:107349540
@SQ SN:chr15 LN:102531392
@SQ SN:chr16 LN:90354753
@SQ SN:chr17_ctg5_hap1 LN:1680828
@SQ SN:chr17 LN:81195210

I also tried:

samtools view -bT wg.fa test_sequence.sam > test_sequence.bam

but got the same error.

I am just starting out in this area, and any help would be appreciated.

Thanks,

[wolma]

Hi,
for question specific to samtools this mailing list is a better place to ask:

[email protected]

You need to subscribe to it at:

https://lists.sourceforge.net/lists/listinfo/samtools-help

before you can post.
---------
Regarding your actual question, the sam header you've posted is certainly not valid.
You can actually use samtools view with the -H switch to show just the header of a sam or bam file.
This should also give an error in your case.

[swbarnes2]

Did you look at the options in samtools view? Try

samtools view -bSh test_sequence.sam > test_sequence.bam

[mattbawn]

[I]

Did you look at the options in samtools view? Try

samtools view -bSh test_sequence.sam > test_sequence.bam

[I/]

Hi, Yes I get the same result.

[dpryan]

It's likely that the file somehow got corrupted. Just delete it and remake it.

[wolma]

the header as you've posted it -

@SQ SN:chr10 LN:135534747
@SQ SN:chr11 LN:135006516
@SQ SN:chr11_gl000202_random LN:40103
@SQ SN:chr12 LN:133851895
@SQ SN:chr13 LN:115169878
@SQ SN:chr14 LN:107349540
@SQ SN:chr15 LN:102531392
@SQ SN:chr16 LN:90354753
@SQ SN:chr17_ctg5_hap1 LN:1680828
@SQ SN:chr17 LN:81195210

lacks a @HD header line and @RG read group lines, so I wouldn't expect this to work with samtools.
I don't know why you would get such an incomplete header with bwa though.

Did you provide the full sam header and , if so, how did you generate it ?

[mattbawn]

the header as you've posted it -

@SQ SN:chr10 LN:135534747
@SQ SN:chr11 LN:135006516
@SQ SN:chr11_gl000202_random LN:40103
@SQ SN:chr12 LN:133851895
@SQ SN:chr13 LN:115169878
@SQ SN:chr14 LN:107349540
@SQ SN:chr15 LN:102531392
@SQ SN:chr16 LN:90354753
@SQ SN:chr17_ctg5_hap1 LN:1680828
@SQ SN:chr17 LN:81195210

lacks a @HD header line and @RG read group lines, so I wouldn't expect this to work with samtools.
I don't know why you would get such an incomplete header with bwa though.

Did you provide the full sam header and , if so, how did you generate it ?

I used the bwa command:

bwa mem -aM -t 8 hg19bwaidx eG2013-072pf_RR_FAM-SJDM-3604.fastq > testX_sequence.sam

with and without the -aM flags. and I also tried bwa samse and obtained the same header.

[dpryan]

The @HD line isn't actually required, though it's normally there. Read groups (@RG lines) are also not required and, in fact, are normally only used if you're calling variants or want to keep some metainfo around. In fact, you don't actually need a header in a BAM file at all, though that only makes sense with unmapped reads.

The "[bam_header_read] invalid BAM binary header (this is not a BAM file)." only occurs if the file actually isn't a BAM file, since what causes this is if the first 4 decompressed bytes aren't equal to "BAM\001" or if samtools can't read up to 4 bytes (or if it somehow reads more, though I'm not sure how that'd happen since it's just reading into a buffer and returning up to the first 4 bytes). Often this error occurs if the file somehow got corrupted or is actually a SAM file with the wrong extension.

[mattbawn]

I just found the following on github:

https://github.com/samtools/samtools/issues/42

Your command line takes advantage of GNU extensions to the behaviour of getopt() that allow options and arguments (e.g. filenames) to be intermingled. Mac OS X provides the BSD getopt() which implements only the standard behaviour with all options preceding the arguments.

The easy solution is to write your samtools command portably, with all the options at the start:

samtools view -bS -o file.bam file.sam

The alternative would be to compile samtools yourself using GNU getopt. However this is likely more bother than it's worth, and you would end up with a samtools command that behaves differently (albeit more conveniently!) from the other commands on your OS X machine.

i therefore retried with the following:

samtools view -bS -o testX_sequence.bam testX_sequence.sam

which at the moment seems to have worked.

[dpryan]

That's actually not the source of your problem, since "samtools view -bS test_sequence.sam > test_sequence.bam" wouldn't run into that issue. Remaking the BAM file just deletes it and replaces it with a non-corrupt version, as I had previously suggested.

[wokai001]

The error message "[bam_header_read] invalid BAM binary header (this is not a BAM file)"
is generated when "samopen" calls "bam_header_read". When bam_header_read retrieves the first four bytes from the file and either gets less than 4 Bytes or gets somthing different from "BAM\001", this kind of error is thrown.

Seems as if the file is corrupt.

[wolma]

but if you were right, then why did ALL the previous attempts fail with the same error ?
mattbawn says he tried all of these:

samtools view -bS test_sequence.sam > test_sequence.bam
samtools view -bT wg.fa test_sequence.sam > test_sequence.bam
samtools view -bSh test_sequence.sam > test_sequence.bam

but now
samtools view -bS -o testX_sequence.bam testX_sequence.sam
seems to work.

The only difference I can see is that this time there is no shell-dependent redirection involved.

P.S.: you are right about the header specification. thanks for pointing that out.

[dpryan]

@wolma: To whom are you replying?

My guess is that there's a bad sector on the drive that kept getting used. That's more likely than finding a broken shell these days.

BTW, a quick way to check if a BAM file has any chance of being valid is:

zcat foo.bam | head -c 4 | od -c

If you don't get the following, then samtools will fail:

0000000   B   A   M 001

[wolma]

personally, I favor the less cryptic python version:

import gzip
with gzip.open('test.bam','rb') as i:
    print(i.read(4))

which prints:
b'BAM\x01'

but leaving that aside the appropriate test would then be to just try to write to yet another file using redirection, then to the same file using -o .
I haven't looked into the source code, but potentially it could also have to do with how samtools handles writing to files as opposed to writing to stdout. Maybe the later fails under some exotic combination of compiler/OS/zlib version ?

[dpryan]

I haven't looked into the source code, but potentially it could also have to do with how samtools handles writing to files as opposed to writing to stdout. Maybe the later fails under some exotic combination of compiler/OS/zlib version ?

Perhaps, it's always tough to rule out obscure issues like that. The underlying C code doesn't do anything particularly weird there, so I'd expect a lot of things to break in that case (but who knows).