Error types MUST implement std::error::Error using thiserror and capture necessary context and provide causes.

Magic numbers (bit flags, format constants, operation codes) MUST be defined as named constants, not used as raw literals. BAM flags MUST use named constants (e.g., FLAG_UNMAPPED instead of 0x4). CIGAR operation types MUST use named constants (e.g., CIGAR_M instead of 0).

BAM flags MUST have a BamFlags newtype wrapping u16 with named predicate methods (e.g., is_unmapped(), is_reverse()). The public API (RecordRef) MUST expose BamFlags rather than raw u16 where appropriate.

CIGAR operations MUST have a CigarOpType enum with variants for all 9 SAM-spec operations (M, I, D, N, S, H, P, =, X). The enum MUST provide consumes_ref() and consumes_query() methods. Invalid operation codes MUST be represented as None via from_bam(u8) -> Option<CigarOpType>.

Every platform-specific optimization (SIMD, architecture-specific intrinsics) MUST have a scalar fallback that produces identical results on all platforms. The optimized and fallback paths MUST be tested with property-based tests that verify equivalence for arbitrary inputs.

IndexedReader::open(path) MUST auto-detect the file format by inspecting magic bytes:

After format detection, the reader MUST locate the appropriate index file:

If the format is detected but no matching index is found, the error MUST name the expected index extension and suggest the tool to create it (samtools index for BAM/CRAM, samtools index or tabix for SAM).

The unified reader MUST be an enum dispatching to format-specific readers, not a trait object. This avoids dynamic dispatch overhead in the hot path and keeps the type concrete for fork().

The unified reader MUST expose:

For a given BAM file and its SAM/CRAM representations of the same data, fetch_into for the same region MUST produce records with equivalent logical content: same positions, flags, sequences, qualities, and CIGAR operations. Aux tags MUST contain the same set of tag names and values, but tag ordering and integer type codes (e.g., c vs i for small values) MAY differ between formats because BAM writers choose specific integer widths that SAM text cannot preserve.

Currently RecordStore::push_raw() takes raw BAM bytes. For SAM and CRAM, records arrive as parsed fields rather than BAM binary.

push_fields() and push_raw() MUST produce identical SlimRecord fixed fields and identical name/bases/cigar/qual slab contents for the same logical record. Aux tag slab contents MAY differ in integer type codes (a SAM-derived i:42 may serialize as c while the original BAM used C), since the original BAM writer's type choice is not recoverable from SAM text. Tests MUST verify fixed-field and non-aux-slab equivalence by pushing the same record both ways.

SAM headers are text (the @HD, @SQ, @RG, @PG, @CO lines). CRAM containers embed SAM header text. BamHeader already stores header_text: String and parses @SQ lines for targets.

Indexed random access assumes coordinate-sorted data. If the @HD header line contains SO:unsorted or SO:queryname, the reader MUST return an error explaining that indexed region queries require coordinate-sorted input. SO:coordinate or absent SO (common in older files) MUST be accepted.

BAM fork: shares Arc<BamShared> (index + header), opens fresh File handle. (Already implemented.)

SAM fork: shares Arc holding parsed tabix/CSI index + header, opens fresh File handle for BGZF reading.

CRAM fork: shares Arc holding CRAI index + header, opens fresh File handle for container reading. Each fork MUST have its own FASTA reader (via fasta_reader.fork()) for thread-safe reference lookups. Reference caching (r[cram.perf.reference_caching]) MUST be per-fork, not shared.

The Readers struct bundles an IndexedReader (alignment) with an IndexedFastaReader (reference) in a single type. This eliminates the need for separate FASTA path passing — CRAM can access the reference it needs for sequence reconstruction, and all formats have uniform open/fork semantics.

Readers::open(alignment_path, fasta_path) MUST auto-detect the alignment format (via r[unified.detect_format]), open the appropriate reader, and open the FASTA reader. For CRAM, the fasta_path is passed to IndexedCramReader::open() for sequence reconstruction. For BAM/SAM, the FASTA reader is opened but not used by the alignment reader.

Readers::fork() MUST fork both the alignment reader and the FASTA reader, returning a new Readers with independent I/O handles but shared immutable state (indices, headers). The CRAM fork gets its own FASTA reader via IndexedFastaReader::fork().

Readers MUST expose:

IndexedReader::open(path) MUST continue to work for BAM and SAM files without a FASTA path. CRAM detection in IndexedReader::open() MUST return an error explaining that CRAM requires a reference and suggesting Readers::open() instead. This preserves backward compatibility for code that only needs BAM/SAM.

A BGZF block MUST begin with the gzip magic bytes 1f 8b 08 04 (gzip, DEFLATE, FEXTRA flag set).

The extra field MUST contain a BC subfield (SI1=0x42, SI2=0x43) with SLEN=2 whose value is the total block size minus one (BSIZE). The total compressed block size is BSIZE + 1.

Each block's DEFLATE payload MUST be decompressed independently. The last 8 bytes of the block are the gzip footer: a CRC32 checksum (4 bytes) followed by the uncompressed size ISIZE (4 bytes), both little-endian. Decompression MUST produce exactly ISIZE bytes.

BGZF blocks MUST NOT exceed 65536 bytes uncompressed. If the ISIZE footer field claims a larger value, the reader MUST return an error rather than allocating an unbounded buffer.

An EOF marker block has ISIZE=0. When encountered, the reader MUST signal end-of-stream.

A virtual offset is a 64-bit value where the upper 48 bits encode the compressed block offset in the file and the lower 16 bits encode the byte offset within the uncompressed block.

The reader MUST support seeking to an arbitrary virtual offset by seeking the underlying file to the block offset, decompressing that block, and advancing to the within-block offset.

The reader MUST support reading an exact number of bytes, transparently crossing block boundaries when the requested data spans multiple blocks.

The reader MUST support partial reads (up to N bytes) returning the actual count, returning 0 at EOF.

Decompression SHOULD use the libdeflater crate (Rust bindings to the libdeflate C library) for performance parity with htslib's libdeflate usage. libdeflate provides hardware-accelerated DEFLATE decompression and CRC32 computation.

After decompression, the reader MUST verify the CRC32 checksum of the decompressed data against the expected CRC32 stored in the gzip footer (4 bytes before ISIZE). A mismatch MUST return a ChecksumMismatch error.

The reader SHOULD fast-path standard BGZF headers where XLEN=6 and the BC subfield is at the fixed offset (bytes 12–17 of the 18-byte header). All BAM files produced by samtools, htslib, and Picard use this layout. This avoids allocating an extra-fields buffer for the common case. Non-standard layouts MUST fall back to searching the extra fields for the BC subfield.

When resizing buffers that will be immediately and fully overwritten (by read_exact or decompression), the reader SHOULD use an uninitialized resize to avoid redundant zero-filling. This applies to the decompressed block buffer (~64 KB per block) and the compressed data buffer.

The reader MUST track the current block's compressed file offset for virtual offset calculation. The offset MUST be queried from the stream position before reading each new block, and set directly on seeks.

The following CIGAR operations MUST be supported:

calc_matches_indels(cigar) MUST compute the total number of matching bases (sum of M operation lengths) and the total number of indel bases (sum of I and D operation lengths) from a packed CIGAR array.

A CigarIndex MUST be constructable from a record's CIGAR and start position. It precomputes cumulative reference and query offsets per operation.

CigarIndex::qpos_at(ref_pos) MUST return Some(query_pos) when the reference position falls within a query-consuming operation (M, =, X, I at the boundary), or None when the position falls within a deletion, reference skip, or outside the alignment.

qpos_at MUST return None for positions outside the alignment's reference span. This applies to all mapping variants, including the linear fast-path. Positions before the alignment start or at/after the alignment end MUST return None.

qpos_at MUST produce identical results to htslib's Alignment::qpos() for every reference position within the alignment span. This is verified by comparison tests.

Phred scores MUST be non-negative. from_phred MUST accept only u8 values (0..=255), making negative inputs unrepresentable at the type level.

as_int MUST clamp the result to the range [0, 99]. Values less than or equal to zero (including negative results from floating-point arithmetic) MUST return 0. NaN MUST return 0.

The reader MUST open a BAM file by path, parse its header, and locate the corresponding .bai index file (trying both <path>.bai and <path_without_ext>.bai).

The reader MUST expose the parsed BamHeader for contig name/tid/length lookups.

fetch_into(region, store) MUST clear the store, query the BAI index for chunks overlapping the region, seek through those chunks via BGZF, decode all records, and add records that actually overlap the region to the store (RecordStore).

A record overlaps a region if record.pos <= region.end AND record.end_pos >= region.start. Records outside the region that appear in index chunks MUST be filtered out. This is necessary because the BAM index is approximate — chunks may contain records slightly outside the queried region.

Records in the store MUST be in coordinate-sorted order (by pos, then by end_pos as tiebreaker) after fetch_into, matching the order a sorted BAM file provides.

Unmapped reads (flag 0x4) MUST be excluded from the store during fetch_into. htslib's pileup engine rejects them in bam_plp_push, so including them would inflate depth counts.

Secondary (0x100) and supplementary (0x800) reads MUST be included in the store by default. The caller's pileup filter is responsible for excluding them if desired, matching htslib's behavior.

The BAM index and header MUST be stored behind Arc so that forked readers share a single copy. These structures are read-only after initial parsing and are accessed concurrently without synchronization.

fork() MUST produce a new reader that:

A forked reader MUST produce identical fetch_into results as independently opening the same BAM file. The fork is purely an optimization — it MUST NOT change observable behavior.

Forked readers MUST be fully independent for mutable operations. Seeking or fetching on one reader MUST NOT affect any other reader (original or forked). Each reader owns its own file descriptors and buffer state.

Arc::ptr_eq on the shared state of the original and forked reader MUST return true, confirming they share the same allocation rather than having made a deep copy.

Multiple forks from the same prototype MUST be safe to use concurrently from different threads. Since the shared state is immutable and each fork owns its file handles, no synchronization is needed beyond the Arc reference count.

fetch_into MUST validate that tid fits in i32 and that start/end fit in i64 before using them. If any value exceeds the signed maximum, fetch_into MUST return a BamError::CoordinateOverflow rather than silently truncating via as casts.

A BAI file MUST begin with the magic bytes BAI\1 (0x42, 0x41, 0x49, 0x01). The reader MUST reject files that do not match.

The reader MUST parse the BAI binary format: n_ref references, each containing bins with chunks (pairs of virtual offsets) and a linear index of 16 KiB window offsets.

When decompressing BGZF blocks (e.g. for tabix), the reader MUST validate that bsize is large enough for the block header and footer before indexing into the block. A bsize that is too small MUST return an error rather than panicking.

Given a region (tid, start, end), the index MUST return a list of BGZF chunk pairs (begin, end) that contain all reads overlapping the region. The query MUST use the binning scheme to identify candidate bins and the linear index to skip chunks that cannot contain reads starting before the query start.

The BAI binning scheme uses reg2bin(beg, end) where bin 0 spans the whole reference, bins 1–8 span 64 Mbp each, and leaf bins span 16 KiB. The implementation MUST correctly compute bins for a region and enumerate all overlapping bins at every level.

The query MUST always include bin 0 in the candidate bins for any region query. Omitting it would silently drop reads that straddle 64 Mbp boundaries.

The index SHOULD provide a query_split method that separates distant (level 0–2) chunks from nearby (level 3–5) chunks in the result, returning a QueryChunks struct with nearby and distant fields. This allows callers to handle them differently (e.g. caching, separate I/O). The original query method MUST remain available for callers that don't need the separation.

The query_split method MUST separate chunks into nearby (levels 3–5, loaded via RegionBuf per query) and distant (levels 0–2, suitable for caching). The threshold (level ≤ 2) is based on profiling data showing that levels 0–2 contribute >90% of the file span in typical region queries.

The reader SHOULD maintain a ChunkCache that loads ALL records from bins at levels 0–2 (covering ≥8 Mbp) once per reference sequence (tid) per thread. The BamIndex::distant_chunks method returns all level 0–2 chunks for a tid, unfiltered by region, for use when populating the cache. Per query, matching records (overlapping the query region) are injected from the cache into the RecordStore without any additional I/O. The cache MUST be invalidated when switching to a different reference (tid). This eliminates expensive seeks for distant-bin chunks on every region query.

A BAM file MUST begin with the magic bytes BAM\1 (0x42, 0x41, 0x4d, 0x01) after BGZF decompression. The reader MUST reject files that do not start with this magic.

After the magic, a 4-byte little-endian l_text gives the length of the SAM header text. The reader MUST parse this text to extract @HD, @SQ, @RG, and @PG header lines.

Header length fields (l_text, n_ref, l_name) are stored as i32 in the BAM format. The reader MUST validate that these values are non-negative before using them as sizes. Negative values from corrupt files MUST return a BamHeaderError rather than causing a massive allocation via i32 as usize wrapping.

After the header text, a 4-byte little-endian n_ref gives the number of reference sequences. Each reference is encoded as a 4-byte l_name (including NUL), the name bytes, and a 4-byte l_ref (sequence length). The reader MUST parse all reference sequences.

The header MUST support O(1) lookup from contig name to target ID (tid) and O(1) lookup from tid to contig name and length.

The header MUST provide access to the list of all target names in tid order.

The header MUST provide the length of each target sequence by tid.

A BAM record MUST be decoded from its binary representation: 32 fixed bytes followed by variable-length qname, CIGAR, sequence, quality, and auxiliary data. The 4-byte block_size prefix is read by the caller.

BAM records have a 2 MiB size cap. The reader MUST reject records whose block_size exceeds this limit rather than allocating an unbounded buffer from untrusted input.

Variable-length field offset calculations (var_start, cigar_end, seq_end, qual_end) MUST use checked arithmetic to prevent wrapping on malformed input. Overflow MUST return a decode error.

The decoder MUST extract: tid (i32), pos (i32→i64), mapq (u8), flags (u16), seq_len (u32), n_cigar_ops (u16), qname (NUL-terminated, stored without NUL), packed sequence (4-bit), quality scores, CIGAR operations (packed u32), and raw auxiliary data.

The decoder MUST precompute the reference end position from the CIGAR. Reference-consuming operations are: M(0), D(2), N(3), =(7), X(8). The end position is pos + ref_consumed - 1 (inclusive).

When a read has zero reference-consuming CIGAR operations (e.g., pure soft-clip or insertion-only), the reference span is zero. In this case end_pos MUST equal pos (matching htslib's bam_cigar2rlen which returns 0, placing the read at exactly one position). The implementation MUST NOT compute pos - 1 for such reads.

BAM encodes sequences in 4-bit nibbles: two bases per byte (high nibble first). The standard lookup table maps 0→=, 1→A, 2→C, 4→G, 8→T, 15→N, and other values to IUPAC ambiguity codes.

The record MUST support decoding a single base at a given read position.

The decoder MAY decode the full 4-bit sequence into Base enum values at construction time using SIMD-accelerated bulk decoding (see base_decode.md). Single-base access via seq_at(idx, pos) returns the pre-decoded Base value directly.

is_reverse() MUST return true when the 0x10 bit is set.

is_first_in_template() MUST return true when the 0x40 bit is set.

is_second_in_template() MUST return true when the 0x80 bit is set.

is_unmapped() MUST return true when the 0x4 bit is set.

The record MUST support looking up auxiliary tags by their 2-byte name. Tag types A, c, C, s, S, i, I, f, d, Z, H, and B (typed array) MUST be supported.

The record MUST provide access to raw auxiliary data bytes for efficient filtering without full tag parsing.

The scalar decoder MUST use a 16-entry lookup table mapping 4-bit codes to ASCII bases (=ACMGRSVTWYHKDBN). It MUST support decoding full sequences (all bases) and single bases at arbitrary positions.

For bulk decoding, a 256-entry pair table (DECODE_PAIR[byte] → [high_base, low_base]) SHOULD be used to decode two bases per byte lookup, avoiding repeated shift-and-mask operations.

On platforms with SIMD support (SSSE3 on x86_64, NEON on aarch64), bulk sequence decoding MUST use SIMD table-lookup instructions (pshufb / vqtbl1q_u8) to decode 32 bases per iteration from 16 packed bytes. A scalar tail MUST handle the remaining bytes.

Runtime dispatch MUST select the fastest available decoder: SSSE3 on x86_64 (with feature detection), NEON on aarch64 (always available), scalar fallback on all other platforms.

The scalar encoder MUST use a 256-entry lookup table mapping ASCII bytes to 4-bit codes. Unknown bases MUST map to 15 (N). Two bases MUST be packed per byte (high nibble first).

The SIMD decoder MUST produce byte-identical output to the scalar decoder for all inputs. This MUST be verified by property-based tests covering arbitrary sequence lengths and byte patterns, including SIMD boundary lengths (0, 1, 15, 16, 31, 32, 33, 63, 64, 65, 128).

A 16-entry lookup table MUST map 4-bit BAM nibble codes directly to Base discriminant values: 1→A(65), 2→C(67), 4→G(71), 8→T(84), and all other codes (0, 3, 5, 6, 7, 9–14, 15) → Unknown(78).

Every entry in DECODE_BASE_TYPED (16 entries) and DECODE_PAIR_TYPED (256 × 2 = 512 bytes) MUST be a valid Base discriminant: A(65), C(67), G(71), T(84), or Unknown(78). This MUST be verified by an exhaustive test covering all entries.

decode_bases(packed_seq, seq_len) MUST decode a 4-bit packed BAM sequence into a Vec<Base>. The function SHOULD use the SIMD-accelerated decode path (SSSE3/NEON) with a Base-valued lookup table for identical throughput to the ASCII decoder.

The RecordStore MUST store decoded bases in a dedicated slab (Vec<Base>) separate from the data slab. The data slab continues to hold cigar, qual, and aux data. The seq_at method MUST return Base instead of u8.

PileupAlignment.base MUST be of type Base, not u8. This eliminates the Base::from(u8) conversion in the metrics accumulation hot loop.

Base::from_ascii_vec(Vec<u8>) -> Vec<Base> MUST convert a vector of ASCII bytes to Base values in-place (reusing the allocation). Maps A/a→A, C/c→C, G/g→G, T/t→T, all other bytes→Unknown. The function MUST NOT allocate a new vector — it operates on the input's buffer and transmutes the result. This is safe because Base is repr(u8) and the function only writes valid Base discriminant bytes (65, 67, 71, 78, 84).

The ASCII batch converter MUST use SIMD acceleration (SSSE3 on x86_64, NEON on aarch64) with a scalar fallback, following the same dispatch pattern as decode_bases. This is the required path for all bulk u8→Base conversions outside BAM 4-bit decoding: FASTA reference sequences, SAM text sequences, and any other ASCII-encoded base data.

The SIMD ASCII converter MUST produce identical output to applying Base::from(u8) element-wise. This MUST be verified by property-based tests covering arbitrary byte values and lengths including SIMD boundary lengths (0, 1, 15, 16, 31, 32, 33, 63, 64, 65, 128).

FromStr for Base MUST accept only a single character A, C, G, T, or N (case-insensitive) after trimming whitespace. Multi-character input (after trimming) MUST return Err(BaseError::MultipleChars). Any single character that is not A/C/G/T/N MUST return Err(BaseError::InvalidBase(byte)) where byte is the first byte of the trimmed input. Empty input MUST return Err(BaseError::Empty).

BaseError variants MUST display the actual invalid value. InvalidBase(byte) SHOULD format as "Invalid base: 0x{byte:02x}". Empty SHOULD format as "Empty".

The reader MUST support CRAM v3.0 and v3.1. CRAM v2.0 and v2.1 MAY be supported but are not required — they are rarely encountered in practice. The magic bytes CRAM followed by major version 3 MUST be accepted; major version 2 SHOULD produce a warning; other versions MUST be rejected.

CRAM decoding requires access to the reference genome (FASTA). The reader MUST accept a FASTA reader (the one being built on the other branch) for sequence lookups. Embedded reference slices (where the reference is stored in the CRAM container itself) MUST also be supported as a fallback.

The file MUST start with CRAM (bytes 43 52 41 4d), followed by major version (1 byte) and minor version (1 byte), followed by a 20-byte file ID.

The first container after the file definition MUST be the file header container. It contains a single block (content_type=FILE_HEADER) whose data is:

The file MUST end with an EOF container. The EOF container produced by samtools has length=15 (not 0) and num_records=0. The reader MUST detect the EOF container as any container where num_records=0 and length <= 15. The reader SHOULD verify the EOF marker but MUST NOT fail if it's missing (truncated files should produce a warning, not a hard error).

Each container header contains:

During fetch_into, the reader MUST skip containers whose [alignment_start, alignment_start + alignment_span) range does not overlap the query region. The CRAI index provides byte offsets for this.

Each slice header contains:

When reference_sequence_id == -2, the slice contains records from multiple references. The RI (reference ID) data series MUST be decoded per-record to determine which reference each record aligns to. Multi-ref slices are common in coordinate-sorted CRAM files at chromosome boundaries and for reads with unmapped mates.

When embedded_reference >= 0, the corresponding external data block contains the reference bases for this slice. The reader MUST use these instead of the FASTA reader. This enables CRAM files to be self-contained (produced by samtools view --output-fmt-option embed_ref or cramtools).

Each block contains:

CRC32 MUST be verified on every block in v3.0+. The CRC32 covers all bytes from the start of the block (method byte) through the end of the compressed data, excluding the 4 CRC32 bytes themselves. Mismatches MUST produce an error with the block's content_type and content_id for diagnostics.

ITF8 (up to 32 bits) MUST be decoded as:

LTF8 (up to 64 bits) MUST be decoded with the same leading-bit pattern extended to 9 bytes, with 0xFF prefix indicating 8 continuation bytes.

The core data block (content_type=5) is read as a bit stream. Bits are read MSB-first (most significant bit first) within each byte:

The preservation map MUST be decoded as a series of key-value pairs:

The substitution matrix is 5 bytes, one per reference base in order: A, C, G, T, N. Each byte encodes the ordering of substitution alternatives using 2-bit codes packed MSB-first:

The data series encoding map MUST be decoded as key-value pairs where:

The tag encoding map uses ITF8 keys computed as (char1 << 16) + (char2 << 8) + type_char. Each key maps to an encoding descriptor for that tag type. Tag values are stored as raw BAM-format bytes (without the 3-byte tag+type prefix).

Encoding ID 0 (NULL): no data. Used for absent data series.

Encoding ID 1 (EXTERNAL): reads raw bytes from an external data block identified by content_id. The most common encoding — used for quality scores, bases, tags, etc.

Encoding ID 3 (HUFFMAN): canonical Huffman coding. Parameters: alphabet (itf8 array) and bit lengths (itf8 array). Reads from the core data bit stream.

Encoding ID 4 (BYTE_ARRAY_LEN): a length encoding followed by a value encoding. The length is decoded first (from any encoding), then that many bytes are decoded using the value encoding.

Encoding ID 5 (BYTE_ARRAY_STOP): reads bytes from an external block until a stop byte is encountered. Used for read names (stop=0x00) and variable-length strings.

Encoding ID 6 (BETA): fixed-width integer. Parameters: offset and number of bits. Reads from core bit stream.

Encoding ID 7 (SUBEXP): sub-exponential code. Parameters: offset (itf8) and K (itf8). Reads from core bit stream. Decode procedure:

Encoding ID 9 (GAMMA): Elias gamma code with offset. Parameters: offset (itf8). Reads from core bit stream. Decode procedure:

Method 0 (RAW): uncompressed data.

Method 1 (GZIP): deflate compression. MUST be supported — this is the most common codec.

Method 2 (BZIP2): bzip2 compression. MUST be supported for v3.0 compatibility.

Method 3 (LZMA): LZMA compression. MUST be supported for v3.0 compatibility.

Method 4 (rANS 4×8): rANS entropy coder with 4-way interleaved states and 8-bit renormalization (L=2^23). Supports order-0 (context-free) and order-1 (previous-symbol context) modes. MUST be supported — this is the default codec in CRAM v3.0 writers including samtools.

Method 5 (rANS Nx16): v3.1 codec. rANS with N-way interleaving (N=4 or 32), 16-bit renormalization (L=2^15), and optional transforms controlled by an 8-bit flags byte: ORDER(1), N32(4), STRIPE(8), NoSize(16), CAT(32), RLE(64), PACK(128). Transforms are applied in order: decode entropy → unRLE → unPack. MUST be supported for v3.1 files — samtools uses this for most data blocks in v3.1 output.

Method 8 (tok3): v3.1 read-name tokeniser. MUST be supported for v3.1 files — samtools uses tok3 for read name blocks by default in v3.1 output. Without tok3 support, v3.1 CRAM files produced by samtools view -C cannot be read.

Unknown codec methods MUST produce a clear error naming the method ID and suggesting conversion to BAM (samtools view -b).

Records within a slice MUST be decoded in a strict field order, because the core bit stream is sequential. Reading fields out of order corrupts all subsequent data. The order is:

BF (BAM flags) and CF (CRAM flags) MUST be decoded per record. BAM flags map directly to the u16 flags in SlimRecord. CRAM flags are separate and encode:

AP (alignment position): if the preservation map AP=true (delta mode), each record's position is a delta from the previous record's position within the slice. The reader MUST maintain a running position accumulator and convert to absolute 0-based coordinates (CRAM positions are 1-based). Delta-coded positions can be negative (valid for supplementary alignments in coordinate-sorted CRAM).

RL (read length): the number of bases in the read. Used to allocate sequence and quality arrays. Long reads (PacBio/ONT) can be 10,000-1,000,000+ bases — buffer pre-allocation MUST NOT assume short-read sizes.

RG (read group): an index into the @RG entries in the header. Stored as itf8. The reader MUST convert to a RG:Z:<id> aux tag with the read group ID string.

RN (read name): if preservation map RN=true, names are stored. If RN=false, names are auto-generated. Rastair uses qnames for mate-pair dedup, so RN=false CRAM files MUST produce a warning that dedup will not work correctly. Generated names SHOULD be deterministic (e.g., cram_<record_counter>).

For each record, FN gives the number of read features. For each feature:

Read sequences are NOT stored directly. They MUST be reconstructed from the reference and read features:

Quality scores are decoded from QS (per-base) or QQ (stretch) data series, depending on CRAM flags bit 0x1. If CRAM flags bit 0x8 is set (sequence unknown), quality is unavailable (all 0xFF). Quality length MUST equal read length.

CRAM does not store CIGAR explicitly. The CIGAR MUST be reconstructed from the read features:

When CRAM flags bit 0x2 is set (detached), mate information is stored inline: MF (mate flags), NS (next ref ID), NP (next mate position), TS (template size). These fields are NOT needed for RecordStore (which doesn't store mate info), but the decoder MUST read and discard them to advance the data stream correctly.

When bit 0x2 is clear, the mate is "attached" — NF (next fragment distance) gives the record-skip to the mate within the same slice. The decoder MUST read NF but does not need to resolve the linkage for RecordStore purposes.

TL (tag line index) selects which tag combination this record has from the tag dictionary in the preservation map. For each tag in the combination, the tag encoding map provides the encoding for its value. Tag values MUST be decoded and serialized to BAM binary aux format for storage in the data slab.

The RG data series is separate from aux tags. If a read group is present, the reader MUST emit an RG:Z:<id> aux tag using the read group ID from the header's @RG entries. This tag MUST be included in the aux slab alongside dictionary-decoded tags.

The .crai file is gzip-compressed TAB-delimited text. Each line has 6 fields:

Region queries MUST find all index entries where [start, start+span) overlaps the query region for the given reference ID. The reader then seeks to the container byte offset, skips to the slice offset, and decodes the slice.

CRAI entries with alignment_span == 0 (but alignment_start > 0) indicate unknown extent — this occurs when samtools writes CRAM with embedded references or when the span is not computed. These entries MUST be included in query results whenever entry.alignment_start < query_end, since the slice may contain records at any position past the start. Treating span=0 as "zero-width" would silently drop all records in that slice.

Multi-ref slices produce multiple index entries (one per reference they span). A query for one reference may hit a multi-ref slice that also contains records from other references. The reader MUST filter records by reference ID after decoding.

Unmapped records (ref ID = -1) have start=0, span=0. They are only returned if explicitly queried. For fetch_into, unmapped records MUST be skipped (same as BAM).

If the reference MD5 in the slice header does not match the MD5 of the FASTA sequence for that region, the reader MUST return an error. The MD5 is computed over the uppercase reference sequence for the exact range [alignment_start, alignment_start + alignment_span), with no newlines or other formatting. If the range extends beyond the reference length, this is a file/reference mismatch error.

If the FASTA reader cannot provide the reference sequence for a slice's region (e.g., contig not in FASTA), the reader MUST return an error unless the slice has an embedded reference. The error SHOULD mention REF_PATH/REF_CACHE as alternatives.

When CRAM flags bit 0x8 is set, the read's sequence is unknown (* in SAM). seq_len MUST be 0 and no bases are pushed to the bases slab.

Slices with 0 records MUST be silently skipped.

Delta-coded positions can produce negative values if records are not strictly sorted. The decoder MUST handle this gracefully (negative deltas are valid for supplementary alignments within a coordinate-sorted CRAM).

Unmapped reads (flag 0x4) may appear in the last containers of the file, in slices with reference_sequence_id = -1. Their bases are stored in the BA data series (no reference reconstruction). fetch_into MUST skip them, same as BAM.

PacBio and ONT CRAM files contain reads of 10,000-1,000,000+ bases. Read feature lists can have thousands of entries, and quality/sequence arrays are correspondingly large. Buffer pre-allocation MUST scale with RL, not assume short-read sizes.

Query coordinates passed to fetch_into are u64, but internal CRAM positions are i32/i64. When converting u64::MAX to i64 for overlap filtering, the cast wraps to -1, which silently filters out all records. The reader MUST clamp coordinates: end.min(i64::MAX as u64) as i64.

The rANS 4x8 frequency table uses run-length encoding where a symbol counter sym: u8 is incremented for each run element. When sym == 255, sym += 1 overflows. The reader MUST use wrapping arithmetic (wrapping_add(1)) for this counter — the overflow is harmless because the loop terminates before using the overflowed value.

Frequency normalization (normalize_frequencies) sums all 256 frequency entries and doubles the sum in a loop. Both the sum and the doubling MUST use checked arithmetic to detect overflow from malformed frequency tables. Overflow MUST produce an error, not silent wraparound.

The uint7 variable-length integer decoding loop MUST be bounded to at most 5 iterations. A u32 can hold at most 32 bits; 5 iterations of 7-bit groups produce 35 bits, which is the maximum meaningful input. More than 5 continuation bytes indicate malformed data and MUST produce an error.

Length fields decoded from ITF8 (e.g., num_content_ids, num_blocks, alignment_span) may be negative when interpreted as i32. Before using such values as usize for allocation or iteration, the reader MUST validate they are non-negative via i32::try_from or equivalent. Negative values MUST produce an error, not wrap to huge usize values causing OOM.

The TokenReader::get() and get_mut() methods MUST return the same reader for every TokenType variant. In particular, DZLen MUST map to dz_len_reader in both methods. A mismatch causes the dup-copy path (which uses get()) to read from the wrong stream.

The name_count field in tok3 headers comes from untrusted data. Allocation based on name_count MUST be bounded to a reasonable limit (e.g., 10,000,000 names). The per-name token vector MUST be grown dynamically rather than pre-allocated to a fixed size, since the number of token positions varies per name.

CRAM random access is slice-granular, not record-granular. A region query may decompress an entire slice (typically 10,000 records) to extract a few records at the edges. This is inherently coarser than BAM's record-level seeking.

Reference sequence lookups happen per-slice (to reconstruct all records in the slice). The reader SHOULD cache the most recently used reference region to avoid repeated FASTA lookups for consecutive slices on the same contig. This cache MUST be per-fork, not shared across threads.

rANS and arithmetic decoders have higher per-byte CPU cost than zlib. The reader SHOULD pre-allocate decode buffers and reuse them across blocks within a slice.

The FAI parser MUST read a .fai file and produce an index with one entry per sequence. Each entry MUST contain: sequence name (String), length (number of bases, u64), byte offset of the first base in the file (u64), bases per line (u64), and bytes per line including line terminator (u64). The parser MUST NOT trim whitespace from field values — tabs are the only delimiter, and sequence names may contain spaces.

Each line MUST have exactly 5 tab-separated fields: NAME, LENGTH, OFFSET, LINEBASES, LINEWIDTH. Lines with fewer or more fields MUST produce an error identifying the malformed line. Empty lines MUST be skipped.

The parser MUST reject entries where linewidth < linebases (line terminator would be negative), linebases == 0 (would cause division by zero in offset calculations), or length == 0 (empty sequences are not useful for region queries). Note: linewidth == linebases is valid — it means lines have no trailing newline character, which occurs when a sequence is stored on a single line without a mid-sequence newline.

The index MUST support O(1) lookup of a sequence entry by name. Duplicate sequence names MUST produce an error at parse time.

For a FASTA file at path P, the reader MUST look for the index at P.fai (i.e., appending .fai to the full path including any .gz extension). This matches samtools/htslib convention: ref.fa → ref.fa.fai, ref.fa.gz → ref.fa.gz.fai.

When the FAI index file is not found, the reader MUST produce a clear error message that names the expected index path and suggests running samtools faidx <fasta_path> to create it. The reader MUST NOT attempt to build the index itself. Rationale: auto-indexing hides user errors (wrong FASTA passed silently), fails on read-only shared reference directories (common on HPC clusters), introduces non-deterministic file creation and race conditions, and takes 10–30s for a human genome which makes the tool appear to hang with no feedback.

To find the byte offset of base at 0-based position pos within a sequence, the reader MUST compute: offset + (pos / linebases) * linewidth + (pos % linebases), where offset, linebases, and linewidth come from the FAI entry. This accounts for newline characters inserted at line boundaries.

fetch_seq(name, start, stop) MUST use 0-based half-open coordinates [start, stop). A request for (name, 0, 4) returns the first 4 bases of the named sequence.

The reader MUST verify that start < stop (zero-length and reversed ranges are invalid) and stop <= sequence_length. Out-of-bounds requests MUST produce an error naming the sequence, the requested range, and the actual sequence length.

Returned sequence bytes MUST be uppercased. FASTA files may contain lowercase bases (indicating soft-masking), but Seqair always works with uppercase bases.

When reading from the file, the reader MUST strip line terminators (LF \n and CR \r) from the raw bytes. The returned sequence MUST contain only base characters.

When the requested sequence name is not in the index, the reader MUST produce an error listing the requested name. It SHOULD include available sequence names if the index is small (fewer than 20 entries) to help diagnose typos.

fetch_seq and fetch_seq_into MUST return raw uppercase ASCII bytes (Vec<u8> / &mut Vec<u8>), not Base values. CRAM reference MD5 verification computes checksums on the exact byte sequence from the FASTA file — converting to Base and back is lossy because IUPAC ambiguity codes (M, R, S, etc.) collapse to Unknown. Any u8→Base conversion MUST happen downstream, after CRAM verification is complete.

fetch_seq_into(name, start, stop, buf) MUST accept a caller-provided &mut Vec<u8> buffer, clear it, and write the fetched sequence into it. This avoids per-call allocation when fetching many segments. The fetch_seq convenience method MAY allocate a new Vec<u8> internally.

For uncompressed FASTA files, the reader MUST seek to the computed byte offset and read the required byte range in a single I/O operation, then strip newlines and uppercase the result.

The raw read size MUST account for embedded newlines. For a fetch of n bases starting at position pos, the number of raw bytes to read is: end_byte_offset - start_byte_offset where both offsets are computed via fasta.index.offset_calculation.

The reader MUST auto-detect bgzip compression by reading the first 18 bytes of the file and checking for: (1) gzip magic bytes 1f 8b, (2) DEFLATE compression method 08, (3) FEXTRA flag 04, and (4) the BC subfield identifier at bytes 12–13. All four checks MUST pass to identify a file as BGZF. Plain gzip files with the FEXTRA flag but without the BC subfield MUST NOT be misidentified as BGZF.

For bgzip FASTA files, a GZI index (.gzi) MUST be present alongside the FASTA. The GZI index maps compressed block offsets to uncompressed byte offsets, enabling random access into the compressed data.

When the GZI index is not found for a bgzip FASTA, the reader MUST produce a clear error naming the expected GZI path and suggesting samtools faidx <fasta_path> (which creates both .fai and .gzi).

The GZI parser MUST read a binary file containing: an 8-byte little-endian entry count (u64), followed by that many pairs of 8-byte little-endian values: (compressed_offset: u64, uncompressed_offset: u64). Each pair records where a BGZF block starts in the compressed file and what uncompressed byte offset that corresponds to. The entry count multiplication MUST use checked arithmetic to prevent overflow on corrupted files.

GZI entries MUST be sorted in strictly increasing order by uncompressed_offset. The parser MUST reject files with unsorted or duplicate entries, as translate() relies on binary search.

Given an uncompressed byte offset, the GZI index MUST be able to find the BGZF block that contains it: find the last entry whose uncompressed_offset <= target, then the compressed block starts at that entry's compressed_offset, and the within-block offset is target - uncompressed_offset.

The within-block offset MUST fit in u16 (0–65535), matching the BGZF specification's maximum uncompressed block size of 65,536 bytes. If the computed within-block offset exceeds u16::MAX, the reader MUST return an error indicating the GZI index may be corrupt or missing entries. The reader MUST NOT silently truncate the offset.

The reader MUST decompress BGZF blocks using the same libdeflater-based decompression as the BAM reader, including CRC32 verification.

For a fetch spanning multiple BGZF blocks, the reader MUST read and decompress blocks sequentially from the starting block until all requested bytes are obtained. Since FASTA fetches are contiguous byte ranges (unlike BAM which has disjoint index chunks), a simple sequential scan suffices.

The FAI index and GZI index (if applicable) MUST be stored behind Arc so forked readers share a single parsed copy. These structures are read-only after initial parsing.

fork() MUST produce a new reader that shares the same Arc as the source (no re-parsing) and opens a fresh File handle. Each fork MUST own its own read buffer.

A forked reader MUST produce byte-identical results to independently opening the same FASTA file.

Forked readers MUST be fully independent for mutable operations. Seeking on one MUST NOT affect any other.

All FASTA-specific errors MUST contain context (file paths, sequence names, positions) sufficient to diagnose the problem without a debugger. The implementation MUST NOT use panic!, unwrap(), expect(), or unreachable!() — all error paths MUST propagate errors via Result.

The SAM reader MUST open a bgzf-compressed SAM file, parse its header, and locate a tabix (.tbi) or CSI (.csi) index file. If neither index exists, return an error suggesting samtools index or tabix -p sam.

If the file does not start with BGZF magic (1f 8b with BC subfield), the reader MUST reject it with an error explaining that plain SAM files must be bgzf-compressed for indexed access (bgzip file.sam).

The SAM header consists of all lines starting with @ at the beginning of the file. The reader MUST:

If no @SQ lines are present, the reader MUST return an error. Without @SQ lines, target name→tid mapping is impossible and fetch_into cannot work.

The full header text (all @ lines joined with newlines) MUST be stored in BamHeader::header_text for downstream use (e.g., constructing htslib Header for BAM writing).

Each alignment line MUST be split on TAB characters into the 11 mandatory fields plus optional tag fields. The parser MUST handle:

SAM POS is 1-based; the internal representation is 0-based (matching BAM). The parser MUST subtract 1 from POS. A POS of 0 in SAM means unmapped — after subtraction this becomes -1 (as i64). Unmapped records are filtered by FLAG 0x4 before reaching this point, so -1 positions should not appear in the RecordStore.

The CIGAR string MUST be parsed into the same packed u32 format used by BAM: each operation is (length << 4) | op_code, where op_codes are M=0, I=1, D=2, N=3, S=4, H=5, P=6, =7, X=8. The string * means CIGAR unavailable (0 operations, and end_pos = pos).

SEQ characters MUST be decoded to Base values: A→A, C→C, G→G, T→T, all others (N, IUPAC codes)→Unknown. The string * means sequence unavailable (length 0).

QUAL characters MUST have 33 subtracted (Phred+33 encoding). The string * means quality unavailable; all quality bytes MUST be set to 0xFF (255). QUAL length MUST equal SEQ length when both are present.

Optional fields after the 11th column follow the format TAG:TYPE:VALUE. They MUST be serialized to BAM binary aux format for storage in the data slab:

Malformed aux tag values MUST return an error, not silently default to zero. For integer (i) typed tags, if the value cannot be parsed as an integer, the reader MUST return a SamRecordError::InvalidAuxValue error.

SAM aux integer values (type i) are serialized into the smallest BAM integer type that fits, using the ranges defined in r[sam.record.aux_tags]. Values that exceed the u32 range (i.e., val > u32::MAX or val < i32::MIN when interpreted as i64) cannot be represented in any BAM integer type and MUST return a SamRecordError::AuxIntOutOfRange error instead of silently truncating.

fetch_into(tid, start, end, store) MUST:

The overlap filter uses half-open intervals (0-based). end_pos MUST be computed from the parsed CIGAR. When CIGAR is * (unavailable), end_pos = pos (point record).

A record overlaps [start, end) iff pos < end && end_pos > start. Records where pos == end or end_pos == start do NOT overlap and MUST be filtered out.

Records MUST be added to the store in coordinate-sorted order. Since bgzf-compressed indexed SAM files are coordinate-sorted by construction (indexing requires it), the natural parse order is correct.

The reader MUST support tabix (.tbi) indexes. Tabix is a generic index for bgzf-compressed, TAB-delimited files. For SAM, the format code is 1 (not 0/generic). The column configuration (seq_col=3, beg_col=4, end_col=0) is implicit for format=1 and may be ignored by some implementations; the reader SHOULD NOT reject a tabix file based on stored column values when format=1 is present.

Index file lookup order: <path>.tbi, <path>.bai. If neither exists, return a clear error. CSI (.csi) lookup should be added when CSI support is implemented.

When SEQ is *, the record has no sequence data. seq_len MUST be 0 and no bases are pushed to the bases slab. This is valid per the SAM spec for secondary alignments that omit sequence.

When QUAL is *, all quality bytes MUST be 0xFF. This matches BAM's representation of unavailable quality.

BAM encodes CIGARs with >65535 operations using a CG:B:I aux tag and a placeholder CIGAR. SAM has no such limitation — the CIGAR string can be arbitrarily long. The parser MUST handle CIGARs of any length.

Blank lines (empty or whitespace-only) between alignment records MUST be silently skipped.

RNAME of * means unmapped. Combined with FLAG 0x4, these records MUST be filtered out by fetch_into (same as BAM's unmapped filtering). RNAME * without FLAG 0x4 is malformed but SHOULD be treated as unmapped.

A single SAM alignment line may span multiple BGZF blocks. The reader MUST buffer partial lines across block boundaries and only parse complete lines (terminated by \n). The parser MUST also strip \r at line boundaries to handle \r\n (Windows-style) line endings that may appear in files from mixed-platform pipelines.

Integer parsing functions (parse_u8, parse_u16, parse_u32) MUST reject empty input by returning None. An empty byte slice is not a valid integer representation and MUST NOT silently produce zero.

Like BAM, the SAM reader SHOULD use RegionBuf-style bulk I/O: read all compressed bytes for the region in one large I/O operation, then decompress and parse from memory. This is critical for NFS/Lustre performance.

Tabix (.tbi) files are BGZF-compressed. The parser MUST decompress via BGZF before parsing the binary format. Reading raw bytes without decompression produces garbage.

Tabix files start with magic bytes TBI\x01 (4 bytes), after BGZF decompression.

After magic:

Since tabix uses identical binning parameters to BAI (min_shift=14, depth=5), the index data after the tabix header MUST be parsed by reusing BamIndex::from_reader() (or equivalent). This avoids duplicating the bin/chunk/linear-index parsing code. The tabix header is read and validated first, then the remaining bytes are handed to the BAI parser.

After the header, the index data follows the same structure as BAI:

Bin ID 37450 is a pseudo-bin that stores metadata (virtual offset range of all records for a reference). Its chunks MUST NOT be included in query results. The reader MUST recognize and skip this bin during region queries.

Region queries use the same reg2bins algorithm as BAI (with min_shift=14, depth=5). The bin offset for level l (counting from root, l=0=root) is ((1 << 3*l) - 1) / 7, and the bin ID for position p at level l is:

The existing BamIndex::query() returns Vec<Chunk> with begin/end virtual offsets. The tabix and CSI parsers MUST produce the same Vec<Chunk> output so that RegionBuf::load() can be reused unchanged. The binning parameters differ but the output format is identical.

Chunk and VirtualOffset types from the BAM index module MUST be reused by tabix/CSI indexes.

A reference with no records has no bins and no linear index entries. Queries for such references MUST return an empty chunk list (not an error).

On construction, the store MUST pre-allocate all three vectors with estimated capacities to avoid reallocation during push_raw. Estimates SHOULD be based on the number of compressed bytes loaded for the region: a typical BAM compression ratio of ~3:1, an average record size of ~400 bytes, ~25 bytes per read name, and ~375 bytes per record of non-name data. The store MUST accept a byte-count hint for this purpose.

push_raw(raw_bytes) MUST decode a BAM record's fixed fields and append its variable-length data directly into the slabs, without allocating intermediate Box<[u8]> or Vec<u8> per field. The 4-bit packed sequence MUST be decoded (via SIMD when available) to Base enum values in the bases slab.

Slab offsets (name_off, bases_off, data_off) MUST be checked for u32 overflow before storing. If any slab exceeds u32::MAX bytes, the store MUST return an error rather than silently truncating the offset.

push_fields(...) MUST accept pre-parsed record fields for SAM and CRAM readers: pos, end_pos, flags, mapq, matching_bases, indel_bases, qname bytes, CIGAR as packed BAM-format u32 ops, sequence as &[Base], quality bytes, and aux tag bytes in BAM binary format. This avoids encoding to BAM binary only to immediately decode it again. SAM and CRAM parsers convert their native representations to these types before pushing.

The store MUST provide methods to access each variable-length field for a given record index: qname(idx), cigar(idx), seq(idx), seq_at(idx, pos), qual(idx), aux(idx). These return borrowed slices into the slabs.

Clearing the store for a new region MUST retain the allocated capacity of all four vectors, so that subsequent regions reuse the same memory without reallocation.

All records in the store MUST remain valid and accessible until the store is cleared.

The pileup engine MUST NOT use Rc for record sharing. PileupAlignment carries pre-extracted flat fields (base, qual, mapq, flags, etc.), and ActiveRecord stores a record index into the store. The store is either owned by or borrowed by the engine for the duration of iteration.

Before reading, index chunks MUST be merged into non-overlapping byte ranges sorted by file offset. Overlapping or adjacent chunks MUST be coalesced into a single range. The end of each range MUST extend past the last chunk's end block offset by at least one maximum BGZF block size (64 KiB) to ensure the final block is fully included in the buffer.

RegionBuf MUST NOT include bin 0 chunks in its bulk read. Bin 0 chunks are at distant file offsets and contain very few reads; including them would cause an expensive additional seek. Bin 0 records are handled separately via the bin 0 cache (see bam_index.md).

RegionBuf::load MUST accept a seekable reader and a slice of index chunks. It MUST merge chunks, then perform one seek + read_exact per merged range to load all compressed bytes into a contiguous in-memory buffer. On typical BAM files this results in a single read of a few hundred KB to a few MB per region.

When given zero chunks, load MUST return an empty buffer that immediately signals EOF on any read.

The region buffer MUST support seeking to a virtual offset within the pre-loaded data. Since disjoint ranges are concatenated, the buffer MUST maintain a range map that translates file offsets to buffer positions. A seek to a file offset that falls within any loaded range MUST succeed; a seek to an offset in a gap between ranges or before/after the loaded data MUST return an error.

Block decompression from the in-memory buffer MUST follow the same BGZF format rules as the file-based reader: magic validation, BSIZE extraction, DEFLATE decompression, and CRC32 verification.

The fast-path header parsing (XLEN=6, BC at fixed offset) MUST be used when applicable, with fallback to searching extra fields.

The region buffer MUST support reading an exact number of bytes, transparently crossing block boundaries, with the same API as BgzfReader::read_exact_into.

The Drop implementation MUST never panic. Gap-size calculations between ranges MUST use saturating arithmetic to handle overlapping or malformed range boundaries safely.

The region buffer MUST track virtual offsets correctly: the block offset corresponds to the file position of the current block (not the buffer position), so that virtual offset comparisons with the BAM index chunk boundaries work correctly.

Both IndexedBamReader::fetch_into and IndexedSamReader::fetch_into MUST use RegionBuf to bulk-read compressed data before record decoding. The file-based BgzfReader is retained for header parsing and other sequential access at file open time.

CRAM does NOT use RegionBuf. CRAM containers are not BGZF-compressed — the reader seeks to container byte offsets directly via the CRAI index and reads container data with plain File::read_exact. Block-level decompression (gzip, rANS, etc.) happens after reading the container into memory.

Overlapping pair deduplication MUST be opt-in via set_dedup_overlapping(). When not enabled, the pileup engine MUST yield all alignments without deduplication.

Mates are detected by matching qnames. When set_dedup_overlapping() is called, the engine MUST build a mate index that maps each record to its mate's index (if both are in the arena). Records with no mate in the arena are unaffected.

Only the first two records sharing a qname are treated as mates. If three or more records share the same qname (supplementary alignments, etc.), only the first pair is linked; additional records are treated as unpaired for dedup purposes.

When both mates show the same base at a position, the first encountered (by arena order) is kept and the other is removed.

When mates show different bases at a position, the first-in-template read is kept and the second-in-template read is removed. This preserves the primary observation.

Deduplication MUST be applied per-position, not globally. A read that is removed at one overlapping position MUST still appear at positions where its mate is not active. Both mates contribute normally outside their overlap region.

Deduplication MUST be applied after the read filter and after qpos computation. A read that was excluded by the filter MUST NOT participate in dedup. A read that has no qpos at a position (deletion/skip) MUST NOT participate in dedup at that position.

Deduplication MUST be applied before max_depth truncation. The depth after dedup may be below max_depth even if the raw depth exceeded it.

The engine MUST iterate over every reference position from the region start to the region end (inclusive), yielding a PileupColumn at each position where at least one read is active.

At each position, the active set contains all records whose alignment span [pos, end_pos] includes the current position. Records MUST enter the active set when the current position reaches their pos and MUST be evicted when the current position exceeds their end_pos.

Each PileupColumn MUST provide:

The engine MUST support a per-read filter function that is evaluated once per record when it enters the active set (NOT re-evaluated at each position). Records that fail the filter MUST be excluded from all subsequent columns.

The engine MUST support a maximum depth setting. When more records pass the filter than max_depth allows, excess records MUST be dropped (preferring records already in the active set).

Max depth MUST be enforced per-position, not globally at arena-load time. A read that is rejected at one high-coverage position may still be included at an adjacent lower-coverage position. This matches htslib's behavior where MAXCNT is checked against the current column's depth.

For each alignment in a column, the engine MUST provide qpos: the position within the read's query sequence that aligns to the current reference position. This MUST be computed using the CigarIndex.

When the current reference position falls within a deletion or reference skip in the read's CIGAR, that alignment MUST be excluded from the column (no qpos exists).

Positions with zero active reads MUST be skipped (no column yielded). When there is a gap between reads, the engine MUST jump to the next read's start position rather than iterating through empty positions one by one.

Reads with the unmapped flag (0x4) MUST be excluded from the pileup. htslib's bam_plp_push rejects unmapped reads before they enter the buffer.

Reads with zero reference-consuming CIGAR operations (pure soft-clip, insertion-only) have end_pos == pos and MUST appear in exactly one pileup column at their pos. They MUST NOT be skipped or cause off-by-one errors.

When a reference position falls within the soft-clipped portion of a read's alignment, the read is NOT active at that position (soft clips do not consume reference). The read only becomes active at its first reference-consuming position.

For any BAM file and region, the pileup engine MUST produce identical columns to htslib's bam_plp_auto-based pileup: same positions, same depth at each position, and same qpos values per alignment. This is verified by comparison tests.

The pileup engine allocates a Vec<PileupAlignment> per column using Vec::with_capacity(active.len()). Reusing a buffer via clone was measured to be slower than fresh allocation due to the copy cost of entries. The allocator efficiently reuses recently-freed blocks of the same size.

When a record enters the pileup active set, the engine MUST retrieve the record data once and cache what it needs (cigar, flags, strand, mapq, etc.) in the ActiveRecord — NOT perform repeated store lookups for the same record in the hot loop.

Since BAM records arrive in coordinate-sorted order from fetch_into, the pileup engine MUST NOT sort record indices. It MUST iterate records in arena order directly.

Building a CigarIndex MUST NOT clone the CIGAR bytes via .to_vec(). The CIGAR bytes live in the arena slab for the lifetime of the region; the CigarIndex MUST borrow or reference them by offset, not own a copy.

Matches and indels counts MUST be computed once per record during decode and stored in BamRecord, NOT recomputed from CIGAR at every pileup position.

For records with more than 4 CIGAR operations, qpos_at SHOULD use binary search on ref_start instead of linear scan. For 1–4 ops, linear scan is acceptable.

The RecordStore MUST support a capacity hint so that the first region's allocation can be pre-sized based on an estimate (e.g., from the BAI index chunk sizes), avoiding repeated reallocation during the first fetch_into. This is provided by RecordStore::with_byte_hint().