sort error: "found in binary header but not text header."

[map2085]

samtools 1.3:

samtools sort -T sort.tmp -@ 1 -O bam -o out.bam in.bam

where in.bam has 28 million reads and 94647 reference sequences in the header (all Ensembl mouse transcripts). This is just standard RNA-seq.

Output:

[bam_sort_core] merging from 11 files...
...
[E::trans_tbl_add_sq] @sq SN (ENSMUST00000164489) found in binary header but not text header.
....
[E::hts_idx_push] unsorted positions

samtools reports an error message as above ("found in binary header but not text header") for every single reference sequence in the BAM header - in this case, all 94647 Ensembl mouse transcripts.

This does NOT happen in samtools-1.0

This error is thrown for BAMs consisting of only paired-end reads, and for BAMs consisting of only single-end reads.

[daviesrob]

The code that makes this check has changed quite a lot since samtools-1.0.
Would it be possible to make a cut-down version of your input bam file that we could take a look at?

[map2085]

I also have many small-yield experiments (< 1million reads) for which the same code executed smoothly without error.

For my other datasets with many reads (>20 million), this error happens every time, to every dataset.

When i downsample the bam file from one of my large datasets to 10 million reads, sort works fine, without error, in version 1.3

Thus, it seems to be a problem on large BAMs (> 20 million reads) with large numbers of reference headers (90k mouse transcripts). And it happens with every such BAM file, as per my experience.

[daviesrob]

That's a pity, it would be easier if the files were small. Are any of your datasets available for download?

[map2085]

still unpublished, unfortunately.

i've reproduced the problem with a publicly avaialble Ecoli dataset. BAM size = 1.3 Gb. how can i send this?

[map2085]

SRR794835 available from ENA. no read processing necessary.
aligned to the transcriptome fasta attached.
expression calculated with RSEM.
attempted to sort transcript.bam

Escherichia_coli_str_k_12_substr_mg1655.GCA_000005845.2.26.dna.genome.fa--Escherichia_coli_str_k_12_substr_mg1655.GCA_000005845.2.26.gtf.rsem.mrna.transcripts.fa.zip

[daviesrob]

How did you align it and make the bam file?

[jmarshall]

I've reproduced this with a BAM file with no text @SQ headers, so a fix will be forthcoming that will hopefully fix this for your files too.

Such BAM files are valid but uncommon. What tool did you use to produce your .bam file?

[map2085]

I aligned the reads with GEM mapper, nothing fancy. I sorted this by coordinate successfully.

Then I passed the bam through RSEM (rsem-calculate-expression). This gives a BAM file with a few extra fields in each alignment in the BAM file.

I confirmed that:

1. I successfully samtools sort (by coordinate) the GEM BAM file,
2. but samtools sort failed (as described above) to sort the RSEM BAM file. These files have the same exact alignment entries (bar the extra user-tags in the RSEM bam), albeit in different order.

The headers are almost identical:

samtools view -H gem.bam

@hd VN:1.3 SO:coordinate
@sq SN:AAC73112 LN:66
@sq SN:AAC73113 LN:2463
@sq SN:AAC73114 LN:933
.....
@sq SN:b4708 LN:174
@sq SN:b4712 LN:84
@sq SN:b4713 LN:84
@rg ID:0 PG:GEM PL:ILLUMINA SM:0
@pg ID:GEM PN:gem-2-sam VN:1.423

samtools view -H rsem.bam

@hd VN:1.4 SO:unknown
@pg ID:RSEM
@sq SN:AAC73112 LN:66
@sq SN:AAC73113 LN:2463
@sq SN:AAC73114 LN:933
....
@sq SN:b4708 LN:174
@sq SN:b4712 LN:84
@sq SN:b4713 LN:84

gem.header.txt
rsem.header.txt

(So the rsem bam is "missing" a Read Group tag in the header?)

When I put the RSEM header onto the GEM bam file and do samtools sort, I get the following error:

[bam_sort_core] merging from 7 files...
[bam_translate] RG tag "0" on read "SRR794835.19965709" encountered with no corresponding entry in header, tag lost. Unknown tags are only reported once per input file for each tag ID.
[bam_translate] RG tag "0" on read "SRR794835.9783492" encountered with no corresponding entry in header, tag lost. Unknown tags are only reported once per input file for each tag ID.
...

[jmarshall]

RSEM writes BAM files with mostly no text headers (or so it appears looking at the source code briefly), which produces the effects you've seen.

When you output SAM with samtools view, samtools appends basic @SQ headers if there aren't already any(!). Thus your samtools view -H rsem.bam shows a file without explicit textual @SQ headers (as those shown are synthetic ones appended last) and without @RG headers (leading to the tag lost problems you saw when you grafted the RSEM header on to the other file).

So the fix in f63a282 should have fixed this for your files.

[iamjli]

For those looking for a workaround in v1.3, the following worked for me:

samtools view -H file.bam | samtools reheader - file.bam | samtools sort - -o file.sorted.bam