gwseq-io 0.2.1

//! Reference sequences, for rebuilding the bases of mapped reads.
//!
//! Format reference: `docs/cram_format_v3.1.md` §11.
//!
//! CRAM stores a mapped read as its *differences* from the reference, so the
//! matching stretches — which is nearly all of them — are not in the file. This
//! module is where they come from.
//!
//! # What needs one, and what does not
//!
//! Only `SEQ`. Coordinates, CIGARs, flags, names, mates, tags, soft-clipped and
//! inserted bases are either stored verbatim or derivable from the read
//! features alone, so a coverage or depth read never touches this module. A
//! reader with no reference reports why through `CramReader::reference_error`
//! and reads sequences as `N`, the same shape `BamReader::index_error` already
//! has for a missing index.
//!
//! # Where a reference comes from
//!
//! In order: the `reference` open option; the slice's own embedded reference
//! block, when it carries one; the `UR` field of the header's `@SQ` line; and
//! htslib's `REF_CACHE` / `REF_PATH` conventions, which every CRAM-producing
//! pipeline already sets. The last of those looks a sequence up by the `M5`
//! checksum the header states — the string, not one computed here.
//!
//! # On not checking M5
//!
//! §11 asks readers to check the `M5` checksums and report mismatches. This one
//! does not, and says so rather than implying otherwise: verifying a checksum
//! means reading and hashing the whole chromosome — 195 Mbp for human chr1 —
//! at open, for every reference in the header. That is a different cost from
//! opening a file. What the checksum is *used* for here is finding a sequence
//! in a reference cache, which needs the header's string and no hashing.

use std::collections::HashMap;
use std::sync::Arc;

use bytes::Bytes;
use parking_lot::Mutex;

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

/// The smallest window worth reading.
///
/// The caller asks for exactly what a slice covers, so this is a floor rather
/// than a size: it keeps a small request from costing a read per record, while
/// a large one is served at its own size instead of being rounded up.
///
/// It used to round every request up to eight megabases from its start, which
/// is a fine trade for a walk down a chromosome and a poor one for the
/// library's main use — scattered loci, where every locus lands in a different
/// window and pays a whole one. Measured over 500 random kilobase loci that was
/// about eight milliseconds each, most of the difference between reading a CRAM
/// with a reference and without.
const MIN_WINDOW: i64 = 1 << 20;

/// How many windows to keep, at least.
///
/// One would do for a single-reference slice, which walks forward through one
/// chromosome. A **multi-reference** slice is why this is more than one: it
/// holds the tail of one reference and the head of the next, and its records
/// alternate between them at the boundary. With a single window that is a
/// reload per record.
///
/// The real floor, though, is the worker count: `parallel` workers on
/// different loci want `parallel` windows, and with fewer they evict each
/// other's on every locus. [`ReferenceSource::with_slots`] raises it.
const WINDOWS: usize = 3;

/// One line of a `.fai`: where a sequence sits and how it is wrapped.
#[derive(Debug, Clone, Copy)]
struct FaiEntry {
    length: i64,
    /// Byte offset of the sequence's first base, in the *uncompressed* file.
    offset: u64,
    /// Bases per line, and bytes per line including the line terminator.
    line_bases: i64,
    line_width: i64,
}

/// A bgzip index: where each block starts, compressed and uncompressed.
///
/// Written by `bgzip -i` beside a `.gz`, and the only way to reach the middle
/// of a compressed FASTA without inflating everything before it.
#[derive(Debug, Default)]
struct Gzi {
    /// `(compressed offset, uncompressed offset)`, ascending, with the implicit
    /// first entry `(0, 0)` included.
    blocks: Vec<(u64, u64)>,
}

impl Gzi {
    fn parse(data: &[u8], path: &str) -> Result<Self> {
        if data.len() < 8 {
            return Err(Error::format(
                path,
                "a gzi index too short to hold its count",
            ));
        }
        let count = u64::from_le_bytes(data[..8].try_into().expect("eight bytes"));
        // Checked, not because a real index is ever this large but because the
        // count is a number a file chose: `8 + count * 16` overflows for a
        // count the first eight bytes of anything at all can spell.
        let wanted = usize::try_from(count)
            .ok()
            .and_then(|count| count.checked_mul(16))
            .and_then(|bytes| bytes.checked_add(8));
        let Some(wanted) = wanted.filter(|wanted| *wanted <= data.len()) else {
            return Err(Error::format(
                path,
                format!("a gzi index naming {count} blocks in {} bytes", data.len()),
            ));
        };
        let count = count as usize;
        let mut blocks = Vec::with_capacity(count + 1);
        blocks.push((0u64, 0u64));
        for chunk in data[8..wanted].chunks_exact(16) {
            let compressed = u64::from_le_bytes(chunk[..8].try_into().expect("eight bytes"));
            let uncompressed = u64::from_le_bytes(chunk[8..].try_into().expect("eight bytes"));
            blocks.push((compressed, uncompressed));
        }
        Ok(Self { blocks })
    }

    /// The block holding an uncompressed offset, and how far into it that is.
    fn locate(&self, uncompressed: u64) -> (u64, u64) {
        // The last block starting at or before the offset.
        let index = self
            .blocks
            .partition_point(|(_, start)| *start <= uncompressed)
            .saturating_sub(1);
        let (compressed_start, uncompressed_start) = self.blocks[index];
        (compressed_start, uncompressed - uncompressed_start)
    }
}

/// A FASTA, indexed, and bgzip-compressed or not.
#[derive(Debug)]
pub struct Fasta {
    source: Arc<dyn ByteSource>,
    path: String,
    index: HashMap<String, FaiEntry>,
    /// Present when the FASTA is bgzip-compressed.
    gzi: Option<Gzi>,
}

impl Fasta {
    /// Open `path` with its `.fai`, and its `.gzi` when it is compressed.
    pub fn open(path: &str) -> Result<Self> {
        let source = crate::source::open(path, None, None)?;
        let index_path = format!("{path}.fai");
        let index_source = crate::source::open(&index_path, None, None)?;
        let index = parse_fai(&index_source.read_to_end(0)?, &index_path)?;

        // A bgzip-compressed FASTA is a gzip file, and reaching into one needs
        // the block index beside it. A plain FASTA needs neither.
        let head = source.read_at(0, 2)?;
        let gzi = if crate::source::is_gzipped(&head) {
            let gzi_path = format!("{path}.gzi");
            let gzi_source = crate::source::open(&gzi_path, None, None).map_err(|_| {
                Error::format(
                    path,
                    format!(
                        "this reference is compressed and {gzi_path} is not there, so its \
                         middle cannot be reached; `bgzip -r` writes one"
                    ),
                )
            })?;
            Some(Gzi::parse(&gzi_source.read_to_end(0)?, &gzi_path)?)
        } else {
            None
        };

        Ok(Self {
            source,
            path: path.to_string(),
            index,
            gzi,
        })
    }

    pub fn has(&self, name: &str) -> bool {
        self.resolve_name(name).is_some()
    }

    /// The name this FASTA knows a sequence by, given the name the CRAM's
    /// header uses.
    ///
    /// `SN:1` against `>chr1` is the commonest reference mismatch there is,
    /// and matching only exactly makes it silent: every base comes back `N`,
    /// with `reference_error` empty and a `reference_path` set, which is
    /// indistinguishable from a read with no substitutions. The rule is
    /// [`crate::genomic::chr::ChrMap`]'s, so a CRAM and a bigWig resolve the
    /// same names the same way.
    pub fn resolve_name(&self, name: &str) -> Option<&str> {
        if let Some((key, _)) = self.index.get_key_value(name) {
            return Some(key);
        }
        let toggled = match name.get(..3) {
            Some(prefix) if prefix.eq_ignore_ascii_case("chr") => name[3..].to_string(),
            _ => format!("chr{name}"),
        };
        self.index
            .get_key_value(&toggled)
            .map(|(key, _)| key.as_str())
    }

    /// The length the `.fai` gives for a sequence, under either name.
    pub fn length(&self, name: &str) -> Option<i64> {
        let key = self.resolve_name(name)?;
        self.index.get(key).map(|entry| entry.length)
    }

    /// Bases `[start, end)` of `name`, 0-based, upper-cased, newlines removed.
    ///
    /// A range past the end of the sequence is clipped rather than refused —
    /// §11.3 says reads may be mapped beyond the reference, and the record
    /// decoder turns anything missing into `N`.
    pub fn read(&self, name: &str, start: i64, end: i64) -> Result<Vec<u8>> {
        let name = self.resolve_name(name).unwrap_or(name);
        let entry = *self.index.get(name).ok_or_else(|| {
            Error::format(
                &self.path,
                format!("this reference has no sequence named {name}"),
            )
        })?;
        let start = start.clamp(0, entry.length);
        let end = end.clamp(start, entry.length);
        // A line width narrower than its base count, or no bases per line at
        // all, describes no layout at all. Returning an empty sequence turns
        // that into all-N with nothing said, which is the same silent-wrong-
        // answer shape as a reference that does not fit the file.
        if entry.line_bases <= 0 || entry.line_width < entry.line_bases {
            return Err(Error::format(
                &self.path,
                format!(
                    "the index for {name} says {} bases on a line of {} bytes, which describes \
                     no layout",
                    entry.line_bases, entry.line_width
                ),
            ));
        }
        if start == end {
            return Ok(Vec::new());
        }

        // Wrapped lines mean the bases of a range are not contiguous on disk:
        // this is where the first and last of them sit, terminators included.
        let byte_of = |base: i64| -> u64 {
            entry.offset
                + (base / entry.line_bases * entry.line_width + base % entry.line_bases) as u64
        };
        let first = byte_of(start);
        let last = byte_of(end - 1) + 1;
        let raw = self.read_raw(first, (last - first) as usize)?;

        let mut out = Vec::with_capacity((end - start) as usize);
        for byte in raw.iter() {
            // Everything that is not a line terminator is a base.
            if !byte.is_ascii_whitespace() {
                out.push(byte.to_ascii_uppercase());
            }
        }
        out.truncate((end - start) as usize);
        Ok(out)
    }

    /// `len` bytes of the uncompressed file at `offset`.
    fn read_raw(&self, offset: u64, len: usize) -> Result<Bytes> {
        let Some(gzi) = &self.gzi else {
            return self.source.read_at(offset, len);
        };
        use crate::bam::bgzf::{Chunk, VirtualOffset};
        let (begin_block, begin_within) = gzi.locate(offset);
        let (end_block, end_within) = gzi.locate(offset + len as u64);
        // A within-block offset cannot exceed a BGZF block's 64 KiB of output,
        // so both fit the 16 bits a virtual offset gives them.
        let chunk = Chunk {
            begin: VirtualOffset::new(begin_block, begin_within as u16),
            end: VirtualOffset::new(end_block, end_within as u16),
        };
        crate::bam::bgzf::decompress_chunk(self.source.as_ref(), chunk, &self.path)
    }
}

/// Parse a `.fai`: name, length, offset, bases per line, bytes per line.
fn parse_fai(data: &[u8], path: &str) -> Result<HashMap<String, FaiEntry>> {
    let mut out = HashMap::new();
    for (number, line) in data.split(|b| *b == b'\n').enumerate() {
        let line = line.strip_suffix(b"\r").unwrap_or(line);
        if line.is_empty() {
            continue;
        }
        let mut fields = line.split(|b| *b == b'\t');
        let mut next = |what: &str| -> Result<&[u8]> {
            fields.next().ok_or_else(|| {
                Error::format(
                    path,
                    format!("line {} of this index has no {what}", number + 1),
                )
            })
        };
        let name = String::from_utf8_lossy(next("name")?).into_owned();
        let number_at = |field: &[u8], what: &str| -> Result<i64> {
            std::str::from_utf8(field)
                .ok()
                .and_then(|s| s.parse::<i64>().ok())
                .ok_or_else(|| {
                    Error::format(
                        path,
                        format!("line {} of this index has an unreadable {what}", number + 1),
                    )
                })
        };
        let length = number_at(next("length")?, "length")?;
        let offset = number_at(next("offset")?, "offset")?;
        let line_bases = number_at(next("line length")?, "line length")?;
        let line_width = number_at(next("line width")?, "line width")?;
        out.insert(
            name,
            FaiEntry {
                length,
                offset: offset.max(0) as u64,
                line_bases,
                line_width,
            },
        );
    }
    if out.is_empty() {
        return Err(Error::format(path, "this index names no sequences"));
    }
    Ok(out)
}

/// The window of one chromosome currently held.
#[derive(Debug, Default)]
struct Window {
    name: String,
    /// 0-based half-open span of `bases` on that chromosome.
    start: i64,
    end: i64,
    bases: Arc<Vec<u8>>,
}

/// Where a reader gets its reference bases.
///
/// Holds a few windows, because coordinate-sorted slices walk forward through a
/// reference and re-reading a few megabases per slice is the difference between
/// htslib's `REF_CACHE`: one file per sequence, named by its `M5` checksum and
/// holding the bare bases.
///
/// Not a FASTA. There is no header line, no `.fai` and no line wrapping — the
/// whole file *is* the sequence — so it needs its own reader rather than
/// another `Fasta`. Handing the cache directory to `Fasta::open`, which is what
/// this used to do, looks for `<directory>.fai` and fails; the feature was
/// advertised in four places and could not have worked.
///
/// The name-to-checksum map comes from the file's own `@SQ` lines, which is the
/// only thing that connects a sequence name to a file in this layout.
#[derive(Debug)]
pub struct RefCache {
    root: String,
    checksums: HashMap<String, String>,
}

impl RefCache {
    pub fn open(root: &str, sequences: &[(String, Option<String>, Option<String>)]) -> Self {
        Self {
            root: root.to_string(),
            checksums: sequences
                .iter()
                .filter_map(|(name, _, m5)| Some((name.clone(), m5.clone()?)))
                .collect(),
        }
    }

    fn file(&self, name: &str) -> Option<String> {
        let m5 = self.checksums.get(name)?;
        let path = format!("{}/{m5}", self.root);
        std::path::Path::new(&path).exists().then_some(path)
    }

    pub fn has(&self, name: &str) -> bool {
        self.file(name).is_some()
    }

    pub fn length(&self, name: &str) -> Option<i64> {
        let path = self.file(name)?;
        std::fs::metadata(path).ok().map(|m| m.len() as i64)
    }

    /// Bases `[start, end)`, 0-based and upper-cased.
    ///
    /// Clipped past the end as `Fasta::read` is, for the same reason: §11.3
    /// lets a read be mapped beyond the reference.
    pub fn read(&self, name: &str, start: i64, end: i64) -> Result<Vec<u8>> {
        let Some(path) = self.file(name) else {
            return Ok(Vec::new());
        };
        let source = crate::source::open(&path, None, None)?;
        let length = source.len()? as i64;
        let start = start.clamp(0, length);
        let end = end.clamp(start, length);
        if start == end {
            return Ok(Vec::new());
        }
        let bytes = source.read_at(start as u64, (end - start) as usize)?;
        Ok(bytes.iter().map(|b| b.to_ascii_uppercase()).collect())
    }
}

/// Where a `ReferenceSource` reads from.
#[derive(Debug)]
enum Backend {
    Fasta(Fasta),
    Cache(RefCache),
}

/// a query and a wait.
#[derive(Debug)]
pub struct ReferenceSource {
    backend: Backend,
    windows: Mutex<Vec<Window>>,
    slots: usize,
}

impl ReferenceSource {
    pub fn open(path: &str) -> Result<Self> {
        Ok(Self {
            backend: Backend::Fasta(Fasta::open(path)?),
            windows: Mutex::new(Vec::new()),
            slots: WINDOWS,
        })
    }

    /// Open htslib's per-checksum cache, which needs the header's `@SQ` lines
    /// to know which file is which sequence.
    pub fn open_cache(root: &str, sequences: &[(String, Option<String>, Option<String>)]) -> Self {
        Self {
            backend: Backend::Cache(RefCache::open(root, sequences)),
            windows: Mutex::new(Vec::new()),
            slots: WINDOWS,
        }
    }

    /// Keep at least one window per worker.
    ///
    /// Three slots served four workers on four different loci by evicting one
    /// of them on every locus, so each worker re-read a whole window each
    /// time round.
    pub fn with_slots(mut self, slots: usize) -> Self {
        self.slots = self.slots.max(slots);
        self
    }

    pub fn has(&self, name: &str) -> bool {
        match &self.backend {
            Backend::Fasta(fasta) => fasta.has(name),
            Backend::Cache(cache) => cache.has(name),
        }
    }

    /// The declared length of a sequence, for the check at open.
    pub fn length(&self, name: &str) -> Option<i64> {
        match &self.backend {
            Backend::Fasta(fasta) => fasta.length(name),
            Backend::Cache(cache) => cache.length(name),
        }
    }

    fn read(&self, name: &str, start: i64, end: i64) -> Result<Vec<u8>> {
        match &self.backend {
            Backend::Fasta(fasta) => fasta.read(name, start, end),
            Backend::Cache(cache) => cache.read(name, start, end),
        }
    }

    /// Bases `[start, end)` of `name`, 0-based.
    ///
    /// Returns the window they sit in along with where it starts, so the caller
    /// can hand both to the record decoder without a copy. The window is at
    /// least the requested range and usually much more.
    pub fn window(&self, name: &str, start: i64, end: i64) -> Result<(Arc<Vec<u8>>, i64)> {
        // Clamped to the sequence, because a window can never reach past it:
        // `read` clips, so a request for the last `read_length` bases of a
        // chromosome asks for an end no window can satisfy and misses the
        // cache on every record. The records that do this are the ones in the
        // slice at a chromosome's end — which, in a sorted file with
        // multi-reference slices, is a slice per chromosome.
        let end = match self.length(name) {
            Some(length) => end.min(length),
            None => end,
        };
        {
            let mut held = self.windows.lock();
            if let Some(position) = held
                .iter()
                .position(|w| w.name == name && w.start <= start && w.end >= end)
            {
                // Move to the back, so what is evicted is what has gone
                // longest unused rather than what was read longest ago.
                let hit = held.remove(position);
                let out = (hit.bases.clone(), hit.start);
                held.push(hit);
                return Ok(out);
            }
        }
        // Read out to a whole window from the request's start, so a walk
        // forward through a chromosome reloads rarely.
        let from = start.max(0);
        let to = end.max(from + MIN_WINDOW);
        let bases = Arc::new(self.read(name, from, to)?);
        let mut held = self.windows.lock();
        held.push(Window {
            name: name.to_string(),
            start: from,
            end: from + bases.len() as i64,
            bases: bases.clone(),
        });
        while held.len() > self.slots {
            held.remove(0);
        }
        Ok((bases, from))
    }
}

/// Where to look for a reference, resolved at open.
///
/// The variants are the answers, in the order they are tried.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum Reference {
    /// A FASTA the caller named.
    Given(String),
    /// The `UR` field of an `@SQ` line.
    FromHeader(String),
    /// Found in `REF_CACHE` or on `REF_PATH`, by the `M5` the header states.
    FromCache(String),
    /// Nothing resolved, and why.
    None(String),
}

impl Reference {
    pub fn path(&self) -> Option<&str> {
        match self {
            Self::Given(path) | Self::FromHeader(path) | Self::FromCache(path) => Some(path),
            Self::None(_) => None,
        }
    }
}

/// Find a reference for a file whose `@SQ` lines carry `UR` and `M5`.
///
/// `given` is the caller's `reference` option and wins outright — including
/// when it does not exist, because a caller who names a file wants to be told
/// that file is missing rather than to have another quietly used instead.
pub fn resolve(
    given: Option<&str>,
    sequences: &[(String, Option<String>, Option<String>)],
) -> Reference {
    if let Some(path) = given {
        return Reference::Given(path.to_string());
    }

    // The `UR` of the first `@SQ` that has one and that exists.
    for (_, uri, _) in sequences {
        if let Some(uri) = uri {
            // `file:///abs`, `file://host/abs` and the single-slash
            // `file:/abs` that some tools emit all name a local path.
            let path = uri
                .strip_prefix("file://")
                .or_else(|| uri.strip_prefix("file:"))
                .unwrap_or(uri);
            if crate::source::is_url(path) || std::path::Path::new(path).exists() {
                return Reference::FromHeader(path.to_string());
            }
        }
    }

    // htslib's conventions: one file per sequence, named by its checksum.
    let checksums: Vec<&str> = sequences
        .iter()
        .filter_map(|(_, _, m5)| m5.as_deref())
        .collect();
    if !checksums.is_empty() {
        for variable in ["REF_CACHE", "REF_PATH"] {
            let Ok(value) = std::env::var(variable) else {
                continue;
            };
            for root in value.split(':') {
                // htslib spells the layout with `%s` and digit counts; the
                // common case is a plain directory holding one file per
                // checksum, which is what is honoured here.
                let root = root.trim_end_matches("/%s").trim_end_matches('/');
                if root.is_empty() || root.contains('%') {
                    continue;
                }
                if checksums
                    .iter()
                    .all(|m5| std::path::Path::new(&format!("{root}/{m5}")).exists())
                {
                    return Reference::FromCache(root.to_string());
                }
            }
        }
    }

    let has_uri = sequences.iter().any(|(_, uri, _)| uri.is_some());
    Reference::None(if sequences.is_empty() {
        "this file's header names no reference sequences".to_string()
    } else if has_uri {
        let named: Vec<&str> = sequences
            .iter()
            .filter_map(|(_, uri, _)| uri.as_deref())
            .take(1)
            .collect();
        format!(
            "no reference was given, and the one this file's header names ({}) is not there",
            named.first().copied().unwrap_or("")
        )
    } else {
        "no reference was given, and this file's header names none (no UR field on its @SQ lines)"
            .to_string()
    })
}

/// The `.fai` parser, for `crate::fuzz`, which runs every parser in the crate
/// over hostile bytes and is a sibling of this module rather than a descendant.
#[cfg(test)]
pub(crate) fn parse_fai_for_fuzz(data: &[u8], path: &str) -> Result<()> {
    parse_fai(data, path).map(|_| ())
}

/// The `.gzi` parser, likewise.
#[cfg(test)]
pub(crate) fn parse_gzi_for_fuzz(data: &[u8], path: &str) -> Result<()> {
    Gzi::parse(data, path).map(|_| ())
}

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

    const FAI: &[u8] = b"chr1\t10\t6\t5\t6\nchr2\t4\t24\t4\t5\n";

    #[test]
    fn a_fai_parses_into_its_five_columns() {
        let index = parse_fai(FAI, "test.fai").expect("an index");
        let chr1 = index.get("chr1").expect("chr1");
        assert_eq!(chr1.length, 10);
        assert_eq!(chr1.offset, 6);
        assert_eq!(chr1.line_bases, 5);
        assert_eq!(chr1.line_width, 6);
        assert!(index.contains_key("chr2"));
    }

    #[test]
    fn an_unreadable_fai_is_refused_rather_than_half_read() {
        assert!(parse_fai(b"", "test.fai").is_err());
        assert!(parse_fai(b"chr1\tten\t6\t5\t6\n", "test.fai").is_err());
        assert!(parse_fai(b"chr1\t10\t6\n", "test.fai").is_err());
    }

    /// A FASTA written by hand, read back through the wrapping arithmetic —
    /// which is the part of this module that is easy to get wrong by one.
    #[test]
    fn bases_are_read_across_the_lines_they_are_wrapped_onto() {
        // >chr1\nACGTA\nCGTAC\n>chr2\nTTTT\n
        let fasta = b">chr1\nACGTA\nCGTAC\n>chr2\nTTTT\n";
        let index = parse_fai(FAI, "test.fai").expect("an index");
        let source = crate::source::testing::MemorySource::new(fasta.to_vec());
        let reader = Fasta {
            source: Arc::new(source),
            path: "test.fa".to_string(),
            index,
            gzi: None,
        };
        assert_eq!(reader.read("chr1", 0, 10).expect("bases"), b"ACGTACGTAC");
        assert_eq!(reader.read("chr1", 0, 5).expect("bases"), b"ACGTA");
        // A range that straddles the line break is where an off-by-one shows.
        assert_eq!(reader.read("chr1", 3, 7).expect("bases"), b"TACG");
        assert_eq!(reader.read("chr1", 5, 10).expect("bases"), b"CGTAC");
        assert_eq!(reader.read("chr1", 9, 10).expect("bases"), b"C");
        assert_eq!(reader.read("chr2", 0, 4).expect("bases"), b"TTTT");
    }

    #[test]
    fn a_range_past_the_end_of_a_sequence_is_clipped_rather_than_refused() {
        let index = parse_fai(FAI, "test.fai").expect("an index");
        let source = crate::source::testing::MemorySource::new(
            b">chr1\nACGTA\nCGTAC\n>chr2\nTTTT\n".to_vec(),
        );
        let reader = Fasta {
            source: Arc::new(source),
            path: "test.fa".to_string(),
            index,
            gzi: None,
        };
        assert_eq!(reader.read("chr1", 8, 100).expect("bases"), b"AC");
        assert_eq!(reader.read("chr1", 50, 60).expect("bases"), b"");
        assert!(reader.read("chrX", 0, 1).is_err());
    }

    #[test]
    fn a_gzi_locates_an_offset_in_the_block_that_holds_it() {
        let mut data = Vec::new();
        data.extend_from_slice(&2u64.to_le_bytes());
        // Two blocks beyond the implicit first: uncompressed 100 and 250.
        data.extend_from_slice(&40u64.to_le_bytes());
        data.extend_from_slice(&100u64.to_le_bytes());
        data.extend_from_slice(&90u64.to_le_bytes());
        data.extend_from_slice(&250u64.to_le_bytes());
        let gzi = Gzi::parse(&data, "test.gzi").expect("a gzi");
        assert_eq!(gzi.locate(0), (0, 0));
        assert_eq!(gzi.locate(99), (0, 99));
        assert_eq!(gzi.locate(100), (40, 0));
        assert_eq!(gzi.locate(249), (40, 149));
        assert_eq!(gzi.locate(250), (90, 0));
        assert_eq!(gzi.locate(1000), (90, 750));
    }

    #[test]
    fn a_truncated_gzi_is_refused() {
        assert!(Gzi::parse(&[0u8; 4], "test.gzi").is_err());
        let mut data = Vec::new();
        data.extend_from_slice(&5u64.to_le_bytes());
        data.extend_from_slice(&[0u8; 16]);
        assert!(Gzi::parse(&data, "test.gzi").is_err());
    }

    #[test]
    fn a_ref_cache_reads_the_bare_sequence_its_checksum_names() {
        // htslib's layout: one file per `M5`, holding the bases and nothing
        // else — no header line, no `.fai`, no wrapping.
        let dir = std::env::temp_dir().join(format!("gwseq-refcache-{}", std::process::id()));
        std::fs::create_dir_all(&dir).expect("cache dir");
        let m5 = "d41d8cd98f00b204e9800998ecf8427e";
        std::fs::write(dir.join(m5), b"acgtACGTnnnn").expect("cache file");

        let sequences = vec![("chr1".to_string(), None, Some(m5.to_string()))];
        let source = ReferenceSource::open_cache(&dir.to_string_lossy(), &sequences);
        assert!(source.has("chr1"));
        assert_eq!(source.length("chr1"), Some(12));

        let (bases, from) = source.window("chr1", 4, 8).expect("a window");
        assert_eq!(from, 4);
        // Upper-cased, as the FASTA reader upper-cases.
        assert_eq!(&bases[..4], b"ACGT");

        // A sequence the cache does not hold is absent rather than an error.
        assert!(!source.has("chr2"));
        let _ = std::fs::remove_dir_all(&dir);
    }

    #[test]
    fn a_fai_line_layout_that_describes_nothing_is_refused() {
        let dir = std::env::temp_dir().join(format!("gwseq-badfai-{}", std::process::id()));
        std::fs::create_dir_all(&dir).expect("dir");
        let fasta = dir.join("r.fa");
        std::fs::write(&fasta, b">chr1\nACGT\n").expect("fasta");
        // A line width narrower than the bases it claims to hold.
        std::fs::write(dir.join("r.fa.fai"), b"chr1\t4\t6\t4\t2\n").expect("fai");
        let fasta = Fasta::open(&fasta.to_string_lossy()).expect("opens");
        let error = fasta.read("chr1", 0, 4).expect_err("no layout");
        assert!(error.to_string().contains("no layout"), "{error}");
        let _ = std::fs::remove_dir_all(&dir);
    }

    #[test]
    fn a_single_slash_file_uri_is_a_local_path() {
        let dir = std::env::temp_dir().join(format!("gwseq-uri-{}", std::process::id()));
        std::fs::create_dir_all(&dir).expect("dir");
        let fasta = dir.join("r.fa");
        std::fs::write(&fasta, b">chr1\nACGT\n").expect("fasta");
        let uri = format!("file:{}", fasta.to_string_lossy());
        let sequences = vec![("chr1".to_string(), Some(uri), None)];
        assert_eq!(
            resolve(None, &sequences),
            Reference::FromHeader(fasta.to_string_lossy().to_string())
        );
        let _ = std::fs::remove_dir_all(&dir);
    }

    #[test]
    fn an_explicitly_given_reference_wins_over_everything_else() {
        let sequences = vec![(
            "chr1".to_string(),
            Some("/does/not/exist.fa".to_string()),
            Some("abc".to_string()),
        )];
        assert_eq!(
            resolve(Some("/given.fa"), &sequences),
            Reference::Given("/given.fa".to_string())
        );
    }

    #[test]
    fn a_header_with_no_uri_says_so_rather_than_failing_vaguely() {
        let sequences = vec![("chr1".to_string(), None, None)];
        match resolve(None, &sequences) {
            Reference::None(why) => assert!(why.contains("UR"), "{why}"),
            other => panic!("{other:?}"),
        }
        match resolve(None, &[]) {
            Reference::None(why) => assert!(why.contains("names no reference"), "{why}"),
            other => panic!("{other:?}"),
        }
    }

    #[test]
    fn a_header_naming_a_missing_reference_names_it_in_the_reason() {
        let sequences = vec![(
            "chr1".to_string(),
            Some("/nowhere/genome.fa".to_string()),
            None,
        )];
        match resolve(None, &sequences) {
            Reference::None(why) => assert!(why.contains("/nowhere/genome.fa"), "{why}"),
            other => panic!("{other:?}"),
        }
    }
}