`extract` rejects indexed all-unmapped inputs with no reference targets

[SuhasSrinivasan]

Summary

Indexed extract full/calls rejects valid all-unmapped BAM/CRAM files whose header has zero @SQ reference targets. Serial extraction processes the same ten records, but indexed scheduling fails before reaching the command's existing unmapped fallback with should be at least 1 contig.

The empty mapped-coordinate schedule is a valid terminal state for extract because extract already owns a separate unmapped pass. It should not be represented by a synthetic contig or treated as an input error.

Severity

Severity: Medium — scientific completeness and format/mode consistency

Rationale: The failure is explicit rather than silent, but valid records become unusable solely because an index is present. BAM/CRAM and serial/indexed execution modes should not change whether an all-unmapped dataset can be extracted.

User and scientific impact

• Affected workflows are indexed extract full and extract calls on all-unmapped inputs without reference targets.
• Explicit-threshold, no-filter, and automatic modes are affected before their normal unmapped behavior can run.
• Adding or removing a sidecar index changes a valid command from complete output to failure.
• This is relevant to unmapped direct-RNA reads carrying modification tags and to intermediate modBAM/CRAM subsets with no coordinate dictionary.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Reproduced on macOS arm64 with ten-read all-unmapped BAM/CRAM fixtures and BAI/CRAI.

Steps to reproduce

From the repository's modkit package directory, with samtools available:

tmp=$(mktemp -d /tmp/modkit-zero-target.XXXXXX)

cp ../tests/resources/bc_anchored_10_reads.unmapped.bam "$tmp/serial.bam"
cp ../tests/resources/bc_anchored_10_reads.unmapped.bam "$tmp/indexed.bam"
samtools index "$tmp/indexed.bam"

cp ../tests/resources/bc_anchored_10_reads_unmapped.cram "$tmp/serial.cram"
cp ../tests/resources/bc_anchored_10_reads_unmapped.cram "$tmp/indexed.cram"
cp ../tests/resources/bc_anchored_10_reads_unmapped.cram.crai \
  "$tmp/indexed.cram.crai"

for format in bam cram; do
  modkit extract full "$tmp/serial.$format" "$tmp/serial.$format.tsv" \
    --threads 1 --io-threads 1 --no-headers --suppress-progress --force

  modkit extract full "$tmp/indexed.$format" "$tmp/indexed.$format.tsv" \
    --threads 1 --io-threads 1 --no-headers --suppress-progress --force
done

run_calls() {
  label=$1
  input=$2
  shift 2
  modkit extract calls "$input" "$tmp/$label.tsv" \
    --sample-num-reads 10 --threads 1 --io-threads 1 \
    --no-headers --suppress-progress --force "$@"
}

for format in bam cram; do
  run_calls "calls.auto.serial.$format" "$tmp/serial.$format"
  run_calls "calls.auto.indexed.$format" "$tmp/indexed.$format"
  run_calls "calls.explicit.serial.$format" "$tmp/serial.$format" \
    --filter-threshold 0
  run_calls "calls.explicit.indexed.$format" "$tmp/indexed.$format" \
    --filter-threshold 0
  run_calls "calls.no-filter.serial.$format" "$tmp/serial.$format" \
    --no-filtering
  run_calls "calls.no-filter.indexed.$format" "$tmp/indexed.$format" \
    --no-filtering
done

wc -l "$tmp"/*.tsv

Control or independent oracle

The serial and indexed files contain the same ten unmapped records. extract full has an existing unmapped pass and should emit each record exactly once. The serial control emits 218 rows. extract calls --filter-threshold 0 similarly emits 109 rows.

Observed behavior

On installed modkit 0.6.4, serial BAM and CRAM each succeed with ten reads and 218 full rows. Indexed BAM and CRAM both exit 1 before output with:

Error! should be at least 1 contig

The same failure occurs in calls automatic, explicit-threshold, and no-filter modes. Removing the index makes the valid input succeed.

Expected behavior

For extract consumers that already process unmapped records separately, zero raw reference targets should produce a genuine empty mapped terminal schedule followed by exactly one unmapped pass. BAM/CRAM and serial/indexed modes should emit the same rows and modification fields. Consumers whose contract requires mapped coordinates should retain strict rejection.

Root-cause evidence

The generic interval feeder and automatic-threshold schedule constructors require at least one contig. Extract invokes them before its indexed-unmapped fallback. Selector optimization can also erase the distinction between a genuinely target-free header and a nonempty header whose targets were filtered out, so raw target provenance must be recorded before optimization.

Proposed fix scope

• Keep strict feeder constructors as the default.
• Add an explicit opt-in empty-terminal constructor/state only for consumers that already own unmapped processing.
• Preserve whether the raw header had zero targets before region/BED/motif optimization.
• Let extract full/calls and the automatic-threshold caller select this terminal mode, then run the existing unmapped pass once.
• Complete CRAM setup before spawning mapped workers so an empty mapped schedule does not bypass required reader initialization.

Non-goals

• Do not invent a synthetic reference target.
• Do not relax Summary, pileup, hemi, bedMethyl, or other coordinate-required consumers.
• Do not change fixed-count estimands, coordinate-selector exclusion, or coordinate-placed FLAG-4 routing.
• Do not conflate a zero-target header with a nonempty header filtered to zero selected coordinates.

Acceptance criteria

• Indexed and serial BAM/CRAM all-unmapped inputs emit ten unique records exactly once.
• Full output has 218 rows and calls output has 109 rows in explicit and no-filter modes.
• Automatic calls mode samples the unmapped probability population; a genuinely empty probability population follows the separate explicit empty-sample error contract.
• BAM/CRAM row maps, thresholds, flags, coordinates, and no-header byte output agree.
• Region, motif, include-BED, and mapped-only controls do not accidentally include unmapped records.
• Strict Summary and coordinate-only feeder controls continue to reject zero targets.
• Input files and indexes remain unchanged.

Reproduction artifacts

Artifact	SHA-256	Notes
tests/resources/bc_anchored_10_reads.unmapped.bam	e1a802203410fccc9f6b4258cef3d9ca2cf93bb3dfc487f63e2acfa57678e144	Ten all-unmapped records, no reference targets
Generated BAI	a12ca6e7d51835746d5964310ca3559613fc71bf694be0821ca0c80c7ca3a603	Hash from the frozen reproduction environment
tests/resources/bc_anchored_10_reads_unmapped.cram	88642339b210a7627b88aec762054af8cd545b8ac848b3263a8b54f3da803cb5	Equivalent all-unmapped CRAM
tests/resources/bc_anchored_10_reads_unmapped.cram.crai	0f55839b846f8a33dd2ba101d357d231a0c987aea2489d1f949980fa1c2dc23a	CRAI

Related work

• CRAM reference-reader plumbing is a prerequisite for reference-dependent mapped inputs but does not fix a zero-target schedule.
• Empty automatic-threshold semantics define the no-observation case reached after scheduling succeeds.
• Issue Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654 concerns fixed-count selection, not explicit/no-filter zero-target routing.
• True-mapped preflight and coordinate-placed unmapped records are separate predicates and should not be absorbed into this repair.
• Issue #189 reported the same should be at least 1 contig diagnostic for a reference/header mismatch. Here the header genuinely has zero @SQ targets and valid unmapped records; serial extract succeeds, so the repair is a consumer-scoped empty mapped schedule rather than reference/header correction.