CRAM Reader

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:

• header_text_length (i32, fixed-width little-endian): byte length of the SAM header text
• header_text (UTF-8 bytes): the SAM header text (@HD, @SQ, @RG, etc.)

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:

• length (i32, fixed-width little-endian): total byte length of all blocks in the container. This is the only fixed-width field in the header; all subsequent fields use variable-length encoding.
• reference_sequence_id (itf8): reference ID, or -1 (unmapped), or -2 (multi-ref)
• alignment_start (itf8): 1-based start position on reference
• alignment_span (itf8): number of reference bases covered
• num_records (itf8): number of alignment records
• record_counter (ltf8): global 0-based record index (v3.0+)
• bases (ltf8): total bases in container (v3.0+)
• num_blocks (itf8): block count
• landmarks (itf8 array): count (itf8) followed by that many itf8 values; byte offsets from container data start to each slice header
• crc32 (4 bytes, fixed-width little-endian): CRC32 over all preceding header bytes (from length through end of landmarks)

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:

• reference_sequence_id (itf8): ref ID for this slice, -1 unmapped, -2 multi-ref
• alignment_start, alignment_span (itf8): genomic range covered
• num_records (itf8): records in this slice
• record_counter (ltf8): global record counter
• num_blocks (itf8): data blocks in this slice
• block_content_ids (itf8 array): content IDs of the blocks
• embedded_reference (itf8): -1 if no embedded reference, else block content ID
• reference_md5 (16 bytes): MD5 of the reference segment (all zeros if N/A)
• optional tags (BAM-format aux bytes)

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:

• method (byte): compression codec (see codecs section)
• content_type (byte): 0=FILE_HEADER, 1=COMPRESSION_HEADER, 2=SLICE_HEADER, 4=EXTERNAL, 5=CORE
• content_id (itf8): identifies which data series this block holds
• compressed_size, uncompressed_size (itf8)
• data (bytes): compressed payload
• crc32 (4 bytes): CRC32 of the block (v3.0+)

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:

• 0xxxxxxx: 1 byte, value = bits 6..0 (7 bits)
• 10xxxxxx + 1 byte: 2 bytes, value = (b0 & 0x3F) << 8 | b1 (14 bits)
• 110xxxxx + 2 bytes: 3 bytes, value = (b0 & 0x1F) << 16 | b1 << 8 | b2 (21 bits)
• 1110xxxx + 3 bytes: 4 bytes, value = (b0 & 0x0F) << 24 | b1 << 16 | b2 << 8 | b3 (28 bits)
• 1111xxxx + 4 bytes: 5 bytes, value = (b0 & 0x0F) << 28 | b1 << 20 | b2 << 12 | b3 << 4 | (b4 & 0x0F) (32 bits)

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 MUST be read MSB-first (most significant bit first) within each byte:

• Bit offset 0 is the most significant bit (0x80)
• Bit offset 7 is the least significant bit (0x01)
• Reading N bits shifts them out from the current position, crossing byte boundaries as needed

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

• RN (bool): whether read names are preserved (if false, names are generated)
• AP (bool): whether alignment positions are delta-coded (true) or absolute
• RR (bool): whether the reference is required for decoding
• SM (5 bytes): substitution matrix (see r[cram.compression.substitution_matrix])
• TD (bytes): tag dictionary — null-separated lists of 3-byte entries (2-char tag name + 1-char BAM type)

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:

byte = (alt0 << 6) | (alt1 << 4) | (alt2 << 2) | alt3

where each 2-bit value maps to a base: 0=A, 1=C, 2=G, 3=T. When decoding a substitution (BS data series value = code), the read base is \alt[code] from this ordering.

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

• Key: 2-byte data series code
• Value: encoding descriptor (encoding ID + parameters)

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:

1. Read unary: count consecutive 1 bits until a 0 bit is read. Call this count n.
2. If n == 0: read K bits from the stream as the value.
3. If n > 0: read n + K - 1 bits from the stream as remainder. The value is (1 << (n + K - 1)) - (1 << K) + remainder.
4. The final decoded value is value - offset.

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

1. Read unary: count consecutive 0 bits until a 1 bit is read. Call this count n.
2. If n == 0: the value is 0.
3. If n > 0: read n more bits from the stream as suffix. The value is (1 << n) | suffix - 1.
4. The final decoded value is value - offset.

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:

1. BF (BAM flags)
2. CF (CRAM flags)
3. RI (reference ID — only if slice is multi-ref, i.e., ref_seq_id == -2)
4. RL (read length)
5. AP (alignment position — delta or absolute per preservation map)
6. RG (read group)
7. QS (quality scores — only if CRAM flags bit 0x1 is clear, meaning quality stored per-record not per-base)
8. RN (read name — only if preservation map RN=true)
9. Mate information:
   • If CRAM flags bit 0x2 set (detached): MF, NS, NP, TS
   • If CRAM flags bit 0x4 set (mate downstream): NF
10. TL (tag line index)
11. FN (feature count), then for each feature: FC, FP, feature-specific data
12. MQ (mapping quality)
13. QS per-base (only if CRAM flags bit 0x1 is set — quality as array)

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:

• Bit 0x1: quality scores stored as array (vs per-record)
• Bit 0x2: detached (mate info stored inline)
• Bit 0x4: has mate downstream in same slice
• Bit 0x8: sequence not stored (unknown)

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:

1. Decode FC (feature code byte)
2. Decode FP (position within read, delta-coded from the previous feature's position within the same record — first feature's FP is absolute within the read, starting at 1)
3. Decode feature-specific data:

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

1. Fetch the reference sequence for [alignment_start, alignment_start + alignment_span) from the FASTA reader (or embedded reference block).
2. Initialize the read's bases array with RL entries.
3. Walk the reference and read in parallel, applying features in order:
   • Between features: copy reference bases to read bases (these are the matching segments).
   • X (substitution): look up substitution_matrix[ref_base][BS_code] for the read base.
   • I/i (insertion): insert bases from IN/BA into the read (these don't consume reference).
   • D (deletion): skip DL reference bases (these don't consume read).
   • S (soft clip): insert bases from SC into the read (no reference consumption).
   • B (base+quality): replace base and quality at the feature position.
   • N (ref skip): skip RS reference bases.
   • H (hard clip): record length only (no bases in read or reference).
4. Convert reconstructed bases to Base enum values for the bases slab.

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:

• Stretches of matching reference → M ops (or =/X if distinguishing is needed, but M is standard)
• X (substitution) → part of M ops (substitutions are still alignment matches)
• I/i → I ops
• D → D ops
• S → S ops
• N → N ops
• H → H ops
• P → P ops

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:

1. Reference sequence ID (int)
2. Alignment start (int, 1-based matching the container header convention, 0 for unmapped)
3. Alignment span (int, 0 for unmapped)
4. Container byte offset in file (int, 0-based byte offset from start of file)
5. Slice byte offset relative to container data start (int)
6. Slice size in bytes (int)

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.

During sequence reconstruction, if the reference sequence is shorter than the positions accessed (i.e., ref_offset + ref_pos >= reference_seq.len()), the reader MUST log a warning on the first occurrence per reconstruction and fall back to b'N'. This may indicate a reference/CRAM mismatch but can legitimately occur near chromosome ends. The warning MUST be emitted at most once per record to avoid flooding logs.

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.