feat(variant): populate ac_errtype for the truth callset (#192) - #213
Merged
Conversation
ac_errtype was set for the query callset only, so every truth variant carried AC_UNKNOWN and the GE field on a TRUTH sample always rendered '.'. Record it on both callsets. The value keeps one absolute truth-allele-count-then-query-allele-count direction, so the two records of a matched site report it identically; what differs is which genotype supplies which count, since a record's orig_gt is its own callset's call and its calc_gt is the other callset's genotype recovered by alignment. This makes defs.h's existing per-value annotations reachable for the first time: AC_ERR_1_TO_0 and AC_ERR_2_TO_0, annotated TRUTH_FN, only occur on truth records, and AC_ERR_0_TO_1 and AC_ERR_0_TO_2, annotated QUERY_FP, only on query records. The GE header description, already phrased in truth -> query terms, becomes true of both samples unchanged. calc_prec_recall now mirrors onto matched truth variants what it already did for query variants: a TP records the matched query calls on the truth variant's calc_gt, unioned over the sync group, which is the atomic matching unit. set_allele_errtype takes the callset and reads one shared table from the corresponding side. A genotype carrying no diploid allele count now yields AC_UNKNOWN rather than being counted as zero alleles, which the old else-chain did silently for an unrecognized orig_gt. Production cannot reach it, since orig_gts is always heterozygous or homozygous alternate. No counts move: the hand-rolled truth false-negative tally is untouched and owns the genotype error summary. On the chr20 fixture the only output change is the GE subfield of TRUTH samples.
TimD1
added a commit
that referenced
this pull request
Aug 7, 2026
The rebase onto dev brought in set_allele_errtype(vi, query) and its allele_count/ac_errtype_from_counts helpers, plus a new per-haplotype loop in calc_prec_recall and a make_tvars test helper. Three signatures still spoke in raw integers: - allele_count takes gt_t, since orig_gts and calc_gts hold the enum - ac_errtype_from_counts returns ac_errtype_t - dist.cpp's new sync-group loop and make_tvars iterate hap_t and take gt_t dev's rewritten function body and its documentation are kept verbatim; only the types change. chr20 output is byte-identical to the dev tip, so #213's GE change survives intact.
TimD1
added a commit
that referenced
this pull request
Aug 7, 2026
…216) * refactor(defs): delete unreferenced constants and aliased spellings Removes TRUE, PTR_NONE and the five packed-pointer constants, which have no references anywhere; and TYPE_ALL, TYPE_INDEL, ERRTYPE_NE, REF and type_strs2, whose only references were assertions pinning the aliases themselves. TYPES and SWITCHTYPES are kept: the string tables are sized to them. FALSE was serving as both a zero initializer and a bucket index in pass_min_qual, so that is respelled directly. * feat(defs): add EnumArray and EnumRange for enum-keyed containers A scoped enum alone cannot catch table-index confusion, because std::vector accepts any integer subscript. EnumArray keys a fixed-size array to one enum family so a subscript from another family fails to compile; EnumRange iterates a family's contiguous enumerators, replacing size_t index loops over the tables. idx() converts an enumerator to its underlying integer where genuine index arithmetic is unavoidable. * refactor(defs): convert TIME_* to scoped enum timer_t Adds Globals::stage(timer_t) so the sixteen g.timers[TIME_*] subscripts are checked by type. g.timers stays a std::vector<timer> and timer_strs stays a std::vector<std::string>: init_timers accepts arbitrary-length input and appends rather than replaces, and four tests pin that behavior, so a fixed EnumArray would delete tested behavior rather than retype it. Progress banners print the stage ordinal, so those sites take idx() explicitly. * refactor(defs): convert MAT_* to scoped enum mat_t The three-dimensional offs and ptrs arrays become mat_t-keyed, and the loops over them iterate the enum. wf_swg_max_reach indexes a flat buffer as MAT_SUB*z + s2*y + d, which is genuine index arithmetic, so those sites take idx() explicitly. type_strs[mi+1] still borrows the variant type table with a matrix index; the static_assert and mat_to_edittype() helper land once edittype_t exists. * refactor(defs): convert BED_* to scoped enum bedloc_t bedData::contains now returns bedloc_t rather than int, region_strs and the nregions counter become bedloc_t-keyed, and locs holds the enum. The switch over loc drops its default: -Wswitch now reports an unhandled enumerator at build time, which is what the runtime ERROR() existed to catch. Verified by deleting a case and observing the warning. * refactor(defs): convert SWITCHTYPE_* to scoped enum switchtype_t switch_strs becomes switchtype_t-keyed, and the two locals in phase.cpp that accumulate a switch or flip classification hold the enum rather than int. * refactor(defs): convert PTR_* to scoped enum ptr_t Values stay non-contiguous at 1, 2, 4 and 8, so no EnumArray or EnumRange keys on this family. add_variants takes a std::vector<ptr_t> and print_wfa_ptrs takes ptr_t-valued matrices, which is the whole live surface: both functions are uncalled, recorded for follow-up. * refactor(defs): convert ERRTYPE_* to scoped enum errtype_t error_strs becomes errtype_t-keyed, as does the middle axis of the pr_counts query and truth counters, so a quality index or size class can no longer be passed where an error type belongs. errtypes and hap_fields::errtype hold the enum. * refactor(defs): convert AC_ERR_* to scoped enum ac_errtype_t Separates the sentinel from the count, which previously collided at 8: AC_UNKNOWN is an enumerator and AC_ERRTYPES stays a constexpr, with AC_ERRTYPE_SLOTS at 9 so the sentinel has a row of its own. allele_error_counts was sized to the count, leaving AC_UNKNOWN out of bounds and in range only because of an early ERROR(); it is now structurally in range. The genotype-errors TSV names all eight columns explicitly and never iterates rows, so the extra slot cannot reach output. * refactor(defs): convert PHASE_* to scoped enums phase_t and phaseptr_t Fixes the three-valued phase code stored in a bool: block_state holds phase_t, and the two XOR sites and the print_var_sample calls spell the intent as block_state == PHASE_SWAP rather than relying on a narrowing conversion. Only PHASE_ORIG and PHASE_SWAP reach those reads, so the result is unchanged. The phasing DP matrices are phase_t-keyed with phaseptr_t elements, phase ^= 1 becomes other_phase(), and the switch over phases drops its default now that -Wswitch covers it. PHASE_NONE keeps a slot of its own, separate from PHASES. * refactor(defs): convert VARTYPE_* to scoped enum sizeclass_t get_vartype returns the enum, vartype_strs is sizeclass_t-keyed, and both counter families are now typed on every axis: pr_counts is [sizeclass_t][errtype_t][qual] and allele_error_counts is [ac_errtype_t][sizeclass_t]. Neither axis can be subscripted with the other's key, nor with a bare quality index. * refactor(defs): convert GT_* to scoped enum gt_t orig_gts, calc_gts, var_fields::orig_gt, var_fields::calc_gt and the GT_counts histogram all hold the enum, and gt_strs is gt_t-keyed. Eleven enumerators make this the family -Wswitch covers most usefully. * refactor(defs): convert TYPE_* to scoped enum edittype_t types, var_fields::type, qtypes, ttypes and the ntypes histogram hold the enum, and bedData::contains takes it. Replaces the type_strs[mi+1] coupling with mat_to_edittype(), guarded by three static_asserts pinning each matrix one below the edit type it aligns, so drift breaks the build. With type_strs keyed, type_strs[MAT_SUB], error_strs[HAP1] and type_strs[2] are all now compile errors. The switch over variant type drops its default, covered by -Wswitch. * refactor(defs): convert QUERY/TRUTH to scoped enum callset_t callset_strs, callset_vars, samples, filenames and every callset-keyed local in the clustering and phasing paths become callset_t-keyed, as do the supercluster range and split signatures. Two changes go beyond retyping, both because a sentinel stopped being representable. The 'Invalid callset' range check and its two tests are removed: callset_t admits exactly two values, so the check can no longer fail. And var_info gains an explicit found flag, replacing callset_idx == -1; encoding not-found as a negative enumerator would have reintroduced the out-of-range subscript this refactor exists to remove. * fix(defs): close the switch gaps -Wswitch surfaced on the new enums Three switches silently fell through on an enumerator they never named, which -Wswitch reports now that the operands are typed. Each case is spelled out with the behavior it already had, so nothing changes at runtime: - dist.cpp's variant walk omitted TYPE_REF, which is not a stored variant type - print.cpp's INS and DEL matrix switches omit PTR_MAT, which those matrices never carry Two ERROR() calls printed an int8_t-backed enum through %d and now cast. * refactor(defs): convert HAP1/HAP2 and CTG_IDX/SC_IDX to scoped enums The last and widest family. Every per-haplotype container becomes hap_t-keyed: variants, errtypes, sync_group, callq, ref_ed, query_ed, credit, hap_fields, and the clustering locals. var_on_hap, set_var_calcgt_on_hap, print_var_sample, wf_swg_cluster and the graph constructors take hap_t. hap ^ 1 becomes other_hap(), and the phasing walk's qhi ^ swap ^ block ^ flip becomes an explicit 'flip if an odd number hold', which is what it always meant. As with callset, the hap > 1 range checks and their two tests are removed: hap_t admits exactly two values. A test that sized per-haplotype containers with PHASES rather than HAPS is corrected -- that is the cross-family confusion this issue set out to catch, and it was live in the suite. * refactor(globals): make the pipeline timers and their names enum-keyed timer_strs and Globals::timers become timer_t-keyed, closing the last family left as a plain vector. init_timers() drops its parameter and assigns each slot from timer_strs rather than appending. That removes a real footgun rather than only retyping one: the old init_timers appended, so a second call left twelve timers behind while stage(TIME_TOTAL) still read slot five. Assignment makes a second call idempotent, which is what the replacement test pins. The EmptyInput and AppendsNotClears tests go with the behavior they described, both now unrepresentable. timer.h no longer includes globals.h. It referenced nothing from it, and the cycle left 'timer' incomplete inside Globals, which std::vector tolerated and std::array does not. * refactor(dist): give the wavefront offset buffer a typed strided view wf_swg_max_reach strided a flat vector by hand as offs[idx(MAT_SUB)*z + s2*y + d], which was the last real concentration of untyped subscripts: 46 of the 52 in the tree, and nothing stopped idx(HAP1)*z from compiling there. ReachOffsets wraps the same caller-allocated buffer, so the reuse across calls is unchanged, and takes (mat_t, score, diagonal). The cast moves inside the accessor, where the key has already been type-checked. It immediately caught a loop still iterating the matrices as int. Genuine untyped subscripts across src/ drop from 52 to 5: four htslib gt[] reads on a C array, and the one vector whose emptiness signals no-split-found. * refactor(defs): carry #213's ac_errtype work into the typed containers The rebase onto dev brought in set_allele_errtype(vi, query) and its allele_count/ac_errtype_from_counts helpers, plus a new per-haplotype loop in calc_prec_recall and a make_tvars test helper. Three signatures still spoke in raw integers: - allele_count takes gt_t, since orig_gts and calc_gts hold the enum - ac_errtype_from_counts returns ac_errtype_t - dist.cpp's new sync-group loop and make_tvars iterate hap_t and take gt_t dev's rewritten function body and its documentation are kept verbatim; only the types change. chr20 output is byte-identical to the dev tip, so #213's GE change survives intact. * fix(defs): rename timer_t to stage_t, which POSIX already defines timer_t is a POSIX typedef in <sys/types.h>, so 'enum class timer_t' is illegal on glibc and the Linux CI build failed with 'using typedef-name timer_t after enum'. macOS does not pull that declaration into these translation units, which is why it built locally. stage_t is the better name anyway: the accessor was already Globals::stage(). Checked the other fourteen enum names against gcc 13 on glibc rather than assuming -- declaring each alongside the system headers compiles clean, with timer_t as the positive control confirming the probe detects a collision. Also worth noting the failure mode is not always an error: where a stale reference survived, glibc's timer_t resolved silently to void*. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Note
Authorship: the content below was drafted by Claude Opus 5 (an AI coding agent) and
filed via
ghunder @TimD1-bot, a bot account operated by @TimD1. It reflects theagent's analysis, not a statement authored by @TimD1.
Part of #49.
ac_errtypewas set for the query callset only, so every truth variant carriedAC_UNKNOWNandthe
GEfield on aTRUTHsample always rendered.. This records it on both callsets.Why this is needed
Not for symmetry — precision and recall are not symmetric metrics. The reason is that recall is
truth-denominated. Precision's denominator is query variants, recall's is truth variants, so the
per-truth-record decision is where recall gets computed. Genotype agreement, however, is a two-sided
property of a site: you cannot tell whether a genotype is right by looking at one callset's record
alone. The moment recall becomes genotype-aware — which is what the
gmstringency in #195 requires— the truth record needs a fact about the query.
vcfdistalready needs this and currently fakes it. Inwrite_precision_recall(src/print.cpp),a block commented "custom logic for incorrect original allele count" reads the query's
ac_errtypeto decrement a truth TP and increment a truth FN.That patch is also incomplete. It fires from inside the query loop's "variant not present on this
haplotype" branch, so it reaches only
AC_ERR_2_TO_1— a truth1|1called0|1. In the mirrorcase, a truth
0|1called1|1, the query'sorig_gtis on both haplotypes, that branch neverruns, and no correction happens. Recall today therefore penalizes false-heterozygous sites but not
false-homozygous ones. On the chr20 fixture that is 508
TRUTHsample lines corrected against 195left uncorrected.
One further constraint drove where the value is stored. The truth-to-query pairing exists only
inside
emit_sync_group; everywhere downstream the two callsets are walked in separate loops withno link between them. A site-level fact has to be recorded onto the records while that pairing is
still in hand.
Convention
A truth record carries the same value as its paired query record:
AC_ERR_<truth AC>_TO_<query AC>,with the direction absolute rather than relative to the record being written. What differs between
the callsets is which genotype supplies which count — a record's
orig_gtis its own callset'scall, and its
calc_gtis the other callset's genotype recovered by alignment.This makes the per-value annotations already in
src/defs.hreachable for the first time:AC_ERR_1_TO_0andAC_ERR_2_TO_0, annotatedTRUTH_FN, occur only on truth records;AC_ERR_0_TO_1andAC_ERR_0_TO_2, annotatedQUERY_FP, only on query records. TheGEheaderdescription, already phrased in
truth -> queryterms, becomes true of both samples unchanged.AC_UNKNOWNis unreachable on the truth side:parse_variantsonly ever handsadd_varasimple_gtthat is heterozygous or homozygous alternate, and the cluster merge preserves that.Changes
src/variant.cpp— extractedallele_countandac_errtype_from_counts;set_allele_errtype(vi, query)reads one shared table from the side matching the callset. Theprevious 31-line
else-chain collapses to three lines.src/dist.cpp—calc_prec_recall'semit_sync_grouprecords the matched query calls oneach true-positive truth variant's
calc_gt, mirroring what the line above it already does forquery variants. Guarded so the second truth-haplotype pass does not re-set a haplotype the first
recovered, which
set_var_calcgt_on_haptreats as an error.src/phase.cpp—fix_allele_countscomputes the truth value in its existing truth loop,erroring on
AC_UNKNOWNas the query loop does.src/defs.h,src/variant.h, andsrc/phase.h.Behavior change beyond the field itself
The previous
else-chain treated any unrecognizedorig_gtas zero alleles, so a heterozygouscalc_gtagainst a missingorig_gtsilently returnedAC_ERR_1_TO_0— a fabricated falsenegative. It now returns
AC_UNKNOWN, whichfix_allele_countsescalates to an error. Productioncannot reach it, since
orig_gtsis always heterozygous or homozygous alternate, andSetAlleleErrtype.UnparseableGenotypeIsUnknownNotZeroAllelespins the new behavior.No counts move
The hand-rolled truth false-negative tally is untouched and still owns the genotype error summary;
FixAlleleCounts.TruthValuesStayOutOfTheGenotypeErrorSummaryguards that. Populating truthcalc_gtsis invisible to existing readers:var_on_hapdefaults tocalc=false, so graphconstruction and truth counting still read
orig_gts; the truthGTcolumn is a literal stringkeyed off the haplotype index; and
BD/BCcome offcredit[hi][vi].Verification
716 unit tests and 98 pytest tests pass.
A/B against
devon the committed chr20 fixture, both built at-O3withmake cleanbetween:summary.vcf— identical line counts; the only header difference is##CL=(the invoked binary'spath). 703 data lines differ, every one of them changing only the
GEsubfield of theTRUTHsample: 195
.to+and 508.to-. Zero lines where any other subfield, theQUERYsample,or a fixed field moved.
precision-recall.tsv,genotype-errors.tsv,phasing-summary.tsv, andswitchflips.tsv.Checking that the convention holds rather than just that the field is populated: all 446
superclusters carrying a truth
GEalso carry a true-positive query record with the same value.Known limitation
Two oppositely-phased heterozygous query calls in one sync group union to a query allele count of 2.
The sync group is the atomic matching unit, so its query allele count is taken as the union over the
calls it holds; a comment at the call site records this.
Closes #192