This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Read block operation failed with BAM file

I am trying to automate extracting data from an indexed BAM file and I am getting read errors:

[E::bgzf_read_block] Invalid BGZF header at offset 51686517456
[E::bgzf_read] Read block operation failed with error 2 after 0 of 4 bytes

This failure happens after about 12 hours of analysis, near the end of the 52GB BAM file:

$ ll HUDEP.control.DS182418.chr8.bam*
-rw-r--r-- 1 areynolds stamlab 51901713847 Jun 15 12:46 HUDEP.control.DS182418.chr8.bam
-rw-r--r-- 1 areynolds stamlab     7250712 Jun 15 12:59 HUDEP.control.DS182418.chr8.bam.bai

I ran samtools quickcheck on it, which reports no issues:

$ samtools quickcheck -vvv HUDEP.control.DS182418.chr8.bam
verbosity set to 3
checking HUDEP.control.DS182418.chr8.bam
opened HUDEP.control.DS182418.chr8.bam
HUDEP.control.DS182418.chr8.bam is sequence data
HUDEP.control.DS182418.chr8.bam has 25 targets in header.
HUDEP.control.DS182418.chr8.bam has good EOF block.

I can otherwise read data from it with manual samtools and pysam commands to random locations within the BAM file.

$ python
Python 3.8.16 (default, Mar  2 2023, 03:21:46) 
[GCC 11.2.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pysam
>>> print(pysam.__version__)
0.21.0

Remaking the BAM file might be an option but it is very time-consuming. I'd like to understand what went wrong, so I don't have to remake files again.

Is that the only option or are there other tests and fixes I could apply, before this?

bgzip samtools bam pysam htslib
• 5,491 views • link
updated by John Marshall John Marshall • written by Alex Reynolds Alex Reynolds

Indeed, though it is a seek issue when accessing files over HTTP and other network protocols. So, unless there's something you didn't mention and this question is not in fact about accessing files locally, that is unlikely to be the answer to this question.

• 0 views • link
written by John Marshall John Marshall

3 answers

“Invalid BGZF header” means that the index has sent the BGZF reader to a location that is not the start of a BGZF block, so the index is incorrect or not in sync with the BAM file. Hence:

  • When was the index made? It's common to get this error when using an old .bai index file with a newly regenerated BAM file. Try nuking the .bai file and remaking it with samtools index.

  • How was the index made? There have previously been (and currently are, when multithreading) bugs with the on-the-fly index writing. So if you made the index via --write-index --threads …, try nuking the .bai file and remaking it with samtools index.

• 0 views • link
written by John Marshall John Marshall

Good to know, thanks. I am making the BAM file with samtools merge and the index with samtools index, immediately after the BAM file is generated. So far as I can tell, neither case applies.

• 0 views • link
written by Alex Reynolds Alex Reynolds

other tests and fixes I

test the bam is correctly gzip-ped ?

gunzip -t in.bam && echo OK
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Both files test okay:

$ gunzip -t HUDEP.treatment.DS182417.chr8.bam && echo OK
OK
$ gunzip -t HUDEP.control.DS182418.chr8.bam && echo OK
OK

Weird.

• 0 views • link
written by Alex Reynolds Alex Reynolds

this tool may be useful: bamrescue

• 0 views • link
written by size_t • size_t

Log in to answer this question.

More posts like this