fastx/

bgzf.rs

//! BGZF: the block-compressed gzip variant used across the samtools ecosystem.
//!
//! A BGZF file is an ordinary gzip file — any gzip tool can decompress it — made
//! of independent members of at most 64 KiB each, every one carrying its own
//! compressed size in a `BC` extra field. Because the blocks are independent, a
//! reader that knows where they start can seek to any uncompressed position by
//! jumping to the enclosing block and decompressing just that. This is what makes
//! a 3 GB bgzipped reference genome randomly accessible.
//!
//! [`BgzfReader`] implements [`Read`] and [`Seek`] in *uncompressed* coordinates,
//! so anything generic over those two traits — including [`crate::IndexedFasta`] —
//! works on a bgzipped file without knowing it.
//!
//! ```no_run
//! use fastx::bgzf::{BgzfReader, BgzfWriter};
//! use std::io::{Read, Seek, SeekFrom, Write};
//!
//! // Write a seekable gzip file.
//! let mut writer = BgzfWriter::create("ref.fa.gz")?;
//! writer.write_all(b">chr1\nACGT\n")?;
//! writer.finish()?;
//!
//! // Read 4 bytes starting at uncompressed offset 6.
//! let mut reader = BgzfReader::open("ref.fa.gz")?;
//! reader.seek(SeekFrom::Start(6))?;
//! let mut buf = [0u8; 4];
//! reader.read_exact(&mut buf)?;
//! assert_eq!(&buf, b"ACGT");
//! # Ok::<(), fastx::Error>(())
//! ```

use std::fs::File;
use std::io::{self, BufReader, BufWriter, Read, Seek, SeekFrom, Write};
use std::path::{Path, PathBuf};

use crate::error::{Error, Result};
use crate::format::CompressionLevel;

/// Largest uncompressed payload placed in one block.
///
/// The spec caps a whole block at 64 KiB; htslib uses 0xff00 for the payload so
/// that even incompressible data plus headers stays under the limit.
pub const MAX_BLOCK_PAYLOAD: usize = 0xff00;

/// Fixed part of a BGZF gzip header, up to and including `XLEN`.
const HEADER_LEN: usize = 12;
/// The `BC` extra subfield: `SI1`, `SI2`, `SLEN` and `BSIZE`.
const EXTRA_LEN: usize = 6;
/// The gzip trailer: CRC32 and ISIZE.
const TRAILER_LEN: usize = 8;

/// The 28-byte empty block that marks a complete BGZF file.
///
/// Its presence is how tools tell a truncated file from a finished one, so
/// [`BgzfWriter::finish`] always appends it.
pub const EOF_BLOCK: [u8; 28] = [
    0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff, 0x06, 0x00, 0x42, 0x43, 0x02, 0x00,
    0x1b, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
];

/// Header of one BGZF block, as read from the file.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
struct BlockHeader {
    /// Total bytes the block occupies in the file, `BSIZE + 1`.
    compressed_len: usize,
    /// Bytes of the block that precede the deflate payload.
    payload_offset: usize,
}

/// Read and validate one block header from `bytes`.
fn parse_block_header(bytes: &[u8]) -> Result<BlockHeader> {
    if bytes.len() < HEADER_LEN {
        return Err(bgzf_error("block header is truncated"));
    }
    if bytes[0] != 0x1f || bytes[1] != 0x8b {
        return Err(bgzf_error("not a gzip member"));
    }
    if bytes[2] != 8 {
        return Err(bgzf_error("unsupported compression method"));
    }
    if bytes[3] & 0x04 == 0 {
        return Err(bgzf_error(
            "gzip member has no extra field, so it is gzip but not BGZF",
        ));
    }
    let extra_len = u16::from_le_bytes([bytes[10], bytes[11]]) as usize;
    if bytes.len() < HEADER_LEN + extra_len {
        return Err(bgzf_error("extra field is truncated"));
    }
    let extra = &bytes[HEADER_LEN..HEADER_LEN + extra_len];

    // Walk the subfields looking for `BC`; other subfields are legal.
    let mut cursor = 0;
    while cursor + 4 <= extra.len() {
        let si1 = extra[cursor];
        let si2 = extra[cursor + 1];
        let slen = u16::from_le_bytes([extra[cursor + 2], extra[cursor + 3]]) as usize;
        let value = cursor + 4;
        if value + slen > extra.len() {
            return Err(bgzf_error("extra subfield runs past the extra field"));
        }
        if si1 == b'B' && si2 == b'C' {
            if slen != 2 {
                return Err(bgzf_error("BC subfield is not two bytes"));
            }
            let bsize = u16::from_le_bytes([extra[value], extra[value + 1]]) as usize;
            let compressed_len = bsize + 1;
            let overhead = HEADER_LEN + extra_len + TRAILER_LEN;
            if compressed_len <= overhead {
                return Err(bgzf_error("BSIZE is smaller than the block overhead"));
            }
            return Ok(BlockHeader {
                compressed_len,
                payload_offset: HEADER_LEN + extra_len,
            });
        }
        cursor = value + slen;
    }
    Err(bgzf_error("gzip member has no BC extra subfield"))
}

fn bgzf_error(message: &'static str) -> Error {
    Error::Other(format!("BGZF: {message}"))
}

/// Inflate a raw deflate stream of known uncompressed size.
fn inflate(payload: &[u8], expected: usize, out: &mut Vec<u8>) -> Result<()> {
    out.clear();
    out.reserve(expected);
    let mut decoder = flate2::Decompress::new(false);
    decoder
        .decompress_vec(payload, out, flate2::FlushDecompress::Finish)
        .map_err(|e| Error::Other(format!("BGZF: corrupt block: {e}")))?;
    if out.len() != expected {
        return Err(bgzf_error("block size does not match its ISIZE field"));
    }
    Ok(())
}

/// The `.gzi` index that makes a BGZF file seekable by uncompressed offset.
///
/// The on-disk layout is the one `bgzip --index` writes: a little-endian `u64`
/// count followed by that many `(compressed_offset, uncompressed_offset)` pairs.
/// The first block is implicit — it always sits at `(0, 0)` — so a file of *n*
/// blocks yields *n − 1* entries.
#[derive(Debug, Clone, Default, PartialEq, Eq)]
pub struct GziIndex {
    /// Block starts, including the implicit `(0, 0)` first entry.
    blocks: Vec<BlockOffset>,
    /// Total uncompressed size, known only when the index came from scanning a
    /// file or from the writer that produced it. A `.gzi` on disk records block
    /// starts and nothing else, so a parsed index cannot know where the data
    /// ends — and must not guess, or `SeekFrom::End` would land in the middle of
    /// the last block.
    total_uncompressed: Option<u64>,
}

/// Where one block begins, in both coordinate systems.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct BlockOffset {
    /// Byte offset of the block in the compressed file.
    pub compressed: u64,
    /// Byte offset of the block's first byte in the decompressed stream.
    pub uncompressed: u64,
}

impl GziIndex {
    /// Build an index by walking every block header in the file.
    ///
    /// Only the headers are read, not the payloads, so this is fast: it touches
    /// about 18 bytes per 64 KiB of input.
    pub fn build<R: Read + Seek>(mut reader: R) -> Result<GziIndex> {
        reader.seek(SeekFrom::Start(0))?;
        let mut reader = BufReader::with_capacity(64 * 1024, reader);
        let mut blocks = Vec::new();
        let mut compressed = 0u64;
        let mut uncompressed = 0u64;
        let mut header = [0u8; HEADER_LEN + 64];

        loop {
            if !read_exact_or_eof(&mut reader, &mut header[..HEADER_LEN])? {
                break;
            }
            let extra_len = u16::from_le_bytes([header[10], header[11]]) as usize;
            if HEADER_LEN + extra_len > header.len() {
                return Err(bgzf_error("extra field is implausibly large"));
            }
            reader.read_exact(&mut header[HEADER_LEN..HEADER_LEN + extra_len])?;
            let block = parse_block_header(&header[..HEADER_LEN + extra_len])?;

            // Skip the payload and CRC, then read ISIZE.
            let skip = block.compressed_len - block.payload_offset - 4;
            io::copy(&mut reader.by_ref().take(skip as u64), &mut io::sink())?;
            let mut isize_bytes = [0u8; 4];
            reader.read_exact(&mut isize_bytes)?;
            let payload_len = u32::from_le_bytes(isize_bytes) as u64;

            blocks.push(BlockOffset {
                compressed,
                uncompressed,
            });
            compressed += block.compressed_len as u64;
            uncompressed += payload_len;

            // A zero-length block is the EOF marker; nothing follows it.
            if payload_len == 0 {
                break;
            }
        }
        Ok(GziIndex {
            blocks,
            total_uncompressed: Some(uncompressed),
        })
    }

    /// Build an index for a file on disk.
    pub fn build_from_path<P: AsRef<Path>>(path: P) -> Result<GziIndex> {
        let path = path.as_ref();
        GziIndex::build(
            File::open(path).map_err(|e| {
                Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display())))
            })?,
        )
    }

    /// Parse a `.gzi` file.
    pub fn parse<R: Read>(mut reader: R) -> Result<GziIndex> {
        let mut count_bytes = [0u8; 8];
        reader.read_exact(&mut count_bytes)?;
        let count = u64::from_le_bytes(count_bytes);
        // Guard against a corrupt count asking for a terabyte of allocation.
        if count > 1 << 32 {
            return Err(bgzf_error("index claims an implausible number of blocks"));
        }
        // The first block is implicit.
        let mut blocks = Vec::with_capacity(count as usize + 1);
        blocks.push(BlockOffset {
            compressed: 0,
            uncompressed: 0,
        });
        let mut pair = [0u8; 16];
        for _ in 0..count {
            reader.read_exact(&mut pair)?;
            blocks.push(BlockOffset {
                compressed: u64::from_le_bytes(pair[..8].try_into().expect("8 bytes")),
                uncompressed: u64::from_le_bytes(pair[8..].try_into().expect("8 bytes")),
            });
        }
        Ok(GziIndex {
            blocks,
            // A `.gzi` does not record the total size.
            total_uncompressed: None,
        })
    }

    /// Parse a `.gzi` file from disk.
    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<GziIndex> {
        let path = path.as_ref();
        GziIndex::parse(BufReader::new(File::open(path).map_err(|e| {
            Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display())))
        })?))
    }

    /// Serialise in `bgzip --index` format, omitting the implicit first block.
    pub fn write<W: Write>(&self, out: &mut W) -> Result<()> {
        let count = self.blocks.len().saturating_sub(1) as u64;
        out.write_all(&count.to_le_bytes())?;
        for block in self.blocks.iter().skip(1) {
            out.write_all(&block.compressed.to_le_bytes())?;
            out.write_all(&block.uncompressed.to_le_bytes())?;
        }
        Ok(())
    }

    /// Write the index next to the BGZF file as `<path>.gzi`.
    pub fn write_to_path<P: AsRef<Path>>(&self, bgzf_path: P) -> Result<PathBuf> {
        let target = gzi_path(bgzf_path.as_ref());
        let mut file = BufWriter::new(File::create(&target)?);
        self.write(&mut file)?;
        file.flush()?;
        Ok(target)
    }

    /// Block starts, in file order, including the implicit first one.
    pub fn blocks(&self) -> &[BlockOffset] {
        &self.blocks
    }

    /// Number of blocks, including the EOF marker if the file has one.
    pub fn len(&self) -> usize {
        self.blocks.len()
    }

    /// True when the index describes no blocks at all.
    pub fn is_empty(&self) -> bool {
        self.blocks.is_empty()
    }

    /// Total uncompressed size.
    ///
    /// `None` for an index parsed from a `.gzi` file, which records only where
    /// blocks start. Build the index with [`GziIndex::build`] if you need this —
    /// it is a header-only scan, so it is cheap.
    pub fn uncompressed_len(&self) -> Option<u64> {
        self.total_uncompressed
    }

    /// The block containing uncompressed byte `offset`.
    fn block_for(&self, offset: u64) -> Option<BlockOffset> {
        // The last block whose uncompressed start is <= offset.
        match self
            .blocks
            .binary_search_by(|b| b.uncompressed.cmp(&offset))
        {
            Ok(i) => Some(self.blocks[i]),
            Err(0) => None,
            Err(i) => Some(self.blocks[i - 1]),
        }
    }
}

/// The conventional index path for a BGZF file: `<path>.gzi`.
pub fn gzi_path(bgzf: &Path) -> PathBuf {
    let mut name = bgzf.as_os_str().to_os_string();
    name.push(".gzi");
    PathBuf::from(name)
}

/// True when the first bytes look like a BGZF block rather than plain gzip.
///
/// ```
/// # use fastx::bgzf::{is_bgzf, EOF_BLOCK};
/// assert!(is_bgzf(&EOF_BLOCK));
/// assert!(!is_bgzf(&[0x1f, 0x8b, 0x08, 0x00]));  // gzip without FEXTRA
/// assert!(!is_bgzf(b"ACGT"));
/// ```
pub fn is_bgzf(bytes: &[u8]) -> bool {
    parse_block_header(bytes).is_ok()
}

/// A BGZF reader that seeks in uncompressed coordinates.
///
/// Sequential reads need no index. [`Seek`] does: build one with
/// [`GziIndex::build`] (fast, header-only) or load a `.gzi` file.
pub struct BgzfReader<R: Read + Seek> {
    inner: R,
    index: Option<GziIndex>,
    /// Decompressed contents of the block currently in hand.
    block: Vec<u8>,
    /// Read cursor within `block`.
    block_pos: usize,
    /// Uncompressed offset at which `block` begins.
    block_start: u64,
    /// Compressed offset of the next block to read.
    next_compressed: u64,
    eof: bool,
    /// Scratch buffer for one compressed block.
    raw: Vec<u8>,
}

impl BgzfReader<File> {
    /// Open a BGZF file, loading `<path>.gzi` if it is present.
    ///
    /// Without a `.gzi` the reader still works sequentially, and [`Seek`] will
    /// build an index on first use.
    pub fn open<P: AsRef<Path>>(path: P) -> Result<BgzfReader<File>> {
        let path = path.as_ref();
        let file = File::open(path)
            .map_err(|e| Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display()))))?;
        let index = match GziIndex::from_path(gzi_path(path)) {
            Ok(index) => Some(index),
            Err(Error::Io(e)) if e.kind() == io::ErrorKind::NotFound => None,
            Err(e) => return Err(e),
        };
        let mut reader = BgzfReader::new(file)?;
        reader.index = index;
        Ok(reader)
    }
}

impl<R: Read + Seek> BgzfReader<R> {
    /// Wrap a seekable reader, checking that it really is BGZF.
    pub fn new(mut inner: R) -> Result<BgzfReader<R>> {
        inner.seek(SeekFrom::Start(0))?;
        let mut probe = [0u8; HEADER_LEN + 64];
        let read = read_up_to(&mut inner, &mut probe)?;
        if read == 0 {
            return Err(bgzf_error("file is empty"));
        }
        parse_block_header(&probe[..read])?;
        inner.seek(SeekFrom::Start(0))?;
        Ok(BgzfReader {
            inner,
            index: None,
            block: Vec::new(),
            block_pos: 0,
            block_start: 0,
            next_compressed: 0,
            eof: false,
            raw: Vec::new(),
        })
    }

    /// Attach an index, enabling [`Seek`].
    pub fn with_index(mut self, index: GziIndex) -> Self {
        self.index = Some(index);
        self
    }

    /// The index in use, if any.
    pub fn index(&self) -> Option<&GziIndex> {
        self.index.as_ref()
    }

    /// Current uncompressed position.
    pub fn position(&self) -> u64 {
        self.block_start + self.block_pos as u64
    }

    /// Unwrap the underlying reader.
    pub fn into_inner(self) -> R {
        self.inner
    }

    /// Decompress the block at compressed offset `at`, which becomes current.
    fn load_block_at(&mut self, at: u64, uncompressed_start: u64) -> Result<()> {
        self.inner.seek(SeekFrom::Start(at))?;
        self.next_compressed = at;
        self.block_start = uncompressed_start;
        self.block_pos = 0;
        self.block.clear();
        self.eof = false;
        self.read_next_block()
    }

    /// Read and decompress the block at `self.next_compressed`.
    fn read_next_block(&mut self) -> Result<()> {
        let mut header = [0u8; HEADER_LEN + 64];
        if !read_exact_or_eof(&mut self.inner, &mut header[..HEADER_LEN])? {
            self.eof = true;
            self.block.clear();
            self.block_pos = 0;
            return Ok(());
        }
        let extra_len = u16::from_le_bytes([header[10], header[11]]) as usize;
        if HEADER_LEN + extra_len > header.len() {
            return Err(bgzf_error("extra field is implausibly large"));
        }
        self.inner
            .read_exact(&mut header[HEADER_LEN..HEADER_LEN + extra_len])?;
        let block = parse_block_header(&header[..HEADER_LEN + extra_len])?;

        let payload_len = block.compressed_len - block.payload_offset - TRAILER_LEN;
        self.raw.resize(payload_len, 0);
        self.inner.read_exact(&mut self.raw)?;
        let mut trailer = [0u8; TRAILER_LEN];
        self.inner.read_exact(&mut trailer)?;
        let expected_crc = u32::from_le_bytes(trailer[..4].try_into().expect("4 bytes"));
        let expected_len = u32::from_le_bytes(trailer[4..].try_into().expect("4 bytes")) as usize;

        let mut decompressed = std::mem::take(&mut self.block);
        let result = inflate(&self.raw, expected_len, &mut decompressed);
        self.block = decompressed;
        result?;

        let mut crc = flate2::Crc::new();
        crc.update(&self.block);
        if crc.sum() != expected_crc {
            return Err(bgzf_error("block CRC32 does not match"));
        }

        self.block_pos = 0;
        self.next_compressed += block.compressed_len as u64;
        // The EOF marker decompresses to nothing; treat it as end of stream.
        if self.block.is_empty() {
            self.eof = true;
        }
        Ok(())
    }

    /// Make sure the current block has unread bytes, or set `eof`.
    fn fill(&mut self) -> Result<()> {
        while !self.eof && self.block_pos == self.block.len() {
            let consumed = self.block.len() as u64;
            self.block_start += consumed;
            self.read_next_block()?;
        }
        Ok(())
    }

    /// The index, built on demand if it was not supplied.
    fn ensure_index(&mut self) -> Result<()> {
        if self.index.is_none() {
            let saved = self.inner.stream_position()?;
            let index = GziIndex::build(&mut self.inner)?;
            self.inner.seek(SeekFrom::Start(saved))?;
            self.index = Some(index);
        }
        Ok(())
    }
}

impl<R: Read + Seek> Read for BgzfReader<R> {
    fn read(&mut self, buf: &mut [u8]) -> io::Result<usize> {
        if buf.is_empty() {
            return Ok(0);
        }
        self.fill()?;
        if self.eof && self.block_pos == self.block.len() {
            return Ok(0);
        }
        let available = &self.block[self.block_pos..];
        let take = available.len().min(buf.len());
        buf[..take].copy_from_slice(&available[..take]);
        self.block_pos += take;
        Ok(take)
    }
}

impl<R: Read + Seek> Seek for BgzfReader<R> {
    /// Seek in *uncompressed* coordinates.
    ///
    /// `SeekFrom::End` needs the total uncompressed size, so it requires an
    /// index that covers the whole file.
    fn seek(&mut self, pos: SeekFrom) -> io::Result<u64> {
        self.ensure_index()?;

        let target = match pos {
            SeekFrom::Start(offset) => offset,
            SeekFrom::Current(delta) => add_signed(self.position(), delta)?,
            SeekFrom::End(delta) => {
                let end = self
                    .index
                    .as_ref()
                    .and_then(|index| index.uncompressed_len())
                    .ok_or_else(|| {
                        io::Error::new(
                            io::ErrorKind::InvalidInput,
                            "BGZF: this index records block starts only, so the end of the \
                             stream is unknown; rebuild it with GziIndex::build",
                        )
                    })?;
                add_signed(end, delta)?
            }
        };

        // Staying inside the current block is the common case for short hops.
        let within = target.checked_sub(self.block_start);
        if let Some(within) = within {
            if !self.block.is_empty() && within <= self.block.len() as u64 {
                self.block_pos = within as usize;
                return Ok(target);
            }
        }

        // `BlockOffset` is `Copy`, so the borrow of the index ends here and the
        // load below needs no clone of it.
        let block = self
            .index
            .as_ref()
            .and_then(|index| index.block_for(target))
            .ok_or_else(|| {
                io::Error::new(
                    io::ErrorKind::InvalidInput,
                    "BGZF: no block covers that offset",
                )
            })?;
        self.load_block_at(block.compressed, block.uncompressed)?;
        let within = (target - block.uncompressed) as usize;
        if within > self.block.len() {
            // Past the end of the data: leave the cursor at the end.
            self.block_pos = self.block.len();
            self.eof = true;
        } else {
            self.block_pos = within;
        }
        Ok(target)
    }
}

fn add_signed(base: u64, delta: i64) -> io::Result<u64> {
    let result = if delta >= 0 {
        base.checked_add(delta as u64)
    } else {
        base.checked_sub(delta.unsigned_abs())
    };
    result.ok_or_else(|| {
        io::Error::new(
            io::ErrorKind::InvalidInput,
            "BGZF: seek would leave the file",
        )
    })
}

/// A writer that emits BGZF blocks.
///
/// Output is valid gzip, so any gzip tool can read it, and it is seekable by
/// anything that understands BGZF. Call [`BgzfWriter::finish`] to append the
/// EOF marker; without it the file looks truncated to samtools.
pub struct BgzfWriter<W: Write> {
    /// `None` once `finish` has handed the writer back. Wrapped in an `Option`
    /// only because a type with a `Drop` impl cannot give a field away.
    inner: Option<W>,
    buffer: Vec<u8>,
    level: CompressionLevel,
    /// Block starts, recorded so that an index can be written afterwards.
    blocks: Vec<BlockOffset>,
    compressed: u64,
    uncompressed: u64,
    finished: bool,
}

impl BgzfWriter<BufWriter<File>> {
    /// Create a BGZF file.
    pub fn create<P: AsRef<Path>>(path: P) -> Result<BgzfWriter<BufWriter<File>>> {
        let path = path.as_ref();
        let file = File::create(path)
            .map_err(|e| Error::Io(io::Error::new(e.kind(), format!("{}: {e}", path.display()))))?;
        Ok(BgzfWriter::new(BufWriter::with_capacity(128 * 1024, file)))
    }
}

impl<W: Write> BgzfWriter<W> {
    /// Wrap a writer, compressing at the default level.
    pub fn new(inner: W) -> BgzfWriter<W> {
        BgzfWriter::with_level(inner, CompressionLevel::default())
    }

    /// Wrap a writer, compressing at `level`.
    pub fn with_level(inner: W, level: CompressionLevel) -> BgzfWriter<W> {
        BgzfWriter {
            inner: Some(inner),
            buffer: Vec::with_capacity(MAX_BLOCK_PAYLOAD),
            level,
            blocks: vec![BlockOffset {
                compressed: 0,
                uncompressed: 0,
            }],
            compressed: 0,
            uncompressed: 0,
            finished: false,
        }
    }

    /// The index describing the blocks written *so far*.
    ///
    /// Buffered data that has not been compressed yet is not in it, so calling
    /// this before [`BgzfWriter::finish_with_index`] gives an index that stops
    /// short of the end of the file. Prefer `finish_with_index`, which cannot be
    /// wrong.
    pub fn index(&self) -> GziIndex {
        GziIndex {
            blocks: self.blocks.clone(),
            total_uncompressed: Some(self.uncompressed),
        }
    }

    /// The underlying writer, or an error once `finish` has taken it.
    fn sink(&mut self) -> io::Result<&mut W> {
        self.inner.as_mut().ok_or_else(|| {
            io::Error::new(
                io::ErrorKind::BrokenPipe,
                "BGZF: writer was already finished",
            )
        })
    }

    /// Compress and emit whatever is buffered as one block.
    fn flush_block(&mut self) -> io::Result<()> {
        if self.buffer.is_empty() {
            return Ok(());
        }
        let payload = deflate_raw(&self.buffer, self.level)?;
        let block_len = HEADER_LEN + EXTRA_LEN + payload.len() + TRAILER_LEN;
        if block_len > u16::MAX as usize + 1 {
            // Cannot happen with MAX_BLOCK_PAYLOAD, but a wrong constant here
            // would produce files other tools silently misread.
            return Err(io::Error::new(
                io::ErrorKind::InvalidData,
                "BGZF: block would exceed 64 KiB",
            ));
        }

        let mut header = [0u8; HEADER_LEN + EXTRA_LEN];
        header[0] = 0x1f;
        header[1] = 0x8b;
        header[2] = 8; // deflate
        header[3] = 4; // FEXTRA
                       // MTIME stays zero: reproducible output matters more than a timestamp.
        header[9] = 0xff; // unknown OS
        header[10..12].copy_from_slice(&(EXTRA_LEN as u16).to_le_bytes());
        header[12] = b'B';
        header[13] = b'C';
        header[14..16].copy_from_slice(&2u16.to_le_bytes());
        header[16..18].copy_from_slice(&((block_len - 1) as u16).to_le_bytes());

        let mut crc = flate2::Crc::new();
        crc.update(&self.buffer);
        let checksum = crc.sum().to_le_bytes();
        let payload_len = (self.buffer.len() as u32).to_le_bytes();

        let sink = self.sink()?;
        sink.write_all(&header)?;
        sink.write_all(&payload)?;
        sink.write_all(&checksum)?;
        sink.write_all(&payload_len)?;

        self.compressed += block_len as u64;
        self.uncompressed += self.buffer.len() as u64;
        self.blocks.push(BlockOffset {
            compressed: self.compressed,
            uncompressed: self.uncompressed,
        });
        self.buffer.clear();
        Ok(())
    }

    /// Flush the pending block, append the EOF marker and hand the writer back.
    ///
    /// Dropping the writer does the same on a best-effort basis, but only
    /// `finish` reports a failure.
    pub fn finish(self) -> Result<W> {
        self.finish_with_index().map(|(inner, _)| inner)
    }

    /// Finish, and return the complete index alongside the writer.
    ///
    /// This is the safe way to obtain a `.gzi`: every block has been emitted by
    /// the time the index is taken, so it cannot be missing the tail.
    ///
    /// ```no_run
    /// use fastx::bgzf::BgzfWriter;
    /// use std::io::Write;
    ///
    /// let mut writer = BgzfWriter::create("reads.fq.gz")?;
    /// writer.write_all(b"@r\nACGT\n+\nIIII\n")?;
    /// let (_file, index) = writer.finish_with_index()?;
    /// index.write_to_path("reads.fq.gz")?;   // reads.fq.gz.gzi
    /// # Ok::<(), fastx::Error>(())
    /// ```
    pub fn finish_with_index(mut self) -> Result<(W, GziIndex)> {
        self.finish_in_place()?;
        let index = self.index();
        let inner = self
            .inner
            .take()
            .ok_or_else(|| Error::Other("BGZF: writer was already finished".to_string()))?;
        Ok((inner, index))
    }

    fn finish_in_place(&mut self) -> Result<()> {
        if self.finished || self.inner.is_none() {
            return Ok(());
        }
        self.flush_block()?;
        let sink = self.sink()?;
        sink.write_all(&EOF_BLOCK)?;
        sink.flush()?;
        self.compressed += EOF_BLOCK.len() as u64;
        self.finished = true;
        Ok(())
    }
}

impl<W: Write> Write for BgzfWriter<W> {
    fn write(&mut self, buf: &[u8]) -> io::Result<usize> {
        let room = MAX_BLOCK_PAYLOAD - self.buffer.len();
        let take = room.min(buf.len());
        self.buffer.extend_from_slice(&buf[..take]);
        if self.buffer.len() == MAX_BLOCK_PAYLOAD {
            self.flush_block()?;
        }
        Ok(take)
    }

    /// Ends the current block, so the next byte starts a fresh one.
    ///
    /// This is what makes a flush point seekable, and it costs a little
    /// compression, so flush at record boundaries rather than per record.
    fn flush(&mut self) -> io::Result<()> {
        self.flush_block()?;
        self.sink()?.flush()
    }
}

impl<W: Write> Drop for BgzfWriter<W> {
    fn drop(&mut self) {
        // Best effort: a caller who wants to see errors uses finish().
        let _ = self.finish_in_place();
    }
}

/// Deflate with no zlib or gzip wrapper, which is what a gzip member holds.
fn deflate_raw(data: &[u8], level: CompressionLevel) -> io::Result<Vec<u8>> {
    use flate2::write::DeflateEncoder;
    let mut encoder = DeflateEncoder::new(
        Vec::with_capacity(data.len() / 2 + 64),
        flate2::Compression::new(level.0.min(9)),
    );
    encoder.write_all(data)?;
    encoder.finish()
}

/// Read into `buf` fully, or report `false` if the reader was already at EOF.
fn read_exact_or_eof<R: Read>(reader: &mut R, buf: &mut [u8]) -> Result<bool> {
    let mut filled = 0;
    while filled < buf.len() {
        match reader.read(&mut buf[filled..]) {
            Ok(0) if filled == 0 => return Ok(false),
            Ok(0) => return Err(bgzf_error("file ends in the middle of a block")),
            Ok(n) => filled += n,
            Err(e) if e.kind() == io::ErrorKind::Interrupted => continue,
            Err(e) => return Err(Error::Io(e)),
        }
    }
    Ok(true)
}

/// Read as much as is available, up to `buf.len()`.
fn read_up_to<R: Read>(reader: &mut R, buf: &mut [u8]) -> Result<usize> {
    let mut filled = 0;
    while filled < buf.len() {
        match reader.read(&mut buf[filled..]) {
            Ok(0) => break,
            Ok(n) => filled += n,
            Err(e) if e.kind() == io::ErrorKind::Interrupted => continue,
            Err(e) => return Err(Error::Io(e)),
        }
    }
    Ok(filled)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::io::Cursor;

    fn compress(data: &[u8]) -> Vec<u8> {
        let mut writer = BgzfWriter::new(Vec::new());
        writer.write_all(data).unwrap();
        writer.finish().unwrap()
    }

    #[test]
    fn round_trips_through_our_own_reader() {
        for size in [
            0usize,
            1,
            100,
            MAX_BLOCK_PAYLOAD - 1,
            MAX_BLOCK_PAYLOAD,
            MAX_BLOCK_PAYLOAD + 1,
            300_000,
        ] {
            let data: Vec<u8> = (0..size).map(|i| b"ACGTN"[i % 5]).collect();
            let compressed = compress(&data);
            assert!(is_bgzf(&compressed), "size {size} did not produce BGZF");

            let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();
            let mut out = Vec::new();
            reader.read_to_end(&mut out).unwrap();
            assert_eq!(out, data, "size {size}");
        }
    }

    #[test]
    fn output_is_plain_gzip_too() {
        // The whole point of BGZF: ordinary gzip tools must still read it.
        let data: Vec<u8> = (0..200_000).map(|i| b"ACGT"[i % 4]).collect();
        let compressed = compress(&data);
        let mut out = Vec::new();
        flate2::read::MultiGzDecoder::new(&compressed[..])
            .read_to_end(&mut out)
            .unwrap();
        assert_eq!(out, data);
    }

    #[test]
    fn ends_with_the_eof_marker() {
        let compressed = compress(b"ACGT");
        assert_eq!(
            &compressed[compressed.len() - EOF_BLOCK.len()..],
            &EOF_BLOCK
        );
        // An empty file is just the marker.
        assert_eq!(compress(b""), EOF_BLOCK.to_vec());
    }

    #[test]
    fn blocks_stay_within_the_size_limit() {
        // Incompressible data is the case that could overflow a block.
        let data: Vec<u8> = (0..500_000)
            .map(|i| ((i * 2_654_435_761u64 as usize) >> 7) as u8)
            .collect();
        let compressed = compress(&data);
        let index = GziIndex::build(Cursor::new(&compressed)).unwrap();
        for pair in index.blocks().windows(2) {
            let block_len = pair[1].compressed - pair[0].compressed;
            assert!(block_len <= 65_536, "block of {block_len} bytes");
        }
        let mut out = Vec::new();
        BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .read_to_end(&mut out)
            .unwrap();
        assert_eq!(out, data);
    }

    #[test]
    fn index_from_the_writer_matches_a_rescan() {
        let data: Vec<u8> = (0..250_000).map(|i| b"ACGTN"[i % 5]).collect();
        let mut writer = BgzfWriter::new(Vec::new());
        writer.write_all(&data).unwrap();
        let from_writer = writer.index();
        let compressed = writer.finish().unwrap();

        let from_scan = GziIndex::build(Cursor::new(&compressed)).unwrap();
        // The rescan also sees the EOF block, which the writer's index predates.
        assert_eq!(
            &from_scan.blocks()[..from_writer.len()],
            from_writer.blocks()
        );
        assert_eq!(from_scan.uncompressed_len(), Some(data.len() as u64));
    }

    #[test]
    fn gzi_round_trips_and_omits_the_first_block() {
        let data: Vec<u8> = (0..200_000).map(|i| b"ACGT"[i % 4]).collect();
        let compressed = compress(&data);
        let index = GziIndex::build(Cursor::new(&compressed)).unwrap();

        let mut text = Vec::new();
        index.write(&mut text).unwrap();
        // 8-byte count plus 16 bytes per entry, first block implicit.
        assert_eq!(text.len(), 8 + 16 * (index.len() - 1));
        assert_eq!(
            u64::from_le_bytes(text[..8].try_into().unwrap()),
            index.len() as u64 - 1
        );

        let reparsed = GziIndex::parse(&text[..]).unwrap();
        assert_eq!(reparsed.blocks(), index.blocks());
        // A parsed index cannot know where the data ends, and must not pretend.
        assert_eq!(reparsed.uncompressed_len(), None);
        assert_eq!(index.uncompressed_len(), Some(200_000));
    }

    #[test]
    fn seek_from_end_needs_a_scanned_index() {
        let data: Vec<u8> = (0..200_000).map(|i| b"ACGT"[i % 4]).collect();
        let compressed = compress(&data);
        let scanned = GziIndex::build(Cursor::new(&compressed)).unwrap();
        let mut text = Vec::new();
        scanned.write(&mut text).unwrap();
        let parsed = GziIndex::parse(&text[..]).unwrap();

        // With block starts only, End-relative seeks must fail loudly rather
        // than landing at the start of the last block.
        let mut reader = BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .with_index(parsed);
        assert!(reader.seek(SeekFrom::End(0)).is_err());
        // Absolute seeks still work.
        reader.seek(SeekFrom::Start(199_998)).unwrap();
        let mut tail = Vec::new();
        reader.read_to_end(&mut tail).unwrap();
        assert_eq!(tail, &data[199_998..]);

        let mut reader = BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .with_index(scanned);
        assert_eq!(reader.seek(SeekFrom::End(0)).unwrap(), 200_000);
    }

    #[test]
    fn seeks_to_any_offset() {
        let data: Vec<u8> = (0..300_000).map(|i| (i % 251) as u8).collect();
        let compressed = compress(&data);
        let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();

        // Offsets that land in the first block, deep inside, and on boundaries.
        for target in [
            0usize,
            1,
            MAX_BLOCK_PAYLOAD - 1,
            MAX_BLOCK_PAYLOAD,
            MAX_BLOCK_PAYLOAD + 1,
            2 * MAX_BLOCK_PAYLOAD,
            299_999,
        ] {
            reader.seek(SeekFrom::Start(target as u64)).unwrap();
            assert_eq!(reader.position(), target as u64);
            let mut buf = [0u8; 8];
            let want = (data.len() - target).min(buf.len());
            reader.read_exact(&mut buf[..want]).unwrap();
            assert_eq!(&buf[..want], &data[target..target + want], "at {target}");
        }

        // Relative and end-relative seeks.
        reader.seek(SeekFrom::Start(10)).unwrap();
        reader.seek(SeekFrom::Current(5)).unwrap();
        assert_eq!(reader.position(), 15);
        assert_eq!(reader.seek(SeekFrom::End(0)).unwrap(), data.len() as u64);
        let mut rest = Vec::new();
        reader.read_to_end(&mut rest).unwrap();
        assert!(rest.is_empty());

        // Seeking backwards must work as well as forwards.
        reader.seek(SeekFrom::Start(7)).unwrap();
        let mut buf = [0u8; 4];
        reader.read_exact(&mut buf).unwrap();
        assert_eq!(&buf, &data[7..11]);
    }

    #[test]
    fn seek_before_the_start_is_an_error() {
        let compressed = compress(b"ACGT");
        let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();
        assert!(reader.seek(SeekFrom::Current(-1)).is_err());
        assert!(reader.seek(SeekFrom::End(-100)).is_err());
    }

    #[test]
    fn rejects_plain_gzip_and_garbage() {
        // Valid gzip, but no BC extra field: not BGZF.
        let mut plain = Vec::new();
        {
            let mut encoder =
                flate2::write::GzEncoder::new(&mut plain, flate2::Compression::default());
            encoder.write_all(b"ACGT").unwrap();
            encoder.finish().unwrap();
        }
        assert!(!is_bgzf(&plain));
        assert!(BgzfReader::new(Cursor::new(&plain)).is_err());

        assert!(BgzfReader::new(Cursor::new(b"not gzip at all".to_vec())).is_err());
        assert!(BgzfReader::new(Cursor::new(Vec::new())).is_err());
    }

    #[test]
    fn detects_a_corrupt_block() {
        let mut compressed = compress(&vec![b'A'; 5_000]);
        // Flip a byte in the deflate payload; the CRC or the inflate must catch it.
        let victim = HEADER_LEN + EXTRA_LEN + 5;
        compressed[victim] ^= 0xff;
        let mut out = Vec::new();
        let result = BgzfReader::new(Cursor::new(&compressed))
            .unwrap()
            .read_to_end(&mut out);
        assert!(result.is_err(), "corruption went unnoticed");
    }

    #[test]
    fn truncated_file_is_an_error_not_silent_truncation() {
        let compressed = compress(&vec![b'A'; 200_000]);
        let cut = compressed.len() / 2;
        let mut out = Vec::new();
        let result = BgzfReader::new(Cursor::new(compressed[..cut].to_vec()))
            .unwrap()
            .read_to_end(&mut out);
        assert!(result.is_err(), "truncation went unnoticed");
    }

    #[test]
    fn parse_rejects_an_implausible_index() {
        let mut bad = u64::MAX.to_le_bytes().to_vec();
        bad.extend_from_slice(&[0u8; 16]);
        assert!(GziIndex::parse(&bad[..]).is_err());
        assert!(GziIndex::parse(&[0u8; 3][..]).is_err());
    }

    #[test]
    fn flush_starts_a_new_block() {
        let mut writer = BgzfWriter::new(Vec::new());
        writer.write_all(b"first").unwrap();
        writer.flush().unwrap();
        writer.write_all(b"second").unwrap();
        // index() before finishing sees only the flushed block, which is exactly
        // why finish_with_index exists.
        assert_eq!(writer.index().len(), 2);
        let (compressed, index) = writer.finish_with_index().unwrap();

        // Two data blocks, each a seek point.
        assert_eq!(index.len(), 3); // implicit start + two blocks
        assert_eq!(index.blocks()[1].uncompressed, 5);
        assert_eq!(index.uncompressed_len(), Some(11));

        let mut reader = BgzfReader::new(Cursor::new(&compressed)).unwrap();
        reader.seek(SeekFrom::Start(5)).unwrap();
        let mut out = Vec::new();
        reader.read_to_end(&mut out).unwrap();
        assert_eq!(out, b"second");
    }
}