call: add opt-in read-level genotype likelihoods (--read-likelihood)

[benedictpaten]

What this adds

An alternative genotyping model for vg call that scores an explicit P(reads | genotype) rather
than fitting aggregate read depth to a Poisson model, plus a phasing layer over the GBWT haplotype
panel and a compact representation of the phased result.

It also carries long-read support: a --preset ont for the places the short-read defaults give the
wrong answer at 33 kb, phasing taken from the reads rather than the panel, re-genotyping from
that phase, and an --anchors-out file of read/allele partitions for pangenome-guided assembly.

On 43x ONT chr20 the long-read path moves ALL F1 from 0.926 to 0.952, indel GT F1 from 0.749 to
0.837, and switch error from 3.79% to 0.52% — each reproducing on a held-out contig.

It is opt-in behind --read-likelihood. Without that flag none of this code is constructed and
the default caller is untouched. That is a deliberate constraint, not an interim state.

Where the support-based callers ask "how many reads cover this allele", this asks "how well does
each individual read fit each allele", builds a reads × alleles matrix, and scores every allele
combination as a proper likelihood. Allele balance falls out of the 1/|G| mixture weights instead
of the het_bias knob, and multi-allelic sites need no special casing.

---

Trying it out

Reads in memory, small region:

vg call graph.gbz -p CHM13#0#chr20 -s HG002 -d 2 -t 8 \
    --read-likelihood --gam reads.gam \
    --mosaic-out out.mosaic.tsv > out.vcf

Whole contig, reads fetched per site instead of all held in memory. This is the one to use — no
extra dependency, and about twice as fast as the GAF-Base backend:

vg gamsort -i reads.sorted.gam.gai reads.gam > reads.sorted.gam

vg call graph.gbz -p CHM13#0#chr20 -s HG002 -d 2 -t 8 \
    --read-likelihood --gam reads.sorted.gam --gam-index reads.sorted.gam.gai \
    --mosaic-out out.mosaic.tsv > out.vcf

The same thing against a GAF-Base database, which wins on disk but needs gbz-base and gaf-base
on PATH:

gbz-base construct graph.gbz -o graph.gbz.db
gaf-base construct reads.gaf -r graph.gbz -o reads.gaf.db

vg call graph.gbz -p CHM13#0#chr20 -s HG002 -d 2 -t 8 \
    --read-likelihood --gaf-base reads.gaf.db --gbz-base graph.gbz.db \
    --mosaic-out out.mosaic.tsv > out.vcf

Long reads want one extra flag, and assembly anchors are a fourth output:

vg call graph.gbz -p CHM13#0#chr20 -s HG002 -d 2 -t 8 \
    --read-likelihood --preset ont \
    --gaf-base reads.gaf.db --gbz-base graph.gbz.db \
    --anchors-out out.anchors.tsv > out.vcf

Use a GAF the mapper produced, not one converted from a GAM: vg convert -G is not a lossless
round trip (see Testing), so a converted GAF gives slightly different calls. All three of the
short-read invocations above give byte-identical VCFs on the same reads.

No pack file is needed: on a GBZ carrying a haplotype panel this path enumerates alleles from the
panel, and nothing in the pipeline reads support. -d is the ploidy and --ploidy-bed varies it by
region (chrX wants -d 1 plus a PAR BED). vg call --help groups every flag that only acts under
--read-likelihood under one heading, and passing one without it is a hard error naming the
offenders rather than a silently ignored flag.

---

Outputs

VCF

Phased genotypes (0|1) with FORMAT/PS are on by default wherever the linkage layer runs;
--no-phased opts out. Phase comes from the panel rather than from reads spanning sites, so a phase
set is a whole chain rather than a read-length block. Records carry GL/GQ/GQI/GQN/GP.
--min-confidence X marks records below GQN X as FILTER=lowconf rather than dropping them.

The mosaic file (--mosaic-out)

The phased result written as what the sample is, in panel terms: one line per maximal run of one
strand on one panel haplotype. Only switch points are stored, so it is roughly two orders of
magnitude smaller than writing the paths out explicitly. --mosaic-out implies --phased.

#H  contig  strand  fragment  ref_start  ref_end  start_node  end_node  hap_index  haplotype  sites  gbwt_offset
H   chr20   0       0         22         533      229637742   229638340  4          recombination#19  16  19
H   chr20   0       0         945        2268     229638340   229641408  9          recombination#30  14  9

The properties a consumer relies on:

• (contig, strand, fragment) is a path identity. Consecutive rows of one fragment meet at the
  same node, so concatenating them gives one walk. On chr20 that is 7,144 rows in 3 fragments with
  every junction met; on chrX, 5.
• start_node/end_node are ORIENTED node ids, id * 2 + is_reverse, and are the
  authoritative anchors. ref_start/ref_end are advisory — a strand leaves the reference, and two
  segments can share a node while traversing it in opposite directions, which is not a walk.
• haplotype (sample#phase) is the portable name; hap_index is internal to the run and the
  file carries a #haplotype table mapping one to the other.
• hap_index of ref marks a stretch filled with the reference because no panel haplotype could
  be carried across it — walkable like any other row, and a deliberate approximation.
• A fragment ends rather than containing anything that is not a walk: an inversion boundary, or a row
  whose haplotype the graph does not carry. The file is self-describing — every rule above is
  repeated in its own #note header lines.

It can be read back into exact graph paths, and the test suite asserts that it is:

vg paths -x graph.gbz -A > panel.gaf
vg view -g graph.gbz > graph.gfa
python3 test/mosaic_to_path.py --mosaic out.mosaic.tsv --gaf panel.gaf --gfa graph.gfa

That expands every fragment and checks each step is a real edge. It passes on chr20 and chrX.
--no-mosaic-nested merges switches inside nested chains into the enclosing run;
--no-mosaic-patch-gaps and --mosaic-break-unexplained break the walk instead of bridging it.

Full format documentation, including the phasing model it comes from, is in
doc/read-likelihood-genotyping.md.

---

The model

Three terms, all defaulted on the strength of the accuracy harness below, and all documented in full
in src/allele_likelihood.hpp and src/linkage_model.hpp.

Per-read likelihoods. Scored from each read's existing alignment in the graph rather than by
fresh DP (src/allele_likelihood.*), and every genotype enumerated exhaustively — cheap once the
matrix exists, so there is no top_k pruning (src/read_likelihood_caller.*).

--depth-term W (default 0.1). A Poisson term on total read count beside the per-read
likelihoods: W * ln Poisson(N_eff; lambda_G), with N_eff = sum_r (1 - e_r) so the same MAPQ
discount applies to the count as to the reads. It exists because a pure per-read model has no
opinion about absent reads, which is exactly the evidence a homozygous deletion presents.

--linkage-weight W (default 2) and --linkage-prior F (default 5). A PanGenie-shaped
Li–Stephens layer: hidden states are ordered pairs of GBWT panel haplotypes, transitions are
Li–Stephens over reference distance, and the emission is this caller's own ln P(reads | G) rather
than PanGenie's k-mer model. Genotypes come from forward-backward posteriors summed over the states
implying each genotype; the phasing and the mosaic come from a constrained Viterbi over the same
chain. At weight 0 the transition is uniform and the caller is bit-for-bit unchanged.

It defaults on but only where it can work — it needs panel enumeration and --read-likelihood.
Where either is missing the default declines silently and an explicit --linkage-weight is a hard
error, because a user who asks for linkage and does not get it should be told. It declines with a
warning below two panel haplotypes and warns below eight: the weight was tuned on a 34-haplotype
panel and measures as roughly neutral on four.

--linkage-prior matters more than the transition weight. It is the exponent on the
allele-frequency prior the state space implies — several haplotype pairs can spell one genotype, and
summing over them weights genotypes by how many do. It shipped capped at 1 by a guard that read as an
optimisation and silently made every larger value behave as 1. Uncapped it peaks near 5 and inverts
past 8.

On records the linkage pass changes, GQ becomes the phred complement of the posterior, discounted
by the explained-read share and capped at GQI. The cap is not a consistency nicety: the posterior
is computed under a strong prior, so 1 - posterior understates uncertainty exactly where the
per-site evidence was weakest, which is the only place the layer acts. Measured over four datasets it
is worth about +0.003 AUC and 1–2% fewer false calls at matched recall, against +0.0009 for the share
discount alone. It also makes GQ <= GQI hold on every record.

---

Accuracy

Five arms built from one vg binary, scored on both GIAB benchmarks — small variants through
aardvark, structural variants through truvari — across four datasets (chr20 and chr6 against a
4-haplotype and a 34-haplotype graph). One build per matrix, because a table whose rows come from
different builds is the failure mode the harness exists to prevent.

Small-variant genotype F1, current default against --read-likelihood with all terms on:

dataset	poisson (default)	readlik-z
chr20-4hap	0.9355	0.9507
chr20-34hap	0.9107	0.9645
chr6-4hap	0.9461	0.9602
chr6-34hap	0.9297	0.9689

SV F1 on the same runs: 0.4954→0.5016, 0.4535→0.4944, 0.5490→0.5691, 0.4944→0.5268. Better than the
current default on every dataset on both benchmarks, and best of all five arms on every
small-variant class on all four datasets.

chr20 phasing switch error is 3.03% (1,797 switches over 59,280 pairs), in one block spanning the
contig.

Two results kept visible rather than averaged away. readlik-z loses 0.0018 of SV F1 to readlik on
chr20-4hap — about 1.4 events on a 765-event truth set, below what that benchmark resolves. And the
linkage gain is roughly sevenfold larger on the 34-haplotype graphs than the 4-haplotype ones: both
linkage and panel multiplicity are panel-size effects, so these headline numbers belong to rich
panels.

The scope of the evidence, stated plainly. Every number here is HG002 against CHM13, chr20 and
chr6, two graphs — four datasets, one sample. A reviewer may reasonably want a whole-genome run or a
second sample before --linkage-weight ships on by default; it is the newest default and the one
whose benefit is most panel-dependent.

---

Defaults that changed, and where

Panel enumeration is the default under --read-likelihood, and only there. --enumerate-support
opts back out. The proposal was to make -z the default whenever a GBWT is present; measured, that
premise holds for one caller and fails for the other. Small-variant F1 improves in all eight
comparisons, but SV F1 under the Poisson caller falls on all four datasets (as much as
0.4535→0.4391). So the default changed here and nowhere else.

Two guards, both from the same mechanism — panel enumeration cannot spell an allele no haplotype
carries, so its ceiling is the panel's content. It counts haplotype-sense paths and declines
below two, because enumerating from a reference-only GBWT would offer the reference allele at every
site and lose alt recall in silence. A pack passed alongside the default now says it is unused rather
than being quietly ignored.

--read-window's GAF-Base default is 16384, up from 4096. This is the only change here that
moves existing short-read output, and only in floating point — see Performance and memory for the
measurement and for why a wider window is what long reads need.

New flags, all inert unless set: --gap-open / --gap-extend expose the read scorer's gap
penalties, which were unreachable (both scorer constructors already took them; call_main built
both scorers with all-default arguments). --preset ont sets those two plus --mismap-min. See
Long reads.

No pack file on this path. -k was only ever needed because FlowTraversalFinder is driven by
node and edge weights. Checked rather than assumed: FlowCaller's only use of the support finder is
a non-pure virtual returning a member, never overridden. The trap this exposed is worth a reviewer's
eye — SupportBasedSnarlCaller::get_skip_allele_fn() prunes any allele below a support threshold, so
against a zero-support finder it would prune every allele at every site; the caller overrides it
to never skip when support is unavailable, which only takes effect pack-free. Three cases are refused
rather than genotyped against zero support: no -g/-z, -v, and --bottom-up.

---

For reviewers

Scoring from the graph's implied alignment, not DP. In a variation graph the snarl decomposition
already is a multiple alignment: two traversals share the snarl's boundary nodes by construction.
DP could also find a higher-scoring alignment corresponding to no path in the graph, which is not
the quantity wanted — the alleles are paths, and what is needed is P(read | this path).
Counting differences also charges each mismatched base its own quality rather than a length average.
The cost is that it inherits the mapper's placement. Realignment was built twice as a throwaway and
neither pays: row normalisation divides out any improvement common across a site's alleles, and
optimal alignment changes which allele a read prefers 40 times in 91,914.

Rows are normalised by their own maximum, so every component shares a scale. The background
becomes exactly 1, "no valid placement" becomes 0 rather than -inf, the per-read term is
provably in [ln e_r, 0], and no logsumexp is needed anywhere.

The scoring window belongs to the read and the site, never to the allele. Bases an allele cannot
place are charged as insertions rather than omitted; omitting them would make the difference between
two alleles include match_reward × length_difference. There is an assert on this.

Two scoring bugs found and fixed — worth reading before reviewing the model. Reverse-strand reads
were compared against the wrong allele without being reverse-complemented (on the nested SNP fixture
this split perfectly by strand: 81 forward reads preferred the correct allele, 69 reverse ones the
wrong one — roughly half of all reads, at every site). And reads spanning a deletion were discarded
as uninformative, because informativeness asked whether the read touched a non-boundary node: all 150
deletion-carrying reads on the fixture were dropped, so the parent called hom-ref and the het
deletion vanished. The discriminating signal for a deletion is in the edge, not the node set.
Both were invisible to the tests that existed, which is the hole the concordance harness closes.

Two deliberate departures. GL is emitted in the VCF-specified genotype order, which
PoissonSupportSnarlCaller does not follow at 3+ alleles (flagged, not changed there). And the
posterior omits the -ln(candidate count) term, which cancels under a uniform prior but would make
GP vary with the number of alleles under exhaustive enumeration.

⚠️ One part is not opt-in

Nested calling needed a fix in FlowCaller that applies to PoissonSupportSnarlCaller identically.
When a parent allele does not traverse a child snarl, that child is reached by fewer haplotypes than
the parent's ploidy; the code asked for a full-ploidy genotype anyway and overwrote the positions
belonging to non-traversing alleles with star markers. Wrong twice: a site one haplotype reaches is
not diploid, and genotype() returns a sorted multiset with no haplotype identity, so which allele
got discarded was arbitrary. The child is now genotyped at the ploidy that traverses it. On a nested
site reached by one haplotype with 100 reference and 30 alternative reads, GQ goes from 38 to
256.

This touches graph_caller.cpp, which the default caller uses. It changed no existing test outcome,
but it can change --top-down output on other data, deliberately. Reviewers may prefer it split
into its own PR; happy to do that.

---

Long reads

The model was built and tuned on 30x 150 bp Illumina. Long reads turned out to need two things: a
run time that does not scale with read length where it should scale with overlap, and a different
answer to one question the defaults had only ever been asked about short reads.

Run time: 9.1x, and now cheaper per base than short reads

The first ONT chr20 run took 1806.9 s against 169.9 s for short reads on the same chromosome, on
data with 375 times fewer reads. That was not a long-read cost. Every per-read operation was
written as a walk over the whole alignment — which for a 150 bp read is the site overlap and for a
33 kb read is 200 times it — and there were five such walks per (read, site): touches,
get_read_steps twice, site_span_of, and each of the two anchor pins. Roughly 110 billion mapping
visits against 0.4 billion.

WindowedSiteReadSource now builds a per-window index — every mapping as (node, read, mapping),
sorted, with a read-offset table — once per fetch and reused by every site inside the window. Sites
find their reads by binary search instead of asking each read whether it touches, and the four
consumers take mapping indices and offsets straight from it, so each is O(site overlap). The index is
a hint: a source that cannot index cheaply leaves it null and consumers fall back to the walk.

chr20, -t 6, anchors on	wall	peak RSS
starting point	1806.9 s	4.5 GB
flip only the site's stretch of a read	459.5 s	6.8 GB
per-window mapping index	293.6 s	6.9 GB
--read-window 16384	236.4 s	7.5 GB
split queries run concurrently	198.7 s	6.6 GB
short reads, 30x, same options	158.6 s	8.7 GB

Those are 43.1x and 30.3x coverage, so 4.61 s per fold against 5.24 s: long reads now cost less
per unit of sequence than short ones. Every step is byte-identical on ONT.

What is left is the subprocess, not vg: the profile puts 76% of a run in the fetch, 56% blocked in
wait4. More concurrency does not help — -t 8 gives 203.7 s and -t 10 gives 212.2 s against
-t 6's 198.7 — so it is CPU-bound rather than latency-bound, and the remaining lever is a
gbz-base C API rather than anything in this PR.

--preset ont

--preset ont = --gap-open 1 --gap-extend 1 --mismap-min 0.05, applied after the option loop so an
explicit flag wins from either side of it.

Long-read accuracy had exactly one real gap, and it was not where the parameters were. SNV F1 is a
dead heat with short reads (0.98493 against 0.98517) and SVs are better (F1 0.5561 against 0.5318,
recall 0.631 against 0.578). The whole deficit is small-indel precision: 8,164 indel false
positives against 1,649, a factor of 5.

Stratifying aardvark's per-record BD by reference homopolymer run length localises it to one cell:

1 bp indels	ONT FP rate	short-read FP rate
HP 1–2	6.0% (91/1,507)	2.8% (42/1,497)
HP 3–4	8.4% (45/535)	2.5% (13/524)
HP >= 5	44.8% (5,174/11,540)	4.0% (286/7,230)

That cell carries 59.3% of every small-variant false positive at 11x the short-read rate, and
record counts outside it match between the two arms to within 2%.

The mechanism is that score_gap is the one scoring primitive QualAdjAlignmentScorer does not
override — its override list carries score_exact_match, score_mismatch,
score_full_length_bonus and score_partial_alignment, and the sole definition is
EditAlignmentScorer::score_gap. So an indel difference is scored quality-blind, and at the shipped
gap_open = 6 a 1 bp difference is worth 7 units, which at log_base = 1.383325 nats/unit is a
likelihood ratio of 6.2e-5 that no base quality can soften. Right for a read whose errors are
substitutions; wrong for one whose modal error is a single-base homopolymer indel. --mismap-min
damps the same per-read term, and swept, the two axes are almost exactly additive (+0.048 and +0.049
alone, +0.087 together).

chr20, 43x ONT	control	--preset ont
indel GT F1	0.74864	0.81568
SNV GT F1	0.98493	0.98466
ALL GT F1	0.92635	0.94477
indel FP / FN	8164 / 3975	5277 / 3055
SV F1 (precision)	0.5561 (0.497)	0.5641 (0.514)

Gated on criteria fixed before the data was seen: indel F1 gain >= 0.005 (got +0.0670), SNV cost
<= 0.001 (0.00027), reproduces on a held-out contig with the same sign at >= half magnitude
(chr6: +0.0610, 91%, and SNV F1 rises there), survives depth-matching (30.23x: +0.0692, the
same size as at 43x, so it is a read-length effect and not a coverage one), and no SV veto —
precision rises and DEL 1k+ het recall goes 0.303 to 0.333, the one class where ONT trailed.

Not a default: the same values cost 150 bp reads 0.0037 indel F1, because with no homopolymer
false-positive mountain to remove the recall loss dominates. The short-read default is unchanged to
five decimal places on this binary.

Two things a reviewer should know about its limits. Higher floors buy more — 0.10 gives +0.0866, 0.20
gives +0.1015 — at 0.0018 to 0.0096 SNV F1; they are excluded because they break the pre-registered
SNV budget, not because the trade is bad. And the preset is a global compensation for an error
whose rate varies 10–100x with local homopolymer length: dropping the whole HP>=5 cell is
catastrophic (it loses ~6,400 true indels to remove 2,825 false ones, F1 0.800 to 0.630), so the cell
is mostly real variation the caller gets wrong. An oracle removing only its false positives reaches
0.858 from the preset's 0.800, still 8x short reads' remaining headroom there. The targeted fix
is a context-aware gap penalty, which is not in this PR.

Two bugs only long reads could expose

FORMAT/DR was mis-scaled 3.25x. It reads 0.985 on short reads by design and read 0.308 on
ONT, tight and flat across size class. Two errors, both exactly zero at fixed read length, which is
why 151 bp never showed them: lambda = rate * (L + R - 1) needs rate in starts per bp and
local_read_rate counted reads overlapping the window (inflating it by 1 + R/W), and R needs
the population mean while the mean over the reads delivered to a site is size-biased — it estimates
E[L^2]/E[L], measured here as 78,165 bp against a population mean of 33,449. Counting only reads
that begin in the window fixes both with one mechanism. Median DR 0.308 -> 1.013, a 3.29x
correction against the 3.25x predicted; short reads unchanged to five decimals.

GQN was blanked after a linkage move, honestly — the pre-linkage value described the genotype
linkage moved away from — but that left --min-confidence, the caller's only filter, unable to see
the population it most needs to: those records run a 37.8% false-positive rate against 8.6%
overall, 4.4x enriched, holding 13.2% of every small-variant false positive. It is now re-derived
for the settled genotype as the same fraction, with the scale recovered from the record itself
(GQN was gap/achievable and GQI was that gap in phred, so their ratio is the denominator). It is
negative exactly when linkage overrode the reads, so the header now documents [-1,1].
GQN='.' 3,635 -> 438.

Neither moves a scored metric, which is the reason to gate a quality field on a counter.

Phasing from the reads (--read-phasing)

The shipped phase comes from the panel: the linkage layer's Viterbi path names a pair of panel
haplotypes and the order of that pair is the phase. The reads are never asked. At 33 kb that is
leaving the answer on the table — against a 349 bp median heterozygous spacing, 99.96% of
adjacent heterozygous pairs share a read.

Three stages, and adjacency is deliberately not the axis: chain consecutive reliable sites
stepping over the rest; relink each break over a few reliable sites either side; then hang the
unreliable sites off the settled chain, where a misplacement is a flip rather than a switch that
propagates. That shape is what the measurement forced — at the junctions stage 1 breaks, the
adjacent link is 15.06% wrong, the nearest reliable site either side is 0.97% wrong, and
the nearest site regardless of quality is 14.93%, indistinguishable from adjacent. Stepping over
bad sites is the entire effect; reaching further buys nothing.

switch error	panel	reads
chr20 ONT 44x	3.7942%	0.5180%
chr6, held out at the same values	3.1849%	0.3447%

Genotype-neutral by construction — a read's allele responsibility depends on the settled
genotype, not on the phase, so re-orienting relabels slots and nothing else. Verified byte-identical
on GT. Blocks stay chromosome-scale; only the orientation within them changes.

Re-genotyping from that phase (--regenotype)

Once the other sites are phased, the site's per-slot mixture weights can become per-read ones: a
read crossing 40 other heterozygous sites carries a strand log-odds, computed leave-one-out so the
site cannot vote on itself. Applied as a delta to the existing likelihood, so the depth term and
the normalisers cancel and nothing outside the mixture can drift.

	ALL F1	indel F1
chr20 off → on	0.94477 → 0.95151	0.81568 → 0.83725
chr6 held out, off → on	0.95342 → 0.95960	0.83781 → 0.86019

It aims at what is broken: the 6,350 chr20 sites it moves run a 40.6% false-positive rate
against a 6.7% background. A control that shuffles each read's strand sign while keeping |Λ| loses
0.107 F1, so the signal is the phase and not the magnitude. --regeno-temper 0 is exactly the
identity, byte for byte, which is the regression gate.

It runs as an EM over the barrier. --regeno-passes defaults to 2; beyond that the iteration can
enter a limit cycle — chr20 reaches a period-6 cycle at round 8 — which is detected and reported
rather than presented as a fixed point.

Two extensions, one kept and one not, both measured. A nested haploid chain has one slot and
nothing to reweight, so its reads are instead weighted by whether the phase places them on the
chain's strand at all; that is on by default and moves SV F1 (+0.0022 chr20, +0.0006 chr6) with
every other metric identical to five decimals. A two-parameter per-read escape calibrates the
agreement curve much better and makes the caller worse (ALL F1 −0.0014, indel −0.0047); it is
off by default and the reason is recorded in the code — the errors creating the ceiling are
correlated, being mismapping, which the site-level escape already answers.

Assembly anchors (--anchors-out)

For pangenome-guided assembly (PMC12621775):
at each called site, the reads that cross it partitioned by which chosen haplotype traversal they
match, written as zero-length pins — a point between two positions, so there is no anchor
sequence, no minimum length, and no way for members to disagree about what the anchor spells.

A pin is at one end of a boundary node, not at the node: two snarls sharing a boundary node in a
chain pin at opposite ends and never collide, whatever the node's length. Measured on chr20, all
529,304 pins land on 529,304 distinct (node, side) sites; drop the side and 37.3% collide. The
guarantee that no read position appears in two anchors is checked over the whole file and it found a
real defect while it was being built — the entry pin was stepping by read index, which invisibly
crosses a node the read deletes outright and lands on the neighbouring snarl's pin (11 collisions in
4.07 M placements, every one on a 1 bp node). It steps one mapping along the read's walk now.

Confidence is reported per site and per read, and read names are interned into a header table,
which is 32.8% of the file.

slot — which of the two called haplotypes a read is assigned to — took three format versions
to be right, and each bug parsed cleanly and joined wrongly. v2/v3 wrote it in allele order, so it
carried no phase at all. v4 fixed the diploid pair and left the haploid half standing: a nested
haploid chain sits on one strand of a diploid locus, and its single slot was stamped 0 whichever
strand that was, so every .|a site named the wrong haplotype — 3,886 chr20 sites, 2,334 of
which have no VCF line at all. The VCF and the mosaic are byte-identical across that fix, which
is why nothing else could see it. v6 adds a reliability column.

Which anchors to trust, since not all heterozygous anchors are equally good. Measured against
the GIAB truth over 59,919 labelled heterozygous chr20 sites, baseline FP rate 4.49%:

filter	kept	FP rate	FPs removed
reliability >= 9.5	85.3%	1.28%	75.7%
gqn >= 0.5	76.0%	0.73%	87.6%
both, at gqn >= 0.3	81.4%	0.78%	85.8%

Filtering barely costs linkage — with 43.1 heterozygous sites per read, singleton anchors stay at
0–1 at every threshold — but it costs block length, N50 43,774 → 20,891 sites, and the loss is
concentrated at about a dozen junctions around the centromere where the dropped sites were the only
bridge.

---

Performance and memory

The in-memory backend holds every alignment, which caps a run at a small region. --gam-index and
--gaf-base fetch on demand, sharing one windowing layer and differing in a single fetch_span()
primitive. Ordered visiting makes that pay: sorting top-level snarls by node ID and walking them
window by window means each window is fetched once — 920 site queries went from 0% to 98% cache hits,
17 window fetches instead of 920 index scans, 3.03 s → 0.42 s. That ordering is gated on an indexed
source, so the default caller's traversal order is bit-for-bit unchanged.

Measured on 80k reads over 400 kb, single-threaded, all backends producing the same 889 variants:

backend	wall	peak RSS	disk
in-memory GAM	0.73 s	309 MB	7.6 MB
indexed GAM (--gam-index)	0.92 s	56 MB	7.6 MB + index
GAF-Base (--gaf-base)	1.90 s	52 MB	2.8 MB

So --gam-index is the one that matters: it gets the memory reduction, is about twice as fast,
and adds no dependency. GAF-Base ties on memory and wins on disk but is an alternative rather than an
improvement; its dependency is runtime only — it runs gbz-base as a subprocess the way vg already
runs kmc, so nothing links. --read-window (default 256 for --gam-index, 16384 for
--gaf-base) is a genuine trade and deliberately not auto-sized: sparse sites want smaller
windows, so the obvious heuristic is backwards.

The GAF-Base default was 4096 until long reads showed why it should not be. A window has to be much
wider than a read's node-ID span or every read straddling a boundary is fetched twice; at 4096 that
was 1.62 fetches per ONT read. This is the one change in the PR that moves existing short-read
output, and it does so only in floating point: a window wider than max_query_nodes is fetched as
several concurrent queries whose results are concatenated, which reorders the reads a site sees and
so reorders a sum. On chr20 short reads that moves GL in the sixth decimal for 3,402 of 115,410
records and changes no genotype at all, verified record by record; --read-window 4096
reproduces the old output exactly. Gate a window change on genotype identity, not on cmp.

For a whole contig with the linkage layer on, chr20 runs in ~174 s and chrX in ~275 s at -t 5.

---

Testing

• t/18_vg_call.t: plan 405 (from 101 when this opened), on fixtures already in the repo.
  GAF-Base assertions self-skip with the reason in the test name when gbz-base is absent.
• vg test: 900 cases, 12.5 M assertions, including dedicated suites for the linkage model,
  the anchor writer, read phasing and re-genotyping.
• Byte-identity equivalence checks, which are what make a new read source reviewable: calls
  without a pack are byte-identical to calls with one (the direct evidence support is unused rather
  than read as zero); --gam-index matches in-memory; --gaf-base matches the same GAF in memory;
  ordered visiting matches unordered at -t 1 and -t 6; neither thread count nor --read-window
  changes a call. One caveat, precisely because it looks like a backend bug and is not: a GAM
  converted to GAF does not give identical calls, from two vg convert round-trip asymmetries
  (insertion re-partitioning across a node boundary, and a trailing zero-reference-base mapping being
  dropped). Raised separately; the GAF-Base test compares against in-memory GAF or it would be
  testing vg convert.
• Unit tests over hand-built inputs, so each layer is tested without a graph or GBWT in play:
  [allele_likelihood] 32 cases / 257 assertions, [linkage_model] 37 / 666, [site_read_source]
  13 / 39, [read_likelihood_caller] 6 / 38. The linkage suite found three bugs before the layer was
  ever wired up, and pins the two properties it must not lose: at zero weight the posterior is the
  per-site likelihood exactly, and an allele no panel haplotype carries stays callable.
• A mosaic structural fixture: a 396 bp graph with a parent snarl holding two children, an
  inversion over a stretch containing a site of its own, and a haplotype clipped into two GBWT
  fragments. The previous mosaic fixture gave one segment per thread, so no junction was tested at
  all. It immediately earned itself — it caught a containment test comparing node ids, which was
  suppressing 171 of chr20's 213 nested boundaries and passing only because real graphs are
  id-sorted.
• Thread invariance is asserted, not assumed: a windowed HMM parallel over windows is exactly the
  shape that introduces order dependence. It found a real defect — sites were sorted by position
  alone, so equal positions were ordered by whichever thread finished first.

---

Merged with master, and one semantic disagreement it surfaced

Merged origin/master (95 commits). Two conflicts, both substantive:

• Header emission. Master factored the six nesting INFO headers into
  VCFOutputCaller::nesting_info_headers(), now shared with Deconstructor. Took master's refactor.
• emit_variant's tail. Master added suppressed_ref_info: a site that produces no line still
  records its reference interval so its nested children have an RC/RS/RD to inherit,
  deliberately in step with Deconstructor::deconstruct_site. This branch's tail is the
  block-record emission, the wants_line decision, the linkage feed and the line-length warning.
  Both kept — master's recording goes in the else of wants_line, which is this caller's version
  of "produced no line".

One thing a reviewer should look at. Master's new test asserted that an island record whose
parent emitted nothing carries CH=0. This branch computes CH as the greater of the in-VCF
ancestor hops and the record's own gref contig level, because counting only ancestors that happened
to emit a record made a record on a gref fragment indistinguishable from one on the linear
reference — 29,843 of 41,669 off-reference records on a gref-covered chr20, which made the
documented bcftools view -i 'INFO/CH==1' filter select a quarter of what it should. The cost is
that CH >= 1 no longer implies an in-VCF parent, so it no longer implies PS; that coupling was
an accident of counting only emitted ancestors. The test now expects CH=1 and says why, and the
shared header description was updated to describe what the code computes. vg deconstruct is
unaffected — set_gref_levels is called only from call_main, so the floor is inert there and its
14 CH assertions are untouched. If you would rather keep master's definition, this is the one
change here to revert, and it is self-contained.

Linting, and four renamed options

make lint is a CI job of its own and gates the pipeline, and make test depends on it. It could
not run on this branch: scripts/lint.py raises on the first over-length helptext string and 95
had accumulated. Rewrapped by paragraph — never starting a wrapped line with --, which the checker
reads as an option declaration.

Four option names could not be expressed under the checker's column rule, which needs
flag + argument <= 19 characters for the description to land at column 28. All four are new here and
unreleased, so renaming costs no compatibility:

was	is
--anchors-min-read-score	--anchors-min-q
--anchors-end-pin-min-new	--anchors-end-new
--anchors-min-reads	--anchors-reads
--linkage-freq-prior	--linkage-prior

And seven options set a filename directly rather than through the standard checks; --gam,
--gaf-reads, --gam-index, --gaf-base, --gbz-base and --ploidy-bed now use
require_exists(), --dump-likelihoods uses ensure_writable(). That is a behaviour change worth
knowing: a missing file is now refused during option parsing, before any flag combination is
validated. Three tests were relying on the old order and now pass a real placeholder path.

Known limitations

• Reads are treated as independent. Mates overlapping a site are not, so confidence accumulates
  like R rather than √R and GQ is over-confident, increasingly so with depth. GL/GQ are
  useful for ranking, not as calibrated probabilities, and the VCF header says so. The fix (a
  scalar effective-sample-size discount) is additive and does not change the model's shape.
• Phase comes from the panel, not from reads. A phase set is a whole chain, so switch errors are
  panel-limited rather than read-limited.
• Two places the mosaic uses the reference, both stated approximations. A row whose named
  haplotype the graph does not carry across it is rewritten as a reference substitution (181 rows on
  chr20); and a boundary neither the earlier nor the later haplotype can be carried across is filled
  with the reference (36 on chr20, 0 left as a gap). The fill is the weaker of the two — 17 of those
  boundaries sit at over ten times average node density and the worst two at ~10,000× — so it is a
  contiguous path and substantively a guess. Each row says what it filled and how far; the judgement
  is the reader's. --no-mosaic-patch-gaps breaks the walk instead.
• Reference-crossed nested chains are only descended into where a parent candidate traversal reaches
  them, so a chain no candidate crosses is not called.

---

Status

The accuracy question that made this a draft is answered, so the open questions are about scope and
shape rather than whether the model earns its place.

Long-read support is newer than the rest and its scope is worth a separate opinion. --preset ont
is fitted on one sample — chr20 of HG002 at 43x, on a 16-haplotype graph, against a draft
benchmark whose own README flags errors in homopolymers and tandem repeats, which is the exact
epicentre of the finding — and validated on chr6 of the same read set and at matched 30x coverage.
It has never been checked against a second sample, which is already the largest unquantified
risk in the shipped short-read defaults and is now equally true here. A reviewer who would rather
not ship a fitted preset on that evidence has a fair point; --gap-open / --gap-extend stand on
their own and the preset can wait.

The same caveat applies with more force to the two flags the preset now also turns on.
--read-phasing is genotype-neutral and its gain is large and reproduces held out, so it is the
easier of the two. --regenotype changes genotypes, and while every arm measured is positive —
two contigs, two depths, precision and recall both up — it is the newest thing here and the one a
reviewer should push hardest on. Both can be turned off without touching the scorer values:
--preset ont --no-read-phasing --no-regenotype leaves exactly the three-parameter preset the
earlier revision of this description proposed.

SV is where this work is at best neutral, and the PR should not claim otherwise. The
re-genotyping's mixture correction costs SVs 0.0039 at 30x and the haploid weight gives back
0.0020; across chr20, chr6 and 30x the combined effect is +0.0047, +0.0003 and −0.0019. The gain is
at indels, which an earlier version of this plan had written off.

Two things a reviewer may reasonably want changed. The linkage layer's benefit is concentrated on
haplotype-sampled panels of tens of haplotypes and is close to inert on thin ones, so it may belong
behind a panel-size check rather than a bare flag. And most of its gain comes from the
allele-frequency prior rather than the transition model — a per-site prior reaches about 88% of the
genotype-F1 improvement with no forward-backward, no windowing and no retained state — so if the
Viterbi phasing pass that motivates keeping the chain is not wanted, the cheaper thing is defensible.
The phasing and the mosaic are what motivate keeping it.

Happy to be told the shape is wrong.

🤖 Generated with Claude Code

[faithokamoto]

Ah, you've been caught by one of the checks I added to force filename validation:

src/subcommand/call_main.cpp: Use require_exists() or ensure_writable() for standardized file checks: gam_filename = optarg;
src/subcommand/call_main.cpp: Use require_exists() or ensure_writable() for standardized file checks: gaf_filename = optarg;
src/subcommand/call_main.cpp: Use require_exists() or ensure_writable() for standardized file checks: dump_likelihoods_filename = optarg;

basically replace gam_filename = optarg; with gam_filename = require_exists(logger, optarg);

You can run this check by yourself with scripts/check_options.py