feat(variant): retain source records and carry their columns via htslib (#248)

[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 #248. Part of #48. Supersedes the string-based design of #104 and #105, which #239 made
obsolete by switching the writer to htslib.

Based on dev rather than stacked on #247: that PR's only src/ change is bed.cpp, this one
touches variant.{h,cpp}, phase.cpp, and cluster.cpp, so there is no overlap — and a
dev-based branch gets CI. The #246 dependency belongs to the next PR in the chain, which retains
variants vcfdist does not evaluate.

What was wrong

summary.vcf was synthesized from the internal ctgVariants arrays rather than copied from the
input, so its site-level columns were hardcoded placeholders in ctgVariants::set_var_record
(src/variant.cpp:544): ID=., QUAL=., FILTER=PASS, INFO=.. Only CHROM/POS/REF/ALT carried
data, and those came from normalized alleles.

Change

Each variant now holds a std::shared_ptr<bcf1_t> for the record it was parsed from, and htslib
copies the columns onto the output.

• parse_variants bcf_dups each record once, on first use, and shares it across every variant
  derived from it — the two halves of a split complex variant, the two co-located entries of a
  het-alt record. Duplicating lazily means a record every haplotype discards is never copied.
• Shared provenance is now pointer equality, so rec_idx and rec_idxs are removed: they
  existed only as the key into a record-indexed store of reconstructed column text.
• The retained record is never mutated. Every output record is a fresh bcf_dup, which is what
  lets two entries subset one source to different alleles without aliasing.
• The parse header's ownership moves into a shared_ptr, since a retained record's tag IDs index its
  dictionaries and are meaningless without it. ctgVariants gains hdr and callset alongside
  recs; a second constructor takes all three so a creation site cannot forget them.
• summary_vcf_header folds both inputs in with bcf_hdr_merge, after building vcfdist's own
  contigs, FILTER, FORMAT, and samples — that ordering is what makes bcf_hdr_merge's
  destination-wins-on-collision behavior keep vcfdist's contig lengths and its own FORMAT
  cardinalities, both of which its records are written against.

Record construction is ordered around two htslib behaviors

bcf_remove_allele_set (htslib/vcfutils.h) does the allele subsetting, and per its contract
"Number=A,R,G INFO and FORMAT fields will be updated accordingly" — which is what removes any
hand-rolled subsetting. Header cardinality needs no rewriting either: the declarations arrive saying
Number=A/R/G and stay correct, since every output record is biallelic.

Two steps fail silently if reordered:

1. A biallelic placeholder genotype must be written before the removal. bcf_remove_allele_set
   fails on a genotype naming an allele it is removing, which is every het-alt record. The
   placeholder goes on through the record's own 1-sample header, not the output header, for
   reason 2 below; set_record_samples overwrites it with the real genotype afterwards.
2. Every carried FORMAT field must be read before any is written. bcf_update_format_* sets
   n_sample from the header it is given, so the first two-sample write flips the record; a read
   after that computes its value count as nsmpl * nvals against a buffer still holding one sample
   and returns the neighbouring field's bytes. This is the AD=17,7,97 corruption. Every read
   therefore happens against the 1-sample input header, before bcf_translate.

The re-keying reuses the existing pad_per_allele and update_format overloads, dispatching on
bcf_hdr_id2type for BCF_HT_INT, BCF_HT_REAL, and BCF_HT_STR. A FORMAT ID vcfdist writes
itself is not carried — its own declaration wins the header merge, so carrying the input's values
under it would contradict the declaration.

Metrics are unchanged

On the committed chr20 fixture, precision-recall-summary.tsv and phasing-summary.tsv are
byte-identical before and after, and the record count in summary.vcf is unchanged at 72,584.
bcftools view and bcftools stats read the output with no complaint, and leaks --atExit reports
zero leaks on both the chr20 and synthetic runs.

Verification

Each new test was checked to be load-bearing by mutating the implementation and confirming the right
test failed, since the tests for the writer were written alongside code rather than strictly before
it:

Mutation	Caught by
carried fields read against the 2-sample output header	3 tests; AD reports 97,17, floats become 4.36e+24
bcf_remove_allele_set never called	2 tests; AD=17,7,97 and un-subset AC/PL
carried values always written to the query column	CarriedValuesLandInTheOwningSampleColumn

The placeholder-genotype guard is the exception: it cannot be caught locally, because htslib 1.20
subsets a 1|2 genotype without complaint. On the 1.17 that CI pins, removal fails with "Problem
updating genotypes". The guard is kept and the test comment says so. htslib 1.17 has no osx-arm64
build on conda-forge and the release tarball is not reachable from this environment, so CI is the
only place that behavior is exercised.

Two things worth a look

Peak RSS on chr20 rose 142 → 223 MiB (+57%, ~508 bytes per retained input record). The issue
expected retention to be a small fraction of the alignment working set; at this rate a 5M-record WGS
callset would add roughly 2.5 GB. The WGS-scale measurement was explicitly deferred, so this chr20
figure is the only data point — but it is larger than predicted and may deserve action before this
lands.

A matched query/truth pair reports only the query's carried columns. One record is written per
variant, so when both callsets call the same variant the record is built from the query's source and
the truth sample's own ID, QUAL, FILTER, and INFO have no second record to go on; its carried FORMAT
fields report missing. This follows the design's "re-key from 1 sample to 2, non-owning sample
missing", but it is a real gap rather than an obvious consequence, so it is called out here.

Tests

• Unit, 861 total (+19). Retention: one record shared across haplotypes, across the two entries of
  a het-alt, and across both halves of a split complex variant; the retained record keeps every source
  allele. Header merge: input INFO/FORMAT/FILTER declarations survive, vcfdist's contig length and
  FORMAT declarations win a collision, samples stay TRUTH/QUERY, a null input header is skipped.
  Writer: site columns carried, unreported QUAL/FILTER carried as missing, Number=1/fixed/Flag
  carried verbatim, Number=A/R/G subset per allele with the G case asserted at alt_idx 2,
  both htslib ordering invariants, sample-column ownership, and the retained record unmutated by
  writing.
• Integration, +2. A new carried_fields fixture pair on the synthetic contig, because the chr20
  fixture cannot cover allele subsetting — every multiallelic record in it is unphased and so dropped
  at parse time. The second test reads the output back with bcftools, which is the only check that a
  carried field's value count still matches its declared cardinality; bcftools is added to the CI
  apt-get install for it.
• 8 existing integration assertions updated: their fixtures carry QUAL=50, which now reaches the
  output instead of ..

[TimD1]

rather than installing bcftools, can you use the htslib API to try reading in a gtest unit test?

[TimD1]

really? why not remove the additional comment text and use vars[callset][ctg]->hdr? That's much more clear

[TimD1]

please make these std::strings

[TimD1]

question: does this work correctly for complex variants that have been split?

[TimD1]

question: why is this a stand-alone function and not phaseblockData::input_headers()?

[TimD1]

nit: remove this comment

[TimD1]

shorten this, keep on one line like other comments

[TimD1]

ditto

[TimD1]

follow existing comment convention

[TimD1]

question: why was this test modified?

[TimD1]

question: are there existing functions (in htslib or this file) that could be used instead of writing new ones?