refactor(bed)!: stop pruning contigs and reference sequence in intersect_contigs (#246)

[TimD1-bot]

Note

Authorship: the content below was drafted by Claude Opus 5 (an AI coding agent) and
filed via gh under @TimD1-bot, a bot account operated by @TimD1. It reflects the
agent's analysis, not a statement authored by @TimD1.

Addresses #246. Part of #48, and a prerequisite for retaining variants on contigs absent from
the -b BED file (#107).

What was wrong

intersect_contigs() erased four separate things:

Erased	Where
query contigs absent from the BED, with their lengths, observed_ploidies and variants[hap]	src/bed.cpp
the same for truth	src/bed.cpp
reference FASTA sequences absent from the BED	src/bed.cpp
reference sequences absent from the truth VCF (no-BED branch)	src/bed.cpp

A variant on an erased contig could never be reported: the contig was gone from the output
contig list before the writer ran, and its reference sequence had been freed too —
ctgVariants::set_var_record reads ref->fasta to anchor an INS or DEL, so retaining the
contig alone would not have been enough. All four loops are deleted.

Evaluation is unchanged

parse_variants() already discards every variant on a BED-absent contig as BED_OFFCTG
(src/variant.cpp:1177), so a retained contig arrives empty, and clustering, superclustering
and phasing all skip empty contigs. On the committed chr20 fixture,
precision-recall-summary.tsv and phasing-summary.tsv are byte-identical before and after,
as are precision-recall.tsv, query.tsv, truth.tsv, phase-blocks.tsv, switchflips.tsv
and genotype-errors.tsv. The only output difference anywhere is a ##contig header line for
the retained contig in summary.vcf, which carries no record — that header is what #107 will
write into.

Warnings are now scoped to BED contigs

The cross-callset false-positive/false-negative warnings and the observed-ploidy comparison
both used to run after the erase, so they only ever saw evaluated contigs. Left global they
would fire for contigs nothing evaluates, where All query variants on 'chrN' will be false positives is simply untrue. Ploidy is likewise recorded before the BED filter applies
(src/variant.cpp:1025 precedes :1177), so a BED-absent contig carries ploidies no
evaluation consults.

This deviates from #246, which proposed leaving the ploidy comparison global. The same argument
the issue makes for the cross-callset warnings applies to it unchanged.

Two memory defects this would otherwise have exposed

Both stem from the same broken invariant. Pruning plus the injection block used to guarantee the
two callsets carried identical contig lists; without pruning, a one-sided contig outside the BED
survives unpaired, since injection iterates g.bed.contigs alone.

An out-of-bounds read. The ploidy comparison indexed query->observed_ploidies with the
result of a std::find that could fail, yielding size(). Verified with AddressSanitizer —
before the fix the truth-only scenario aborts:

ERROR: AddressSanitizer: heap-buffer-overflow
    #0 intersect_contigs(...) bed.cpp:601
    #1 IntersectContigs_TruthOnlyContigOutsideBed...::TestBody() test_bed.cpp:1162

A plain build reads past the end with no visible symptom, so this had to land with the deletion.

A null dereference (d3f456f). superclusterData builds its contig list as the union of both
callsets (src/cluster.cpp:319-338), then load_and_merge_callset_vars_across_haps indexes
vars[h][ctg] for every contig in it (src/cluster.cpp:97). vars[h] is an unordered_map, so
operator[] on an absent contig inserts a null shared_ptr and ->n dereferences it.

The trigger is narrower than "present in one callset, absent from the other": parse_variants
creates a ctgVariants per contig in the VCF header (src/variant.cpp:869-874) but appends
to contigs only per record (src/variant.cpp:924), so a declared-but-uncalled contig is still
indexable. It has to be absent from the other callset's header — which two callsets from
different pipelines routinely are. On such a fixture the pre-change build exits 0 and pruning-
removed-only exits 139.

tests/unit/src/test_cluster.cpp already documented this hazard on declare_contig(), noting a
contig one callset never declared "is not a state these tests construct" — because
intersect_contigs used to rule it out.

Restoring the invariant, not guarding each consumer

The no-BED branch's cross-injection moves into pair_missing_contigs(), called unconditionally in
both directions after the g.bed_exists branch; those two loops keep only their WARN calls.
ref->lengths.at() stays safe because the reference-FASTA assertion has already rejected any VCF
contig the FASTA lacks.

A paired contig is empty on both sides, so it is skipped everywhere and no metric moves. The
continue guard in the ploidy loop is kept as a bounds guard, now unreachable because pairing
guarantees the lookup succeeds rather than because the BED skip does.

Memory

No peak regression: fastaData loads the entire FASTA in its constructor (src/fasta.h:31-40),
before intersect_contigs runs, so the erase was a post-hoc free rather than a prevention.
Residency during alignment rises by the size of the kept contigs — roughly 3 GB for hg38,
against the 64 GB --max-ram default.

Breaking change

A contig present in either input VCF but absent from the reference FASTA is now a hard error
naming the contig and the callset:

Contig 'chrEBV' found in QUERY VCF but not reference FASTA.

It was previously dropped silently whenever it was also absent from the BED, so a query VCF
carrying decoys or alts (chrEBV, HLA-*) against a chromosome-scoped BED ran and reported
nothing about them. Retaining the contig means its reference sequence is now required, so the
omission has to be reported. Documented in README.md, beside the existing coordinate-sorting
requirement.

The old no-BED message ended Please provide BED file.; that advice is dropped, because
supplying a BED no longer avoids the error. In that branch this replaces a worse failure rather
than introducing one — a query-only contig missing from the FASTA already reached
ref_ptr->lengths.at(ctg) and threw an unhandled std::out_of_range.

Tests

Unit, 810 → 819. In test_bed.cpp:

• a contig outside the BED keeps its entry in contigs, lengths, observed_ploidies and
  variants[hap], and its reference sequence survives
• without a BED, the reference keeps contigs the truth VCF does not carry
• a query contig and a truth contig absent from the reference FASTA each error, naming the
  contig and the callset, with and without a BED
• the cross-callset warnings do not fire for a BED-absent contig
• a BED-absent contig with disagreeing ploidies produces no ploidy warning
• a one-sided BED-absent contig is paired into the other callset, in either direction

In test_cluster.cpp: a superclusterData built from an unpaired BED-absent contig — which
segfaults without pair_missing_contigs(). It calls intersect_contigs and the constructor
together, because neither function alone can show the invariant holds.

IntersectContigs.BedDropsExtraneous asserted the deleted behavior and is replaced by
BedKeepsContigOutsideBed; NobedFastaMissingErrors is renamed and its expected message
updated.

Integration (tests/integration/test-integration.yml), three cases. Two reuse the existing
one_sided_contig_* VCF pair against synthetic.bed (sc1 only) instead of synthetic_2ctg.bed,
so sc2 falls outside the evaluated regions in each direction; each pins that sc2 reaches the
output contig list but contributes no record and no metric. The query-only case contrasts directly
with test_one-sided-contig_query-only, which runs the identical VCFs against a BED covering sc2
and scores those same calls as false positives.

The third adds undeclared_contig_{query,truth}.vcf, where the truth's header omits sc2
altogether, and runs the segfault case end to end — a crash is worth catching in CI. Both new
fixtures are documented in tests/integration/data/README.md and needed git add -f, since a
broad data/ rule in .gitignore covers that directory, as the existing 30 fixtures did.

Verification

• tests/unit/build/test_vcfdist: 819 passed, from a clean build, no compiler warnings
• same suite under -fsanitize=address: 819 passed, no sanitizer reports
• pytest in tests/: 156 passed
• chr20 before/after: all eight metrics files byte-identical, including both acceptance files
• the header-mismatch fixture: exit 139 with the pruning removed alone, exit 0 and output
  identical to the pre-change build with the pairing fix

Every new test was confirmed red first. The two BED-scoped integration cases fail on the
pre-change build specifically on the contig-retention assertion, with every metric assertion
passing — which is what makes the "unchanged by construction" claim checkable rather than
asserted. The superclusterData test was red by segfault, not assertion failure.

Closes #246 is omitted deliberately: this PR targets dev while the default branch is
master, so a closing keyword would not fire. #246 needs closing by hand on merge.

[TimD1-bot]

Claude Opus 5 🤖: The null dereference described in #246's comment was real and
not addressed by the first commit — d3f456f fixes it. I had checked the wrong thing and concluded
it was already covered; recording the correction here since the reasoning matters for the fix.

Why my integration test missed it

test_bed-absent-contig_query-only-not-evaluated sets up what looks like the exact trigger — query
carries sc1+sc2, truth carries sc1, BED covers sc1 alone — and passes. That is because
parse_variants creates a ctgVariants for every contig in the VCF header
(src/variant.cpp:869-874) but appends to contigs only once a record is seen
(src/variant.cpp:924). Both one_sided_contig_* fixtures declare sc2 in their headers, so
truth->variants[hap]["sc2"] exists and the merge finds it.

So the trigger is narrower than "one callset carries a contig the other does not": the contig has to
be absent from the other callset's header. That is not an exotic input — two callsets from
different pipelines routinely disagree on header contigs.

Confirmed with a fixture whose truth header omits sc2:

pre-change build (pruning intact):  exit=0
first commit  (pruning removed):    exit=139   # SIGSEGV during superclustering
with d3f456f:                       exit=0, output identical to the pre-change build

Worth noting that tests/unit/src/test_cluster.cpp already documented this hazard on
declare_contig(), stating that a contig one callset never declared "is not a state these tests
construct" — precisely because intersect_contigs used to rule it out.

The fix

Taken as recommended: restore the invariant centrally rather than guard each consumer. The no-BED
branch's cross-injection moves into pair_missing_contigs(), called unconditionally in both
directions after the g.bed_exists branch; those two loops keep only their WARN calls.
ref->lengths.at() stays safe because the reference-FASTA assertion has already rejected any VCF
contig the FASTA lacks.

One deviation: the continue guard in the ploidy loop is kept, but its comment now says it is
unreachable because pairing guarantees the lookup succeeds — not because the BED skip does. Both
reasons are now true; the pairing one is the stronger.

Verification

• unit: 819 pass (817 → 819), clean build, no warnings
• same suite under -fsanitize=address: 819 pass, no sanitizer reports
• pytest: 156 pass (151 → 156)
• chr20: all eight metrics files still byte-identical to the pre-change baseline

Each new test was confirmed red first — the superclusterData one by segfault, not assertion
failure.