CRAM region queries decoded the whole chromosome; skip non-overlapping containers

[JamesKane]

Every CRAM region query cost the same as reading the entire contig, no matter how small the region.

Measured on a 30× 1kGP CRAM, a one-base query:

contig	length	query cost	same query on a BAM
chr21	45 Mb	20.92 s	4–6 ms
chr1	248 Mb	116.4 s	4–6 ms

Cost tracked contig length exactly (248/45 = 5.5× length, 116/21 = 5.5× time) and was completely independent of region size — a 1 bp and a 1 Mb query on chr21 both took 20.91 s.

Cause

noodles' Query::read_next_container skips a container only when its reference_sequence_id doesn't match — it never consults the requested interval. So it decodes every container of the chromosome and filters records afterwards via intersects(&record, self.interval). Our own for_each had replicated the same reference-sequence-only filter.

The .crai already carries alignment_start and alignment_span per container. The information needed to skip was present and unused. On HG00096, chr21 holds 1,140 containers and a point query needs exactly 1.

Fix

cram_container_offsets selects containers by interval overlap, shared by both read paths. query no longer delegates to noodles — it decodes only containers that can overlap, lazily, one at a time, so an early-stopping caller doesn't pay for the rest.

A container the index cannot place (alignment_start = None) is kept. Skipping is only ever done on positive evidence that a container lies outside the interval, because a wrongly skipped container is reads silently missing from a variant call.

Result

	before	after
chr21, 1 bp	20.92 s	8.3 ms	~2,500×
chr1, 1 bp	116.4 s	113 ms	~1,030×
chr21, 1 Mb	20.91 s	438 ms	~48×
navigator call --contig chr21	>96 min (killed, unfinished)	44 s	~130×

The last row is the one that matters: the whole-chromosome cost was paid once per caller chunk and once per realignment candidate, which is why a single chromosome ran 96 minutes without finishing.

Correctness

Faster-but-lossy would present as a faster caller rather than a broken one, so equivalence is tested rather than assumed:

• cram_query_matches_noodles_query — runs our query and noodles' own Query over the checked-in fixture across three regions, comparing every field of every record.
• container_offsets_select_only_overlapping_containers — pins the selection boundaries: first/last base of a container overlaps, one past either end does not, unbounded intervals keep everything, other references are never selected, unplaceable containers are kept. The fixture is a single container and cannot catch an off-by-one here.
• On the real 11 GB CRAM — VERIFY=1 cram_query_probe compares against noodles record for record: 305 records over 1 kb and 47,898 over 200 kb, byte-identical on name, position, flags and sequence.

cargo clippy --all-targets -- -D warnings clean; 30 test suites pass.

Notes

Ruled out before finding the cause, in case they come up again: not the NAS (204 MB/s, and the process sat at 98% CPU — a local copy behaves identically), not the reference FASTA (45 MB chr21-only reference: 20.93 s vs 21.00 s with the full 3.1 GB), not CRAM 3.1 decode recursion (these are 3.0), not chunk count.

Why it went unnoticed: the ground-truth subject's alignment is a BAM, whose whole genome calls in ~5 minutes. Every CRAM in the workspace — including all 3,216 1kGP alignments — has been paying the whole-contig cost.

Found while trying to validate archaic segment calls against the hmmix 1000G callset, which needs those CRAMs. This unblocks that work: chr21+22 per sample goes from ~3.5 hours to ~90 seconds.

Adds cram_query_probe, the harness these numbers come from, as the regression tool — it attributes cost per phase (open, first query, warm query, bulk region) so a future slowdown names the part responsible.

🤖 Generated with Claude Code