CRAM references are not propagated to all threshold and extract readers

[SuhasSrinivasan]

Summary

Several commands accept or require a reference FASTA for CRAM input but do not attach it to every reader that decodes reference-dependent CRAM records. Legacy threshold readers cannot receive a reference at all in some command paths, and extract full/calls omits the supplied reference from its main serial and indexed body readers.

Equivalent BAM and CRAM inputs can therefore produce different results. In the most serious extract case, modkit 0.6.4 exits successfully after downgrading the CRAM body decode failure and writes zero rows instead of the expected result.

Severity

Severity: High — silent format-dependent data loss

Rationale: A valid reference-dependent CRAM can yield an empty successful result while the equivalent BAM emits hundreds of rows. Other affected paths fail with a misleading zero-contig diagnostic even though the command already has a reference option or should expose one.

User and scientific impact

• Affected paths include reference-capable legacy automatic-threshold consumers and extract full/calls entry, indexed mapped, indexed-unmapped, fixed-count schedule, and truly unindexed automatic-threshold readers.
• BAM/CRAM conversion alone can change output from complete to empty.
• Missing or wrong references may be reported late or unclearly, and outputs can be created before the true input problem is visible.
• Reference-independent all-unmapped CRAM must remain usable without forcing a FASTA.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Reproduced on macOS arm64 with equivalent checked-in BAM/CRAM plus FASTA and BAI/CRAI.

Steps to reproduce

From the repository's modkit package directory:

tmp=$(mktemp -d /tmp/modkit-cram-reference.XXXXXX)
bam=../tests/resources/bc_anchored_10_reads.sorted.bam
cram=../tests/resources/bc_anchored_10_reads.sorted.cram
ref=../tests/resources/CGI_ladder_3.6kb_ref.fa

modkit extract full "$bam" "$tmp/full.bam.tsv" --reference "$ref" \
  --threads 1 --io-threads 1 --no-headers --suppress-progress --force
modkit extract full "$cram" "$tmp/full.cram.tsv" --reference "$ref" \
  --threads 1 --io-threads 1 --no-headers --suppress-progress --force

modkit extract calls "$bam" "$tmp/calls.bam.tsv" --reference "$ref" \
  --filter-threshold 0 --threads 1 --io-threads 1 --no-headers \
  --suppress-progress --force
modkit extract calls "$cram" "$tmp/calls.cram.tsv" --reference "$ref" \
  --filter-threshold 0 --threads 1 --io-threads 1 --no-headers \
  --suppress-progress --force

wc -l "$tmp"/*.tsv

The legacy threshold-reader gap can be seen with call-mods:

modkit call-mods "$cram" "$tmp/call.sam" --output-sam \
  --sampling-frac 1 --seed 7 --threads 1 \
  --sampling-interval-size 20 --suppress-progress

Control or independent oracle

The BAM and CRAM encode the same ten alignments. After normalizing format-specific headers, extract full should emit 218 identical rows and extract calls --filter-threshold 0 should emit 109 identical rows. The supplied FASTA is the reference used to create the CRAM.

Observed behavior

On installed modkit 0.6.4:

Command	BAM status/rows	CRAM status/rows
extract full	0 / 218	0 / 0
extract calls --filter-threshold 0	0 / 109	0 / 0

The CRAM body logs one failed record and returns success. call-mods does not expose a reference argument in 0.6.4; without one, the same CRAM fails during sampling with should be at least 1 contig, while the BAM succeeds with ten records.

Expected behavior

Every reader that may decode reference-based CRAM records must receive the user-supplied FASTA. Equivalent BAM and CRAM must produce identical scientific row/call maps. Missing or wrong references detected by bounded preflight must fail clearly before publishing output; later decoder failures that pass that preflight remain separate lifecycle work. Reference-free header/index probes and genuinely reference-independent CRAM should retain their current behavior.

Root-cause evidence

The reference path is present in command arguments but is dropped at internal reader-construction boundaries. Some legacy sampling APIs have no reference parameter. In extract, initial serial readers, indexed interval readers, indexed-unmapped readers, and schedule constructors create new readers without calling the CRAM reference setter. Header/index probes do not decode record bodies and should remain reference-free.

Proposed fix scope

• Add reference-aware internal sampling/schedule entry points while retaining compatibility wrappers where needed.
• Pass optional references from call-mods/adjust-mods and existing required references from duplex/entropy consumers to every threshold body reader.
• Thread an owned optional reference path through extract full/calls, including serial, mapped indexed, indexed-unmapped, fixed-count schedule, and truly unindexed automatic-threshold readers.
• Detect CRAM by content rather than filename extension.
• Validate a missing/wrong reference before output mutation when bounded preflight is possible.

Non-goals

• Do not redesign fixed-count selection semantics; issue Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654 owns that policy.
• Do not change --ignore-index routing, zero-target scheduling, or mapped/unmapped eligibility.
• General late decoder-error propagation and transactional rollback remain separate lifecycle work; this issue fixes decode failures caused by dropping the supplied reference.

Acceptance criteria

• BAM/CRAM parity holds for extract full (218 rows) and extract calls (109 rows) in automatic, explicit, and no-filter modes.
• Indexed, truly unindexed, stdin, mixed mapped/unmapped, and fixed-count reader constructors are covered.
• call-mods, sampled adjust-mods, pileup-hemi, and entropy threshold readers use the reference where supported.
• Renaming a CRAM with a .bam suffix does not bypass reference setup.
• Missing/wrong-reference fixtures detected by bounded preflight return nonzero and preserve an absent output or existing sentinel; later decoder failures remain separate lifecycle work.
• Reference-independent all-unmapped CRAM remains usable without a FASTA.

Reproduction artifacts

Artifact	SHA-256	Notes
tests/resources/bc_anchored_10_reads.sorted.bam	4441acbe1ad59caf6ab5d56bf18097d8999c08351395219c919cb863652d85f0	BAM control
tests/resources/bc_anchored_10_reads.sorted.cram	ec78f57b6524d66ac9e82e0dd1f4fc24f21a9f5959388ffe2a3d2874bcdf881a	Equivalent reference-dependent CRAM
tests/resources/bc_anchored_10_reads.sorted.cram.crai	783940ec8f285dba850852ba3796fb1ecd27991a2c1b930f78c3fa3ea7b0dad1	CRAI
tests/resources/CGI_ladder_3.6kb_ref.fa	7c02f1e6dc3a3dae186ae47696c822efafe849b8c8526bc131237bf7c3a64570	Correct reference

Related work

• Issue Probability sampling can hang or return partial histograms after reader and worker failures #677 concerns generic fatal probability-worker/reader errors; this issue is specifically reference propagation.
• Issue Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654 concerns fixed-count selection and is not a CRAM decode fix.
• Indexed all-unmapped inputs with zero reference targets require a separate empty-schedule repair even after reference plumbing is correct.
• Issue #584 discusses CRAM reference paths in pileup/container workflows. The reproductions here use local equivalent BAM/CRAM inputs and cover additional legacy-threshold and extract body-reader construction boundaries, including successful empty extract output.