gwseq_io/bam/

record.rs

//! One alignment.
//!
//! The design point is laziness: `cigar`, `sequence`,
//! `qualities` and `tags` are decoded on first read and kept, so an alignment
//! read for its coordinates never pays for the rest of it.
//!
//! The record owns `Bytes` — a refcounted slice of the decompressed BGZF chunk
//! — so there is no per-record copy, and each derived field is a `OnceLock`.
//!
//! `OnceLock`, not the `OnceCell` a single-threaded cache would want: PyO3
//! requires `#[pyclass]` types to be `Send + Sync`, and `OnceCell` is not
//! `Sync`. The alternative — `#[pyclass(unsendable)]` — would pin every
//! alignment to the thread that read it, which breaks handing a list of them
//! to a Python thread pool. The cost is one acquire load per access on a field read a handful of times per
//! record, against a decode that walks the whole record.

use std::sync::{Arc, OnceLock};

use bytes::Bytes;

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

/// Bytes of a record before its read name: `block_size` through `tlen`.
pub const RECORD_HEADER_SIZE: usize = 36;

/// The 36-byte fixed part of a record, decoded eagerly.
#[derive(Debug, Clone, Copy)]
pub struct RecordCore {
    pub ref_id: i32,
    pub pos: i32,
    pub mapq: u8,
    pub bai_bin: u16,
    pub n_cigar_op: u16,
    pub flag: u16,
    pub l_seq: i32,
    pub next_ref_id: i32,
    pub next_pos: i32,
    pub tlen: i32,
}

/// Where each variable-length field sits inside the record's own bytes.
#[derive(Debug, Clone, Copy)]
struct Layout {
    name: usize,
    cigar: usize,
    seq: usize,
    qual: usize,
    tags: usize,
    end: usize,
}

/// An optional field's value, typed as the file types it.
#[derive(Debug, Clone, PartialEq)]
pub enum TagValue {
    Char(char),
    Int(i64),
    Float(f32),
    Str(String),
    IntArray(Vec<i64>),
    FloatArray(Vec<f32>),
}

pub struct BamRecord {
    /// The record's own bytes, from `block_size` to its end.
    raw: Bytes,
    core: RecordCore,
    layout: Layout,
    /// Reference length from the cigar, resolved at decode: the filter needs it
    /// before anything else is built.
    end: i64,
    cigar_ops: Vec<u32>,
    /// Whether the optional fields were kept at all.
    has_tags: bool,
    cigar: OnceLock<String>,
    sequence: OnceLock<String>,
    qualities: OnceLock<String>,
    /// Holds a `Result`: a record with malformed optional fields reads fine and
    /// fails here, which is the documented behaviour.
    tags: OnceLock<Result<Vec<(String, TagValue)>>>,
    /// Resolved reference names, shared with the reader.
    chr_names: Arc<Vec<String>>,
}

impl std::fmt::Debug for BamRecord {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("BamRecord")
            .field("chr", &self.chr())
            .field("start", &self.start())
            .field("end", &self.end)
            .field("flag", &self.core.flag)
            .finish()
    }
}

impl BamRecord {
    fn name_of(&self, index: i32) -> &str {
        if index < 0 {
            return "*";
        }
        self.chr_names
            .get(index as usize)
            .map(String::as_str)
            .unwrap_or("*")
    }

    pub fn chr(&self) -> &str {
        self.name_of(self.core.ref_id)
    }

    pub fn start(&self) -> i64 {
        self.core.pos as i64
    }

    /// 0-based half-open end, derived from the cigar. Equal to `start` for an
    /// alignment covering no reference.
    pub fn end(&self) -> i64 {
        self.end
    }

    pub fn read_name(&self) -> &str {
        let bytes = &self.raw[self.layout.name..self.layout.cigar];
        // NUL-terminated in the file.
        let bytes = bytes.strip_suffix(&[0]).unwrap_or(bytes);
        // A name is ASCII in every file anyone writes. One that is not gives
        // back the part of it that is, rather than the empty string, which is
        // what every such name used to come back as — so two of them collided
        // with each other and with a record that has no name at all. Truncating
        // costs no allocation, which a lossy conversion on a type that exists
        // in millions would.
        match std::str::from_utf8(bytes) {
            Ok(name) => name,
            Err(e) => std::str::from_utf8(&bytes[..e.valid_up_to()]).unwrap_or(""),
        }
    }

    pub fn flag(&self) -> u16 {
        self.core.flag
    }

    pub fn mapping_quality(&self) -> u8 {
        self.core.mapq
    }

    pub fn bai_bin(&self) -> u16 {
        self.core.bai_bin
    }

    pub fn next_chr(&self) -> &str {
        self.name_of(self.core.next_ref_id)
    }

    pub fn next_start(&self) -> i64 {
        self.core.next_pos as i64
    }

    pub fn template_length(&self) -> i64 {
        self.core.tlen as i64
    }

    /// The CIGAR, e.g. `"10S80M10S"`.
    ///
    /// A cigar of more than 65 535 operations does not fit the record's own
    /// field, so it is stored in a `CG` optional field and read from there at
    /// decode; the placeholder the record carries in its place is never
    /// returned, and the `CG` field is left in [`tags`](Self::tags).
    pub fn cigar(&self) -> &str {
        self.cigar.get_or_init(|| decode_cigar(&self.cigar_ops))
    }

    /// Bases of the reference the alignment covers.
    pub fn reference_length(&self) -> i64 {
        reference_length(&self.cigar_ops)
    }

    /// Bases of the read its cigar consumes.
    pub fn query_length(&self) -> i64 {
        query_length(&self.cigar_ops)
    }

    /// SEQ, unpacked from its 4-bit encoding, or `"*"` when the record carries
    /// none.
    pub fn sequence(&self) -> &str {
        self.sequence.get_or_init(|| {
            decode_sequence(
                self.core.l_seq as i64,
                &self.raw[self.layout.seq..self.layout.qual],
            )
        })
    }

    /// QUAL as phred+33, or `"*"` when the record carries no qualities at all.
    /// A record with only some of them missing spells those `"*"` in place.
    pub fn qualities(&self) -> &str {
        self.qualities
            .get_or_init(|| decode_qualities(&self.raw[self.layout.qual..self.layout.tags]))
    }

    /// Optional fields by two-letter tag, in the order the record stores them.
    ///
    /// Empty when the record was read with tags off. A record whose optional
    /// fields are malformed reads fine and fails **here**.
    pub fn tags(&self) -> Result<&[(String, TagValue)]> {
        if !self.has_tags {
            return Ok(&[]);
        }
        let parsed = self
            .tags
            .get_or_init(|| parse_tags(&self.raw[self.layout.tags..self.layout.end]));
        match parsed {
            Ok(tags) => Ok(tags),
            // The cell holds the failure so a second read repeats it rather
            // than re-walking bytes already known to be bad.
            //
            // Named by the alignment rather than by the file: a record does not
            // carry its path — there are millions of them and it would be eight
            // bytes each — and a caller holding one already knows which file it
            // came out of. Which *alignment* is what they do not know, and it is
            // what a `Corrupt` here used to say nothing about: the message read
            // ": unterminated bam string tag RG (at offset 0)".
            Err(e) => Err(Error::corrupt(
                format!("{}:{}-{}", self.chr(), self.start(), self.end()),
                self.core.pos.max(0) as u64,
                format!("alignment {}: {e}", self.read_name()),
            )),
        }
    }

    // The flag bits the API decodes. Finer ones are the caller's to mask off
    // `flag()`.
    pub fn is_paired(&self) -> bool {
        self.core.flag & 0x001 != 0
    }
    pub fn is_proper_pair(&self) -> bool {
        self.core.flag & 0x002 != 0
    }
    pub fn is_mapped(&self) -> bool {
        self.core.flag & 0x004 == 0
    }
    pub fn is_next_mapped(&self) -> bool {
        self.core.flag & 0x008 == 0
    }
    pub fn is_reverse(&self) -> bool {
        self.core.flag & 0x010 != 0
    }
    pub fn is_next_reverse(&self) -> bool {
        self.core.flag & 0x020 != 0
    }
    pub fn is_first_in_pair(&self) -> bool {
        self.core.flag & 0x040 != 0
    }
    pub fn is_last_in_pair(&self) -> bool {
        self.core.flag & 0x080 != 0
    }
    pub fn is_secondary_or_supplementary(&self) -> bool {
        self.core.flag & (0x100 | 0x800) != 0
    }
    pub fn is_failed_qc_or_duplicate(&self) -> bool {
        self.core.flag & (0x200 | 0x400) != 0
    }
}

/// M, D, N, `=`, X consume reference.
fn reference_length(ops: &[u32]) -> i64 {
    ops.iter()
        .filter(|op| matches!(*op & 0xF, 0 | 2 | 3 | 7 | 8))
        .map(|op| (op >> 4) as i64)
        .sum()
}

/// M, I, S, `=`, X consume query.
fn query_length(ops: &[u32]) -> i64 {
    ops.iter()
        .filter(|op| matches!(*op & 0xF, 0 | 1 | 4 | 7 | 8))
        .map(|op| (op >> 4) as i64)
        .sum()
}

fn decode_cigar(ops: &[u32]) -> String {
    const OPS: &[u8] = b"MIDNSHP=X";
    let mut out = String::with_capacity(ops.len() * 4);
    for op in ops {
        use std::fmt::Write as _;
        let _ = write!(out, "{}", op >> 4);
        out.push(match OPS.get((op & 0xF) as usize) {
            Some(c) => *c as char,
            None => '?',
        });
    }
    out
}

fn decode_sequence(l_seq: i64, packed: &[u8]) -> String {
    // SAM writes an absent sequence as "*". An empty string would read back as
    // a sequence of no bases rather than as the absence of one.
    if l_seq <= 0 {
        return "*".to_string();
    }
    const LOOKUP: &[u8] = b"=ACMGRSVTWYHKDBN";
    // Bounded by the read length rather than by the bytes holding it: an odd
    // length leaves the low half of the last byte unused, and counting bytes
    // would decode that padding as a base of its own.
    let len = l_seq.min(packed.len() as i64 * 2) as usize;
    let mut out = String::with_capacity(len);
    for i in 0..len {
        let byte = packed[i / 2];
        let code = if i % 2 == 0 { byte >> 4 } else { byte & 0xF };
        out.push(LOOKUP[code as usize] as char);
    }
    out
}

fn decode_qualities(quals: &[u8]) -> String {
    // A record with no qualities stores 0xFF for every base, and SAM writes
    // that back as a single "*". Rendering it per base would give a string as
    // long as the read and indistinguishable from real scores of phred 9,
    // '*' being 42 = 33 + 9.
    if quals.iter().all(|q| *q == 0xFF) {
        return "*".to_string();
    }
    // Widened before the offset: `q + 33` in `u8` overflows for any quality
    // above 222, which is a panic in a debug build and a wrapped character in a
    // release one. Phred qualities do not go that high, but these bytes come out
    // of a file and nothing in the format stops them.
    quals
        .iter()
        .map(|q| {
            if *q == 0xFF {
                '*'
            } else {
                char::from_u32(*q as u32 + 33).unwrap_or('?')
            }
        })
        .collect()
}

/// Walk the optional fields, in the order the record stores them.
///
/// Every length read out of the block is checked against its end before
/// anything is read with it: these are bytes from a file, and a corrupt block
/// supplies them.
fn parse_tags(raw: &[u8]) -> Result<Vec<(String, TagValue)>> {
    let mut out = Vec::new();
    let mut at = 0usize;
    let size = raw.len();
    // The message alone, without a path or an offset: neither is known here —
    // these bytes are a slice of a record, which is a slice of a decompressed
    // block — and `BamRecord::tags` wraps this in one that names the alignment.
    let bad = |what: String| Error::invalid(what);

    while at < size {
        if at + 3 > size {
            return Err(bad(
                "truncated bam tag (no room for its tag and type)".into()
            ));
        }
        // Exactly the two bytes the tag is, never trimmed: a corrupt block could
        // hold a NUL there, and a shorter tag would then be reported.
        let tag = String::from_utf8_lossy(&raw[at..at + 2]).into_owned();
        let kind = raw[at + 2];
        at += 3;

        let value = match kind {
            b'B' => {
                if at + 5 > size {
                    return Err(bad(format!("truncated bam array tag {tag}")));
                }
                let subtype = raw[at];
                let count = u32::from_le_bytes([raw[at + 1], raw[at + 2], raw[at + 3], raw[at + 4]])
                    as usize;
                at += 5;
                let element = match subtype {
                    b'c' | b'C' => 1usize,
                    b's' | b'S' => 2,
                    b'i' | b'I' | b'f' => 4,
                    other => {
                        return Err(bad(format!(
                            "unsupported bam array tag subtype {}",
                            other as char
                        )))
                    }
                };
                let bytes = count
                    .checked_mul(element)
                    .filter(|n| at + n <= size)
                    .ok_or_else(|| bad(format!("bam array tag {tag} runs past its record")))?;
                let data = &raw[at..at + bytes];
                at += bytes;
                read_array(subtype, data)
            }
            b'Z' | b'H' => {
                // Bounded by the block: an unterminated string would otherwise
                // be scanned to the end of it and past it.
                let start = at;
                while at < size && raw[at] != 0 {
                    at += 1;
                }
                if at >= size {
                    return Err(bad(format!("unterminated bam string tag {tag}")));
                }
                let value = String::from_utf8_lossy(&raw[start..at]).into_owned();
                at += 1;
                TagValue::Str(value)
            }
            other => {
                let width = match other {
                    b'A' | b'c' | b'C' => 1usize,
                    b's' | b'S' => 2,
                    b'i' | b'I' | b'f' => 4,
                    _ => return Err(bad(format!("unsupported tag type {}", other as char))),
                };
                if at + width > size {
                    return Err(bad(format!("bam tag {tag} runs past its record")));
                }
                let data = &raw[at..at + width];
                at += width;
                read_scalar(other, data)
            }
        };
        out.push((tag, value));
    }
    Ok(out)
}

fn read_scalar(kind: u8, data: &[u8]) -> TagValue {
    match kind {
        b'A' => TagValue::Char(data[0] as char),
        b'c' => TagValue::Int(data[0] as i8 as i64),
        b'C' => TagValue::Int(data[0] as i64),
        b's' => TagValue::Int(i16::from_le_bytes([data[0], data[1]]) as i64),
        b'S' => TagValue::Int(u16::from_le_bytes([data[0], data[1]]) as i64),
        b'i' => TagValue::Int(i32::from_le_bytes([data[0], data[1], data[2], data[3]]) as i64),
        b'I' => TagValue::Int(u32::from_le_bytes([data[0], data[1], data[2], data[3]]) as i64),
        _ => TagValue::Float(f32::from_le_bytes([data[0], data[1], data[2], data[3]])),
    }
}

fn read_array(subtype: u8, data: &[u8]) -> TagValue {
    match subtype {
        b'f' => TagValue::FloatArray(
            data.chunks_exact(4)
                .map(|c| f32::from_le_bytes([c[0], c[1], c[2], c[3]]))
                .collect(),
        ),
        b'c' => TagValue::IntArray(data.iter().map(|b| *b as i8 as i64).collect()),
        b'C' => TagValue::IntArray(data.iter().map(|b| *b as i64).collect()),
        b's' => TagValue::IntArray(
            data.chunks_exact(2)
                .map(|c| i16::from_le_bytes([c[0], c[1]]) as i64)
                .collect(),
        ),
        b'S' => TagValue::IntArray(
            data.chunks_exact(2)
                .map(|c| u16::from_le_bytes([c[0], c[1]]) as i64)
                .collect(),
        ),
        b'i' => TagValue::IntArray(
            data.chunks_exact(4)
                .map(|c| i32::from_le_bytes([c[0], c[1], c[2], c[3]]) as i64)
                .collect(),
        ),
        _ => TagValue::IntArray(
            data.chunks_exact(4)
                .map(|c| u32::from_le_bytes([c[0], c[1], c[2], c[3]]) as i64)
                .collect(),
        ),
    }
}

/// Whether a record's cigar is the placeholder standing in for a long one.
///
/// A cigar of more than 65 535 operations does not fit `n_cigar_op`, so the
/// record stores `<l_seq>S<ref_len>N` there and the real operations in a `CG`
/// optional field. Coordinates are unaffected either way, the placeholder's `N`
/// carrying the reference length, so this is only about what the cigar itself
/// reads back as.
fn is_long_cigar_placeholder(ops: &[u32], l_seq: i32) -> bool {
    ops.len() == 2
        && (ops[0] & 0xF) == 4 // S
        && (ops[0] >> 4) as i64 == l_seq as i64
        && (ops[1] & 0xF) == 3 // N
}

/// The operations a record's `CG:B,I` field holds, or none when it carries no
/// usable one.
///
/// Malformed optional fields are not an error here: a record reads fine and
/// fails when its tags are looked at, and a long cigar is not the place to break
/// that — the placeholder is kept instead, which is what the file says.
///
/// The `CG` field is left where it is rather than dropped once read, unlike
/// htslib, which strips it. The tags of an alignment are the ones the file
/// stores, and hiding one here would make `tags` disagree with the record.
fn read_long_cigar(tag_bytes: &[u8]) -> Vec<u32> {
    let Ok(tags) = parse_tags(tag_bytes) else {
        return Vec::new();
    };
    for (tag, value) in tags {
        if tag != "CG" {
            continue;
        }
        if let TagValue::IntArray(values) = value {
            return values.into_iter().map(|v| v as u32).collect();
        }
    }
    Vec::new()
}

/// The default `filter=True`: drop unmapped alignments, improper pairs,
/// secondary and supplementary records, and anything marked as failing QC or as
/// a duplicate.
#[derive(Debug, Clone, Copy)]
pub struct RecordFilter {
    pub enabled: bool,
}

impl Default for RecordFilter {
    fn default() -> Self {
        Self { enabled: true }
    }
}

/// Which alignments a read keeps: the region, and optionally the flag rules.
#[derive(Debug, Clone, Copy)]
pub struct EntryFilter {
    /// Reference index, or `None` to accept any.
    pub chr_index: Option<i32>,
    pub start: i64,
    /// `None` to accept any position.
    pub end: Option<i64>,
    pub standard_flags: bool,
}

impl EntryFilter {
    pub fn accepts(&self, chr_index: i32, start: i64, end: i64, flag: u16) -> bool {
        if let Some(wanted) = self.chr_index {
            if chr_index != wanted {
                return false;
            }
        }
        if let Some(region_end) = self.end {
            // An empty region overlaps nothing, wherever the alignment sits.
            if region_end <= self.start {
                return false;
            }
            // An alignment whose cigar consumes no reference — a placed unmapped
            // read, or one cigar'd entirely to insertions — has end == start and
            // would otherwise be reported by a window strictly containing its
            // position but not by one starting on it. htslib's `bam_endpos`
            // treats a reference length of 0 as 1.
            if start >= region_end || end.max(start + 1) <= self.start {
                return false;
            }
        }
        if self.standard_flags {
            if flag & 0x004 != 0 {
                return false; // unmapped
            }
            if flag & 0x001 != 0 && flag & 0x002 == 0 {
                return false; // paired but not properly
            }
            if flag & (0x100 | 0x800) != 0 {
                return false; // secondary or supplementary
            }
            if flag & (0x200 | 0x400) != 0 {
                return false; // failed QC or duplicate
            }
        }
        true
    }
}

/// Read the alignment records packed into a block of decompressed BAM data.
///
/// `block` must start at a record boundary — which is what a virtual offset
/// names. Every length read out of the file is checked against the end of the
/// record it belongs to before anything is read with it.
pub fn decode_block(
    block: &Bytes,
    parse_tags_flag: bool,
    filter: &EntryFilter,
    chr_names: &Arc<Vec<String>>,
    path: &str,
) -> Result<Vec<BamRecord>> {
    let size = block.len();
    let mut at = 0usize;
    let mut out = Vec::new();

    while at < size {
        if at + 4 > size {
            return Err(Error::corrupt(
                path,
                at as u64,
                "truncated bam record (no room for its length)",
            ));
        }
        let block_size =
            u32::from_le_bytes([block[at], block[at + 1], block[at + 2], block[at + 3]]) as usize;
        let record_end = at + 4 + block_size;
        if block_size < RECORD_HEADER_SIZE - 4 || record_end > size {
            return Err(Error::corrupt(
                path,
                at as u64,
                format!(
                    "truncated bam record (declares {block_size} bytes, {} left in the block)",
                    size - at - 4
                ),
            ));
        }

        let r = &block[at..record_end];
        let i32_at = |o: usize| i32::from_le_bytes([r[o], r[o + 1], r[o + 2], r[o + 3]]);
        let u16_at = |o: usize| u16::from_le_bytes([r[o], r[o + 1]]);
        let core = RecordCore {
            ref_id: i32_at(4),
            pos: i32_at(8),
            mapq: r[13],
            bai_bin: u16_at(14),
            n_cigar_op: u16_at(16),
            flag: u16_at(18),
            l_seq: i32_at(20),
            next_ref_id: i32_at(24),
            next_pos: i32_at(28),
            tlen: i32_at(32),
        };
        let l_read_name = r[12] as usize;

        // The variable-length parts follow one another, so checking where the
        // last of them ends checks all of them.
        let name = RECORD_HEADER_SIZE;
        let cigar = name + l_read_name;
        let seq = cigar + core.n_cigar_op as usize * 4;
        let qual = seq + (core.l_seq as usize).div_ceil(2);
        let tags = qual + core.l_seq.max(0) as usize;
        if l_read_name < 1 || core.l_seq < 0 || tags > r.len() {
            return Err(Error::corrupt(
                path,
                at as u64,
                format!(
                    "bam record at {at} declares fields that do not fit its {block_size} bytes"
                ),
            ));
        }
        let layout = Layout {
            name,
            cigar,
            seq,
            qual,
            tags,
            end: r.len(),
        };

        let mut ops: Vec<u32> = r[cigar..seq]
            .chunks_exact(4)
            .map(|c| u32::from_le_bytes([c[0], c[1], c[2], c[3]]))
            .collect();
        let mut end = core.pos as i64 + reference_length(&ops);

        // Tested before the sequence and the tags are built, which is what makes
        // skipping an alignment cheaper than keeping one.
        if !filter.accepts(core.ref_id, core.pos as i64, end, core.flag) {
            at = record_end;
            continue;
        }

        // Read after the filter, and only for the records that carry one: the
        // placeholder gives the same start and end as the cigar it stands for,
        // so an alignment the request does not want costs nothing here. `end` is
        // taken again from the real operations so the two cannot disagree.
        if is_long_cigar_placeholder(&ops, core.l_seq) {
            let real = read_long_cigar(&r[tags..]);
            if !real.is_empty() {
                end = core.pos as i64 + reference_length(&real);
                ops = real;
            }
        }

        out.push(BamRecord {
            raw: block.slice(at..record_end),
            core,
            layout,
            end,
            cigar_ops: ops,
            has_tags: parse_tags_flag,
            cigar: OnceLock::new(),
            sequence: OnceLock::new(),
            qualities: OnceLock::new(),
            tags: OnceLock::new(),
            chr_names: chr_names.clone(),
        });
        at = record_end;
    }
    Ok(out)
}

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

    fn names() -> Arc<Vec<String>> {
        Arc::new(vec!["chr1".to_string(), "chr2".to_string()])
    }

    /// One record, laid out as the format lays one out.
    #[allow(clippy::too_many_arguments)]
    fn record(
        ref_id: i32,
        pos: i32,
        flag: u16,
        name: &str,
        cigar: &[u32],
        seq: &[u8],
        quals: &[u8],
        tags: &[u8],
    ) -> Vec<u8> {
        let l_seq = seq.len() as i32;
        let packed: Vec<u8> = seq
            .chunks(2)
            .map(|pair| {
                let hi = base_code(pair[0]);
                let lo = pair.get(1).map(|b| base_code(*b)).unwrap_or(0);
                (hi << 4) | lo
            })
            .collect();

        let mut body = Vec::new();
        body.extend_from_slice(&ref_id.to_le_bytes());
        body.extend_from_slice(&pos.to_le_bytes());
        body.push(name.len() as u8 + 1);
        body.push(60); // mapq
        body.extend_from_slice(&0u16.to_le_bytes()); // bin
        body.extend_from_slice(&(cigar.len() as u16).to_le_bytes());
        body.extend_from_slice(&flag.to_le_bytes());
        body.extend_from_slice(&l_seq.to_le_bytes());
        body.extend_from_slice(&(-1i32).to_le_bytes()); // next_ref_id
        body.extend_from_slice(&(-1i32).to_le_bytes()); // next_pos
        body.extend_from_slice(&0i32.to_le_bytes()); // tlen
        body.extend_from_slice(name.as_bytes());
        body.push(0);
        for op in cigar {
            body.extend_from_slice(&op.to_le_bytes());
        }
        body.extend_from_slice(&packed);
        body.extend_from_slice(quals);
        body.extend_from_slice(tags);

        let mut out = ((body.len()) as u32).to_le_bytes().to_vec();
        out.extend_from_slice(&body);
        out
    }

    fn base_code(base: u8) -> u8 {
        b"=ACMGRSVTWYHKDBN"
            .iter()
            .position(|b| *b == base)
            .unwrap_or(15) as u8
    }

    fn op(len: u32, kind: u32) -> u32 {
        (len << 4) | kind
    }

    fn any() -> EntryFilter {
        EntryFilter {
            chr_index: None,
            start: 0,
            end: None,
            standard_flags: false,
        }
    }

    fn decode(bytes: Vec<u8>, tags: bool) -> Vec<BamRecord> {
        decode_block(&Bytes::from(bytes), tags, &any(), &names(), "test").unwrap()
    }

    #[test]
    fn a_record_decodes_to_its_documented_fields() {
        let bytes = record(
            0,
            100,
            0x10,
            "read1",
            &[op(5, 4), op(80, 0)],
            b"ACGTA",
            &[30; 5],
            &[],
        );
        let records = decode(bytes, false);
        assert_eq!(records.len(), 1);
        let r = &records[0];
        assert_eq!(r.chr(), "chr1");
        assert_eq!(r.start(), 100);
        assert_eq!(r.end(), 180); // 80M consumes reference; 5S does not
        assert_eq!(r.read_name(), "read1");
        assert_eq!(r.cigar(), "5S80M");
        assert_eq!(r.sequence(), "ACGTA");
        // Phred 30 + 33 = 63 = '?'.
        assert_eq!(r.qualities(), "?????");
        assert_eq!(r.mapping_quality(), 60);
        assert_eq!(r.reference_length(), 80);
        assert_eq!(r.query_length(), 85);
        assert!(r.is_reverse() && !r.is_paired());
        assert_eq!(r.next_chr(), "*");
    }

    #[test]
    fn an_absent_sequence_and_absent_qualities_read_as_a_star() {
        let bytes = record(0, 10, 0, "r", &[op(10, 0)], b"", &[], &[]);
        let records = decode(bytes, false);
        assert_eq!(records[0].sequence(), "*");
        assert_eq!(records[0].qualities(), "*");

        // All-0xFF qualities are how the format says "none".
        let bytes = record(0, 10, 0, "r", &[op(4, 0)], b"ACGT", &[0xFF; 4], &[]);
        assert_eq!(decode(bytes, false)[0].qualities(), "*");
    }

    #[test]
    fn one_missing_quality_among_real_ones_is_a_star_in_place() {
        let bytes = record(
            0,
            10,
            0,
            "r",
            &[op(4, 0)],
            b"ACGT",
            &[30, 0xFF, 30, 30],
            &[],
        );
        let quals = decode(bytes, false)[0].qualities().to_string();
        assert_eq!(quals.chars().nth(1), Some('*'));
        assert_eq!(quals.len(), 4);
    }

    #[test]
    fn an_odd_length_sequence_does_not_decode_its_padding() {
        let bytes = record(0, 10, 0, "r", &[op(3, 0)], b"ACG", &[30; 3], &[]);
        assert_eq!(decode(bytes, false)[0].sequence(), "ACG");
    }

    #[test]
    fn tags_decode_to_their_declared_types_in_file_order() {
        let mut tags = Vec::new();
        tags.extend_from_slice(b"NMi");
        tags.extend_from_slice(&3i32.to_le_bytes());
        tags.extend_from_slice(b"RGZgroup1\0");
        tags.extend_from_slice(b"XAA");
        tags.push(b'x');
        tags.extend_from_slice(b"XFf");
        tags.extend_from_slice(&1.5f32.to_le_bytes());
        tags.extend_from_slice(b"XBBc");
        tags.extend_from_slice(&3u32.to_le_bytes());
        tags.extend_from_slice(&[1u8, 2, 253]);

        let bytes = record(0, 10, 0, "r", &[op(4, 0)], b"ACGT", &[30; 4], &tags);
        let records = decode(bytes, true);
        let got = records[0].tags().unwrap();
        assert_eq!(got[0], ("NM".into(), TagValue::Int(3)));
        assert_eq!(got[1], ("RG".into(), TagValue::Str("group1".into())));
        assert_eq!(got[2], ("XA".into(), TagValue::Char('x')));
        assert_eq!(got[3], ("XF".into(), TagValue::Float(1.5)));
        // 'c' is signed, so 253 reads back as -3.
        assert_eq!(got[4], ("XB".into(), TagValue::IntArray(vec![1, 2, -3])));
    }

    #[test]
    fn tags_off_means_no_tags_rather_than_an_error() {
        let mut tags = b"NMi".to_vec();
        tags.extend_from_slice(&3i32.to_le_bytes());
        let bytes = record(0, 10, 0, "r", &[op(4, 0)], b"ACGT", &[30; 4], &tags);
        assert!(decode(bytes, false)[0].tags().unwrap().is_empty());
    }

    #[test]
    fn a_record_with_malformed_tags_reads_fine_and_fails_at_tags() {
        // A 'Z' string with no terminator.
        let tags = b"RGZunterminated".to_vec();
        let bytes = record(0, 10, 0, "r", &[op(4, 0)], b"ACGT", &[30; 4], &tags);
        let records = decode(bytes, true);
        // Everything else about the record still works.
        assert_eq!(records[0].start(), 10);
        assert_eq!(records[0].sequence(), "ACGT");
        let err = records[0].tags().unwrap_err().to_string();
        assert!(err.contains("unterminated"), "{err}");
    }

    #[test]
    fn a_long_cigar_is_read_from_its_cg_tag_and_the_tag_is_kept() {
        // The placeholder: <l_seq>S<ref_len>N, with the real ops in CG:B,I.
        let real = [op(4, 0), op(6, 2)]; // 4M6D — 10 reference bases
        let mut tags = b"CGBI".to_vec();
        tags.extend_from_slice(&(real.len() as u32).to_le_bytes());
        for o in real {
            tags.extend_from_slice(&o.to_le_bytes());
        }
        let placeholder = [op(4, 4), op(10, 3)]; // 4S10N
        let bytes = record(0, 50, 0, "r", &placeholder, b"ACGT", &[30; 4], &tags);

        let records = decode(bytes, true);
        assert_eq!(records[0].cigar(), "4M6D", "the placeholder was returned");
        assert_eq!(records[0].end(), 60);
        // The CG field is left where it is, unlike htslib which strips it.
        assert!(records[0]
            .tags()
            .unwrap()
            .iter()
            .any(|(tag, _)| tag == "CG"));
    }

    #[test]
    fn the_standard_filter_drops_what_it_documents() {
        let filter = EntryFilter {
            chr_index: None,
            start: 0,
            end: None,
            standard_flags: true,
        };
        assert!(filter.accepts(0, 10, 20, 0x002)); // proper pair, mapped
        assert!(filter.accepts(0, 10, 20, 0)); // unpaired, mapped
        assert!(!filter.accepts(0, 10, 20, 0x004)); // unmapped
        assert!(!filter.accepts(0, 10, 20, 0x001)); // paired, not properly
        assert!(!filter.accepts(0, 10, 20, 0x100)); // secondary
        assert!(!filter.accepts(0, 10, 20, 0x800)); // supplementary
        assert!(!filter.accepts(0, 10, 20, 0x200)); // failed QC
        assert!(!filter.accepts(0, 10, 20, 0x400)); // duplicate
    }

    #[test]
    fn an_alignment_covering_no_reference_still_overlaps_a_window_on_it() {
        // end == start, as a placed unmapped read has. htslib's bam_endpos
        // treats a reference length of 0 as 1, and so does this.
        let filter = EntryFilter {
            chr_index: Some(0),
            start: 100,
            end: Some(200),
            standard_flags: false,
        };
        assert!(filter.accepts(0, 100, 100, 0), "a window starting on it");
        assert!(filter.accepts(0, 150, 150, 0));
        assert!(!filter.accepts(0, 99, 99, 0));
        assert!(!filter.accepts(0, 200, 200, 0));
    }

    #[test]
    fn an_empty_region_overlaps_nothing() {
        let filter = EntryFilter {
            chr_index: Some(0),
            start: 100,
            end: Some(100),
            standard_flags: false,
        };
        assert!(!filter.accepts(0, 100, 200, 0));
    }

    #[test]
    fn a_truncated_record_is_corrupt_not_a_panic() {
        let mut bytes = record(0, 10, 0, "r", &[op(4, 0)], b"ACGT", &[30; 4], &[]);
        bytes.truncate(bytes.len() - 6);
        let err = decode_block(&Bytes::from(bytes), false, &any(), &names(), "test")
            .unwrap_err()
            .to_string();
        assert!(err.contains("truncated bam record"), "{err}");
    }

    #[test]
    fn a_record_declaring_fields_past_its_own_length_is_refused() {
        let mut bytes = record(0, 10, 0, "r", &[op(4, 0)], b"ACGT", &[30; 4], &[]);
        // Claim 500 cigar operations in a record that holds one.
        bytes[16 + 4] = 244;
        bytes[17 + 4] = 1;
        let err = decode_block(&Bytes::from(bytes), false, &any(), &names(), "test")
            .unwrap_err()
            .to_string();
        assert!(err.contains("do not fit"), "{err}");
    }

    #[test]
    fn several_records_in_one_block_all_decode() {
        let mut bytes = record(0, 10, 0, "a", &[op(4, 0)], b"ACGT", &[30; 4], &[]);
        bytes.extend(record(0, 20, 0, "b", &[op(4, 0)], b"TGCA", &[31; 4], &[]));
        bytes.extend(record(1, 30, 0, "c", &[op(4, 0)], b"GGGG", &[32; 4], &[]));
        let records = decode(bytes, false);
        assert_eq!(records.len(), 3);
        assert_eq!(records[2].chr(), "chr2");
        assert_eq!(records[1].read_name(), "b");
    }
}