Conversation
- no real changes to analysis, just minor visual updates
- simplified, no more plotting SNP, INDEL, or LARGE - stats written to json file, added script for calc
- changed to at most 2 threads per supercluster - started high-level reorganization
- this occurs after clustering but before superclustering - this will allow graph generation, with genotype information - introduced add_callset_vars() - added variantData->nc, moved ctg_variants to callset_vars
- updated Graph definition/structure with additional info
- always select INS first if multiple variants occur at same position
- main `calc_prec_recall_aln()` function is complete - added +1 to last query graph node length - added full alignment graph printing for debugging - truth string generation checks if variant occurs on hap
- truth sequence is now a path through a graph - this should enable much simple logic for parsing sync groups etc
- incorrect, but everything completes - still need to uncomment and update initial phaseset printing
- FP counting is still wrong: if calc_gt is 0|0, it's a FP
write_summary_vcf skipped every contig the query does not call on, so a truth-only contig's false negatives were absent from summary.vcf while being classified and counted correctly in truth.tsv and precision-recall-summary.tsv. The contig header line was still emitted, so the file declared a contig it carried no records for. The skip was load-bearing: the flip/swap setup that followed indexed the query lanes unconditionally, and the phase-block advance read phase_blocks[phase_block+1], which is past the end when the query has no variants on the contig and the vector holds only its single past-the-end entry. Drop the skip and derive the query phasing state defensively instead: - Initialize phase_block, block_state, and flip_error to their defaults rather than reading variant 0. The in-loop update already recomputes all three whenever the query has a variant, so the pre-loop computation was redundant duplication; a contig with no query variants now keeps the defaults, leaving truth haplotypes unswapped. - Bound the phase-block advance by the size of phase_blocks. - Narrow phase to the block that uses it. Verified: without the bound, dropping the skip trades the omission for an ASan heap-buffer-overflow at the advance. With it, chr20 output is byte-identical apart from the ##CL= binary path (100,207 records) and ASan is clean. Both one-sided-contig integration tests now pin summary.vcf, and the data README no longer documents the gap. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
add_var clamped var_qual to g.max_qual but stored gt_qual verbatim, so the two quality fields obeyed different rules and variant.h's documented range was false for any input whose GQ exceeded the cap. gt_quals has no readers: the only reads are in cluster.cpp, which passes the value straight back into add_var when variants are re-added during clustering. Every published output uses var_quals instead -- write_vcf declares only GT and PS in FORMAT and emits var_quals as QUAL, the TSV writers use var_quals, and dist.cpp thresholds and sweeps on var_quals. So the clamp cannot change results. Confirmed on the chr20 fixture, whose truth VCF carries GQ up to 866: all 11 substantive output files are byte-identical before and after. Only parameters.tsv and summary.vcf's ##CL= line differ, and solely by the -p output prefix echoed back. Also widen both field comments in variant.h from the hardcoded 0-60 to the --max-qual cap; 60 is only the default, so 0-60 was already wrong for var_quals whenever -mq was passed. No floor at 0 -- AddVar.QualNegative deliberately pins that negatives pass through, and flooring is a separate decision. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…) (#183) * refactor(print): extract compute_pr_f1 from write_precision_recall (#94) The precision/recall/F1 arithmetic sat inline in write_precision_recall, duplicated across its two loops — the per-quality curve and the NONE/BEST summary. Testing it meant running the whole file writer and diffing TSVs. Move it into compute_pr_f1(query_tp, query_fp, truth_tp, truth_fn), returning a prec_recall_f1 struct. All four counts are parameters because precision keys off the query total and recall off the truth total; a three-argument (tp, fp, fn) signature cannot express that. Both call sites now bind the struct's fields to the local names the formatting code below already uses, so the printf blocks are at zero diff. The two copies were not quite identical: the curve loop guarded F1 with `(precision+recall)` and the summary loop with `precision+recall > 0`. The single implementation keeps `> 0`. The two agree everywhere except a negative sum, which needs a negative recall, which needs the AC_ERR_2_TO_1 correction to have driven truth TP below zero (see the AcErr2To1DecrementsTruthTp tally test). In that case the curve loop used to emit a negative F1 score; it now emits 0, matching what the summary row already reported for the same counts. max_f1_score starts at 0 and only moves on a strict increase, so no negative score was ever selectable as BEST. Seven cases in test_print.cpp cover the arithmetic: - an empty query set scores precision 1, and an empty truth set recall 1, rather than dividing by zero - both callsets empty scores 1/1/1 - all-wrong query and all-missed truth gives 0/0/0, the zero-denominator branch - prec 1.0 with recall 0.5 gives F1 0.667 - asymmetric counts (query 1 TP / 9 FP against truth 1 TP / 1 FN) pin precision to the query total and recall to the truth total, which the symmetric cases cannot - a negative truth TP yields a negative recall and an F1 clamped to 0 Output-preserving: a chr20 fixture run before and after produces byte-identical precision-recall.tsv and precision-recall-summary.tsv. 568 unit tests and 65 pytest tests pass; doxygen is warning-free. * style(print): put prec_recall_f1 field comments at line ends --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
region_strs space-padded three of its four entries to a common seven-character width, and print.cpp writes those strings verbatim into the LOCATION column of query.tsv and truth.tsv. Every row therefore read `INSIDE ` with a trailing space, so any consumer comparing exactly against the documented value failed. The padding was vestigial. region_strs has three use sites and none of them is aligned console output: two are the LOCATION field itself (print.cpp:613 and print.cpp:654) and the third reads only region_strs.size(), to size the nregions counter in parse_variants. The per-region console summaries later in that function print hardcoded strings rather than indexing the table, so the width served nothing but the machine-read TSVs. Also switch OFF CTG to OFF_CTG. Every docs snapshot back to docs/v2.3.3/09-Outputs.md documents the domain as INSIDE/OUTSIDE/BORDER/OFF_CTG, but the code emitted a space, which would split the value in any whitespace-delimited reader. If aligned console output is wanted later, a %-7s at the call site gets it without baking whitespace into the data. Only INSIDE is reachable today -- variants outside the BED regions are discarded during parsing -- so the emitted change is confined to that one value. The other three become reachable once unevaluated variants are retained (#48). StringTables.RegionPadded becomes StringTables.RegionUnpadded: the equal-width loop is replaced with one asserting no entry contains a space, which is the property the TSV consumers actually need. On the chr20 fixture, query.tsv and truth.tsv are byte-identical to a dev run after stripping the trailing space, and the other eleven output files are byte-identical outright, apart from runtime.tsv timings and the invocation paths echoed into parameters.tsv and summary.vcf's ##CL=. 569 unit tests and 71 pytest tests pass; doxygen is warning-free. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
Both files were defined as the inputs optionally standardized by --realign. That flag and the VCF-normalization path are gone in v3, so the files are an unnormalized passthrough of the parsed input, written before any evaluation stage runs and carrying only PS:GT fields that summary.vcf already reports per callset sample. Removing them leaves variantData::write_vcf() and variantData::print_variant() with no caller, so both go too. The unit harness reached write_vcf() to build ParseResult::out_vcf; its position assertions now read the parsed variants directly, and wrote_pos() is renamed kept_pos() to match. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
) * refactor(dist): narrow the insertion-leap rule to variant and bypass nodes The truth-edge rule suppressing any edge that leaps a zero-width insertion locus exempted an endpoint on the reasoning that a zero-width endpoint IS the insertion's own alt or bypass node. Test that directly instead, so a zero-width node that is neither does not exempt the edge. No behavior change today: every zero-width truth node the constructor can build is already a variant alt or bypass node, since a pre-variant reference node is only emitted when var_pos > ref_pos and the trailing remainder node always spans at least one base. Prerequisite for the synthetic contig-start entry node, which is the first zero-width node that is neither. Refs #177 * fix(dist): seed the alignment at a zero-width contig-start entry node (#177) calc_prec_recall_aln() seeds Dijkstra at a hardcoded idx4 start(0, 0, 0, 0), which is correct only because ref_beg = min_pos - 1 makes node 0 a one-base left flank away from the contig start. At position 0 no flank can exist, so node 0 was the variant alt node itself: the aligner entered mid-fork and could never take the parallel allele. A matching query variant at position 0 scored FP, and its truth counterpart, whose reference-allele bypass node was equally unreachable, was left unlabeled and reported as "Unknown error type". Emit a synthetic zero-width reference node at the window start instead, whenever a side would otherwise open on a variant node. Its sentinel cell is all the flank node ever contributed to the origin, so the aligner, the backtrack, the origin sentinel and the bypass toll are unchanged. Away from position 0 the flank node is emitted first and no synthetic node is added. Refs #177 --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…201) The alignment graph's trailing node took its sequence from a substr that clamps at the contig end, but its end coordinate from an unclamped ref_end + 1. At a contig end the two disagreed: on a 10 bp contig with a SNP on the final base the trailing node spanned [2, 4) while holding zero bases, and graph->ref was two characters that the coordinate 4 overran. Clamp the window once, in win_end, and use it for the trailing query node, the trailing truth node, this->truth and this->ref. Every node's coordinate span now equals the length of the sequence it holds, with zero-width insertion variant and bypass nodes the deliberate exception. At a contig end the trailing node becomes a zero-width sink, the mirror of the zero-width entry node #177 added at position 0. No behavior change: nothing reads the trailing node's end coordinate. Both alignment endpoints are built from qseqs/tseqs lengths rather than coordinates, edit distances go through the already-clamped graph->ref and graph->truth, the trailing node is the sink so its end coordinate is never matched against any node's begin coordinate, and the insertion-leap rule additionally requires tidxs >= 0 || tskips >= 0, which the trailing node never satisfies. Its sequence was already empty at a contig end, so the endpoint cell was already the sentinel. A chr20 before/after run is byte-identical across every scored output file, which covers the loop restructure; the fixture's last variant sits 158 kb from the contig end, so the clamp itself is covered by the new tests rather than by it. Retire the TODO asking why the +1 is needed "compared to generate_ptrs_strs()". That function was removed in 232036e and reached exactly as far, via `for (int ref_pos = beg_pos; ref_pos <= end_pos; )` and `ref_end = end_pos+1`, so there was never an asymmetry between them to explain. The remaining +1 is the window's right flank and is left alone: unlike the coordinate, the trailing node's sequence is read, so shortening it would be a scoring change. Add the coverage the issue asks for, mirroring #177's: five GraphCtor tests over a variant on a contig's final base and on the second-to-last (all five fail without the clamp), three PrecRecall tests pinning TP, FN and FP there, and a contig_end_snp integration scenario placing a SNP on the final base of sc1. Refs #189 Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…int_* formatters (#182) (#202) Adds 96 cases to tests/unit/src/test_globals.cpp, the last major gap in the unit suite. No source file is modified. parse_args (83 cases): the argc < 4 short-circuit group, distinguished by output rather than exit code alone; the three mandatory positionals; the verbosity pre-pass; every optional flag's happy, missing-value, bad-numeric, and bound-violation cases; the main-loop -h/--version/-ci paths that do not exit; unrecognized options; and the four post-loop cross-field checks. set_thread_ram_steps (7 cases) is driven directly rather than through an argv array, since the constructor and parse_args both call it. print_version, print_usage and print_citation (6 cases) assert on captured stdout, including that the usage text interpolates live settings and omits the four flags whose usage lines are commented out. Three harness pieces, all local to the test file: - ArgsFixture writes an openable query and truth VCF and points the reference at the committed tiny.fasta, so a case states only the flags it varies. - parse() copies the arguments into mutable buffers, since parse_args takes char**, and closes the reference FASTA that parse_args opens and never closes. - parse_showing_stdout() redirects stdout onto stderr inside the death-test child, because a death-test matcher only sees stderr while print_usage, print_version and print_citation write to stdout. Without it the argc < 4 group could assert only an exit code, making -v, -h and -ci indistinguishable. Two cases document current behavior rather than enforcing it, each marked in-file with the issue that will change it: a trailing comma in --filter appends an empty filter name (#199), and -b's catch block reports an invalid filename for what is actually invalid file contents (#200). A third pins that the missing-value error for --verbosity is unreachable, since the pre-pass loop stops at argc-2. Closes #182 Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
add_var took 27 positional parameters, 17 of them defaulted, so inserting a parameter anywhere but the end silently rebound every argument after it. #116 had already shipped that bug once: a caller's supercluster argument became rec_idx and compiled clean under -Wall -Wextra. Every field now travels in a var_fields aggregate, matched by designator instead of position, with the 12 per-haplotype fields nested as hap[HAPS] so a hap1/hap2 transposition is no longer expressible. Members with no default are required: omitting one is a build failure under the newly added -Werror=missing-field-initializers (GCC; clang does not diagnose it). Call sites name the type -- add_var(var_fields{...}) rather than add_var({...}) -- because GCC 13.3, the compiler on ubuntu-24.04 and so on CI, rejects a bare designated-initializer list as a function argument once any member is initialized from a non-constant expression. Naming the type sidesteps that and still enforces the required members. -std=c++20 does not help, so the C++17 standard level is unchanged. The four cluster.cpp merge sites hand-copied 13 parallel vectors each, which is the shape that produced the #116 miswiring; they now read one var_fields through the new ctgVariants::get_var accessor. Pure refactor: on the chr20 fixture every output is unchanged, the two files recording the invocation aside. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…nate (#200) (#203) * fix(bed): report the file, line, and field for a malformed BED coordinate (#200) std::stoi on the start and stop columns threw std::invalid_argument out of the bedData constructor, where -b's catch in parse_args reported it as "Invalid BED filename provided" and discarded the exception. The filename was fine; the contents were not, and neither the offending line nor the column survived. Every other BED failure -- an unopenable path, and the flipped, zero-length, unsorted and overlapping cases in check() -- calls ERROR(), which exits rather than throws, so the catch was dead for all of them. Parse each coordinate through a helper that names the field, its 1-based line and the BED file, building the diagnostic where that context exists. The catch in parse_args now reports e.what() like the -p and -f branches, leaving it a formatting net rather than a mislabel. The helper also requires the whole field be numeric, so a partially-numeric coordinate such as "8bp" is rejected instead of silently truncated to 8. Refs #200 * test(bed): pin the malformed BED coordinate error routes (#200) Five constructor cases cover what reaches the new coordinate parser: a non-numeric field, a partially-numeric one, a value too large for int, an absent third column, and a trailing blank line. Each asserts the field and line number in the message, so a regression to a filename-shaped diagnostic fails here. Refs #200 * test(globals): assert the new BED coordinate message in the -b handler tests (#200) BedMalformedLineErrors pinned "Invalid BED filename provided" as current-behavior documentation, with a note to update it here. A malformed coordinate is now reported by the constructor, so the assertion becomes an enforcing one against the field-and-line message, and both tests' comments drop the mislabeled-handler framing: the catch block is unreachable for every BED input, not just for an unopenable path. Refs #200 --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
… list (#199) (#205) * fix(globals): drop empty filter names from a comma-separated --filter list (#199) parse_args tested filters_ss.good() before the getline that consumed the field, so the read past the last name left filter empty and the unconditional push_back appended it anyway: -f PASS, yielded {"PASS", ""}. The empty name never resolved to a FILTER id, so parse_variants warned once per callset about a filter nobody asked for. Drive the loop off getline's result and skip empty fields, which covers the leading (,PASS) and interior (PASS,,LowQual) cases too. An argument that names no filter at all now errors, since an empty g.filters means keep every variant, the opposite of the request. * test(integration): pin that a trailing --filter comma adds no filter (#199) Runs the swallowed_snps fixture with -f "PASS,". Every record in both callsets is PASS, so the results must match the run without the flag and stderr must stay free of the 'Filter '' not found' warning the empty field used to produce. Fails on the pre-fix binary. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…loidy (#102) (#204) Ploidy was one value per contig, inferred from the first record's GT length and only enforced thereafter. On a chrX carrying both PAR (diploid) and non-PAR (haploid) calls every output genotype was forced into whichever shape the first record happened to have, and the disagreement warning was suppressed by name for contigs ending in X -- so chrX was silently wrong while chrY warned on every record after the first. The ploidies vector added by #101 makes mixed ploidy within a contig representable, so the per-contig value can go. variantData::ploidy becomes observed_ploidies, a vector<set<int>> parallel to contigs and lengths. The set-once-then-enforce block in parse_variants collapses to a single insert, and wrong_ploidy_total and its summary WARN are deleted along with the chrX exemption. The insert stays ahead of the record-dropping logic, so every record in the input contributes its ploidy even if it is later dropped as REF, .|., oversized, or overlapping: the comparison below catches a whole-callset mistake, which should be reported even when the offending contig's variants are all filtered out. intersect_contigs compares the two observed-ploidy sets rather than two scalars. This also fixes a spurious warning: a contig injected empty was given ploidy 0, which compared unequal against the other callset's real ploidy, so a BED contig present in truth but absent from query emitted "has ploidy 2 and ... has ploidy 0". BedAddsEmptyContig triggered that warning and passed anyway, capturing stderr without asserting on it. A contig that observed no ploidy at all is now skipped rather than treated as a disagreement, and the ploidy-inheritance copy that existed only to keep the scalar comparison quiet is gone. superclusterData::ploidy and phaseblockData::ploidy carried that scalar to the summary VCF writer and nothing else, so both are removed. ploidy= leaves the ##contig line -- it is not a VCF-spec contig attribute, and each record's GT now carries its own. The five GT ternaries in write_summary_vcf go through one hap_gt() helper reading each callset's own per-variant ploidy, so a haploid truth against a diploid query renders 1 against 0|1 instead of forcing both into one shape. Ploidy 0, the unknown sentinel, still falls to the diploid branch, matching what the old ploidy == 1 test did for an unset contig. simple_gt is deliberately unchanged. Extending it to take GT_ALT1 for haploid records would look tidier, but var_on_hap returns true for both haplotypes on GT_ALT1 (variant.cpp:231, :233), so a haploid variant would begin counting twice, and set_allele_errtype has no GT_ALT1 branch (variant.cpp:152-177), so haploid variants would fall through to AC_UNKNOWN. A separate vector needs no audit of the clustering, alignment, or phasing consumers of the genotype enum. Half-call rendering is accepted as-is: simple_gt becomes GT_ALT1_REF for a half call, so a 1|. record round-trips through summary.vcf as 1|0. Changing it would need either a fourth simple_gt value or the half call carried alongside ploidies, both of which complicate the ploidy design for a rendering detail. Recorded so the 1|0 is not later read as an oversight. No unit test reached write_summary_vcf before this, so test_phase.cpp gains a write-and-read helper composed from the existing make_phaseblockData and TempDir/read_text pieces. Six tests were written first and each watched fail for the right reason -- haploid records rendering 1|0, the spurious ploidy-0 warning, and ploidy= in the header -- before any production edit. The chrX case is asserted with each ordering first, since the bug being removed was "whichever record came first wins". On the chr20 fixture, all eleven output files are byte-identical to a dev run except summary.vcf, whose only change is ##contig losing ploidy=2; the ##CL=, parameters.tsv, and runtime.tsv diffs are invocation paths and timings. 578 unit tests and 71 pytest tests pass from a clean build, with no new warnings under -Wall -Wextra. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
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. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…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>
add_variants had no callers anywhere in src/ or tests/unit/src/. It built variantData entries from a WFA backtracking-pointer CIGAR, a path nothing in the pipeline takes. print_wfa_ptrs and the ptr_t enum are kept: the WFA matrix printer is worth retaining for debugging, and it is now the only consumer of ptr_t. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
orig_gt and calc_gt read as a before-and-after pair on one sample: the original genotype and the recalculated one. That reading is natural and wrong. They are genotypes of different callsets. On a query record orig_gt is the query's own call, indexed by query haplotype, while calc_gt is the truth genotype recovered by alignment, indexed by *truth* haplotype. On a truth record the roles reverse. That split in index space is load-bearing rather than incidental. matched_gt_is_swapped() exists only to reconcile the two, and matched_hi = is_swapped ? other_hap(hi) : hi is how every per-haplotype array gets subscripted downstream. "calc" named the process, which distinguishes nothing here since every field in ctgVariants is computed from something, and said nothing about whose haplotypes the bits index. orig/matched does not read as a temporal sequence, so the before-and-after misreading goes away, and "matched" names the content: the genotype matched on the other side. It stays accurate at the edges, where nothing matched and 0|0 reads correctly. The accessors carrying the old term move with it: calcgt_is_swapped -> matched_gt_is_swapped, set_var_calcgt_on_hap -> set_var_matched_gt_on_hap, and var_on_hap's `calc` parameter -> `matched`. orig_gt is left alone; it is accurate, used far more widely, and was never the confusing half of the pair. Nominal only: 217 insertions against 217 deletions, with no line added or removed. Every chr20 output file is byte-identical to the dev tip, and the unit suite passes unchanged at 726. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
… modes (#206) (#222) parse_args accepted -h, -v and -ci in two places with two meanings apiece, and the ambiguity was real only because both meanings could occur in one invocation. Deciding the mode by argc before any flag is read makes them disjoint, so the overload stays and the ambiguity goes. Informational mode, entered when fewer than the three mandatory arguments are present, accepts only -h/--help, -v/--version and -ci/--citation. None takes a value; each prints and exits 0. Anything else, including --verbosity, warns 'Invalid usage.' and prints usage. Evaluation mode accepts every other flag, each consuming exactly one following token, and rejects the three informational flags by name: Option '-h' is informational only; use it without the mandatory arguments rather than the generic "Unexpected option '%s'", which reads as "no such flag" and misdirects a user who has made a scoping mistake rather than a typo. Removing -n/--no-output-files is what makes "every flag takes a value" exceptionless: it was the only main-loop flag that took none. Its one writer gone, g.write would have been permanently true, so the eight always-taken guards, the write_outputs row in parameters.tsv, and b2s() -- which existed for that one call site -- all go with it. The re-indentation from unwrapping those guards dominates the diff. One behavior change beyond the contract, and it is a bug fix. The verbosity pre-pass looped to argc-2, so i could never equal argc and the missing-value ERROR was unreachable: a trailing -v was silently accepted and the parse succeeded with verbosity unchanged. Starting the loop at 4 is correct now that -v is legal only past the positionals, and it makes the guard fire. parameters.tsv loses its write_outputs row. The file is row-keyed rather than columnar, so nothing else shifts; only a consumer reading that key is affected, and the value would document nothing once -n is gone. The archived docs/v2.* trees still describe -n as those releases shipped it and are left untouched. print_usage is the only live documentation of the flag set: it gains a second usage line for the informational form and moves -h/-v/-ci out of Miscellaneous into their own section. 731 unit tests and 98 pytest cases pass. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…call_wrapper (#221) (#223) thread2 once gated spawning one thread per haplotype for the two calc_prec_recall_aln calls (a7a111d). That block was dropped in cc4a810 when the alignments collapsed into the per-haplotype evaluate_variants loop, leaving a parameter nothing reads and the only -Wunused-parameter warning in dist.cpp. Reinstating the optimization as-is would race: the threaded unit today is evaluate_variants (alignment plus labeling), and set_var_matched_gt_on_hap is a read-modify-write on a shared per-variant slot that ERROR()s on an unexpected transition. It would also gain nothing on realistic input, where every supercluster lands in bucket 0. The old code stays recoverable via git show a7a111d:src/dist.cpp. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
* perf(ci): build vcfdist once per Tests run instead of three times (#225) Of a ~96s Tests job only 4.6s ran tests; ~60s was pytest_sessionstart recompiling both C++ projects serially, and because the unit-test build's clean removed $(SRC)/*.o, the 8 vcfdist translation units were compiled three times per run. Both builds now use identical flags (-g -O3), so they share the src/ objects rather than each owning a differently-flagged copy: src/Makefile drops -pg -O1 (gprof profiling needs -pg added back, per the comment there), the test build no longer cleans $(SRC)/*.o, and conftest returns early on VCFDIST_SKIP_BUILD so the workflow can build both projects once with -j$(nproc). Local test runs keep the rebuild, now parallel. Google Test is linked prebuilt on both platforms: Ubuntu's libgtest-dev ships libgtest{,_main}.a, so CI stops compiling gtest-all.cc, which is what pays for -O3 costing more to compile than the test build's previous -g. The two apt-get install runs are merged, and pip is replaced by uv. Verified in an ubuntu:24.04 container with GCC 13.3: no new warnings, the test build recompiles zero src translation units, 731 unit tests pass, and chr20 evaluates in 1.06s. * fix(ci): pin setup-uv to v9.0.0 setup-uv stopped publishing floating major tags after v7, so @v9 does not resolve and the job fails at Set up job. * fix(ci): let uv install into the runner's system Python The runner's /usr Python is EXTERNALLY-MANAGED; pip only worked because the image sets break-system-packages in /etc/pip.conf, which uv does not read. * fix(ci): run pytest from a uv venv uv refused the runner's system Python as EXTERNALLY-MANAGED, and with --break-system-packages it then hit permission denied on /usr/local/lib/python3.12/dist-packages. A venv sidesteps both. * perf(ci): compile the test-only objects without optimization At -g -O3 the 11 test translation units took 42s on a 4-core runner, making them the job's bottleneck; test_variant.cpp alone took 19s. They are not shared with src/ and the suite runs in ~0.1s, so optimizing them buys nothing. The 8 shared src objects keep -g -O3, which is what object sharing requires. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…ums (#224) * refactor(defs): drop ALT1 from the stored gt_t value names (#219) The evaluation vocabulary never tracked which alternate allele was involved, so ALT1 was misleading in the values a variant's genotype can hold. GT_REF_ALT1 -> GT_REF_ALT GT_ALT1_REF -> GT_ALT_REF GT_ALT1_ALT1 -> GT_ALT_ALT Mechanical: no behavior change. GT_REF_REF already read correctly and keeps its name. The parse-only enumerators keep theirs until the parse enum is split out. * refactor(variant): split the parse-time genotype vocabulary out of gt_t (#219) Add gtparse_t, a 9-value enum used only for the Genotypes: summary, and extract the classification out of parse_vcf into a pure classify_gt() that can be unit tested with a hand-built GT array. The parse enum is allele-index-agnostic: A and B stand for any alternate, so 1|2, 2|1, and 1|3 all report as A/B. With missing alleles handled first the diploid space is exhaustive over four cases on the two allele indices, and the haploid space over missing, zero, and nonzero, which retires GT_OTHER as a classification outcome. Three consequences: - Haploid '.' now has its own bin rather than sharing the diploid './.' one. - multi_total drops 0|2, which is heterozygous and is not split across haplotypes, so it never belonged in a total described as "homozygous and multi-allelic ... split for evaluation". - The heterozygous phasing-imbalance WARN is gone. Its two inputs merge into a single 0/A bin, leaving nothing to compare. No evaluation logic is touched: the per-haplotype loop still derives its genotype directly from bcf_gt_allele(). All ten chr20 fixture output files are byte-identical. * refactor(defs): narrow gt_t to the four genotypes a variant can hold (#219) orig_gts and matched_gts only ever hold GT_REF_REF, GT_REF_ALT, GT_ALT_REF, or GT_ALT_ALT. Dropping the other seven values makes the type match its domain, shrinks gt_strs from 11 entries to 4, and removes code that could not run: - var_on_hap no longer tests GT_ALT1 on both haplotypes. That branch was unreachable, and wrong if reached: a haploid contig populates HAP1 only, so claiming presence on HAP2 would be incorrect. - allele_count's switch is now exhaustive, so it no longer returns -1 for a genotype carrying no diploid allele count. AC_UNKNOWN remains reachable through the both-counts-zero case. The terminal ERROR() in matched_gt_is_swapped and set_var_matched_gt_on_hap are now guards rather than reachable paths; both are kept, since each is the required terminal branch of an if-chain. The nine tests that reached those paths by constructing an out-of-domain genotype go away with the values they used. * test(variant): cover classify_gt directly, including the bins chr20 never hits (#219) The chr20 fixture only carries 0|1, 1|0, 1|1, and 1|2, so the parse-time classification of haploid '.', 0|2, 2|2, 1|3, .|1, and 1|. had no coverage. Now that the classification is a free function it can be called with a hand-built GT array instead of a VCF. Also asserts the two properties the bin merges depend on: phasing cannot change a classification, and the diploid space is exhaustive over every allele pair from missing through 3. * style(variant,phase): drop trailing whitespace on the renamed genotype lines (#219) --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…d file as empty (#227) (#228) * fix(variant): fail on a critical bcf_read instead of reading it as end of file (#227) htslib returns -1 from bcf_read at end of file and < -1 on a critical error, but the loop tested for 0 and treated both alike. A record htslib refused to parse ended the read as though the file had simply finished, so a malformed query VCF was scored as an empty one: exit 0, recall 0, precision 1, and a full precision/recall table indistinguishable from a caller that found nothing. Split the three cases. < -1 is fatal and names the callset, the file, the record number, the position, and the decoded bcf_strerror() text. A non-critical errcode on an otherwise-successful read is warned and counted rather than made fatal, since BCF_ERR_TAG_UNDEF from an undeclared INFO or FORMAT tag is common in VCFs vcfdist evaluates correctly today. One such errcode cannot be recovered from. BCF_ERR_CTG_UNDEF means htslib appended the contig to the header as it parsed, so rec->rid runs one past the ctgnames array captured before the loop, and ctgnames[rec->rid] reads out of bounds. The added test segfaults without the bounds check. * fix(variant): distinguish the four GT read failures (#227) bcf_get_format_int32 returns -2, -3, -4 and 0 for four different reasons, and all four printed "Failed to read QUERY GT at chr1:99", naming neither the cause nor the fix. -3 was the worst of them, because it is not an error. It means the header declares GT but this record's FORMAT column omits it, which is legal VCF: the FORMAT column is per-record and GT is not mandatory. Aborting rejected a spec-conformant file. GT was also the only tag treating -3 as fatal; GQ reads it as quality 0 and PS counts it and continues. The record is now dropped and counted, the same treatment a stated .|. no-call already gets. It is not given an assumed genotype: unlike GQ's 0 and PS's one-phase-set-per-contig, there is no honest default, and inventing an allele the record never stated would score it as a TP or FN. The skip precedes both classify_gt() and the observed_ploidies insert, so a dropped record reaches neither. -2 stays fatal, since the spec fixes GT's type as String and no retry is possible the way GQ retries as float, but it now says the header is at fault. -4 and 0 keep one generic fallback message and get no test, because no input reaches them: a record malformed enough for -4 fails inside bcf_read first, and 0 is not a documented success return. * test(integration): assert a malformed query VCF exits nonzero (#227) The unit death tests prove ERROR() fires. Only the integration harness proves the exit code a calling pipeline observes, which is the symptom the issue is about: the same run previously exited 0 with a full precision/recall table. The fixture needs an explicit force-add, since *.vcf is gitignored and the un-ignore rules cover tests/unit/data/ alone. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
* refactor(phase): emit one summary VCF record per variant (#103) write_summary_vcf() wrote one record per haplotype per variant, so a homozygous call was re-split into two records after the cross-haplotype merge had already collapsed it into a single entry. Each variant is now written once. GT reports the caller's own claim (orig_gt), not vcfdist's inferred calc_gt, so the record does not present an inference as a call. The per-haplotype fields (BD, BC, RD, QD, BK, SG) become comma-separated lists carrying one value per GT allele, in GT allele order; the evaluation lanes are indexed by calc_gt's haplotypes, which calcgt_is_swapped() reports may be the reverse of orig_gt's. A reference allele was never evaluated, so it reports "." in every such field, which keeps the number of TP/FP/FN values emitted exactly as it was. Those fields are declared Number=. rather than Number=P. Number=P is semantically correct but is a VCF 4.4 addition that htslib only supports from 1.23, so consumers on older bcftools or pysam would report a cardinality error. Number=. produces byte-identical records and gives up only the declared cardinality, so the count and order are stated in the field descriptions instead. A het-alt (1|2) record stays two co-located records: parsing splits it into two entries whose alleles normalize independently and may not even share a position, and nothing rejoins them. precision-recall-summary.tsv and every other TSV output are byte-identical across the change on the chr20 fixture; the counting convention is #49's change, not this one. * refactor(variant): address review of the one-record-per-variant writer (#103) - Drop credit_str() and its fixed-size char buffer; std::to_string(float) is specified to render exactly what sprintf("%f") does, so the records are byte-identical and no custom helper is needed. - Add ploidy_t {PLOIDY_HAPLOID = 1, PLOIDY_DIPLOID = 2} to defs.h and type var_fields::ploidy and ctgVariants::ploidies with it. There is no unknown ploidy to represent: a polyploid record errors out at parse time, and a record whose VCF declares no GT tag reports ngt == -1 and is assumed monoploid, so every variant reaching add_var() is called on 1 or 2 haplotypes. An omitted ploidy now defaults to diploid rather than 0. - Extract display_gt(), mapping the caller's own orig_gt and ploidy to the GT string the sample reports. - Name the loop over haplotypes hi, and the lane it resolves to through the matched_gt swap hi_resolved; the two differ only when the genotypes are swapped, which is exactly what the swap test pins. - Use std::string, not const char*, for the shared FORMAT description suffix. BK still reports lm for a credit at or above the threshold. The am tier is part of the match-tier ladder #49 introduces, and adding it here would change the counts this PR holds byte-identical. summary.vcf is byte-identical across these changes on the chr20 fixture, as are all eight TSV outputs against dev. * docs(variant): describe ploidy as allele count, and rename hi_matched (#103) Ploidy counts the alleles a genotype declares, not the haplotypes carrying the variant: a het 0|1 is called on one haplotype but is diploid. Reword ploidy_t and both ploidy fields accordingly, and drop the same phrasing from print_var_sample() and the test builders. Rename hi_resolved to hi_matched, since what it resolves through is the matched_gt swap. The haplotype count now reads straight off the ploidy rather than re-deriving it from a comparison. Comments only, apart from the rename: summary.vcf and all eight TSV outputs are byte-identical on the chr20 fixture. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…#229) (#239) * refactor(phase): write the summary VCF with htslib instead of fprintf (#229) The summary VCF was emitted with fprintf, so every VCF semantic on the output side had to be re-implemented on strings. Build each record as a bcf1_t and write it with bcf_write instead, so htslib owns the encoding. write_summary_vcf() now opens the file with hts_open and writes a header built by summary_vcf_header(); print_var_info/print_var_empty/ print_var_sample become set_var_record(), which sets the fixed fields, and var_sample_fields()/empty_sample_fields(), which return one sample's FORMAT values for set_record_samples() to write. Per-allele values carry htslib's missing sentinels, and a sample of lower ploidy is padded with the end-of-vector marker rather than a shorter rendered list. The output changes in one visible way: BC is a Float, so htslib renders it compactly ("1" and "0.8", not "1.000000" and "0.800000"). Over the 72,585 records of the committed chr20 fixture that is the only field that differs; every other field is byte-identical and every BC value is unchanged at float32 precision. The header also declares PASS before fileDate, since bcf_hdr_init emits it first. Test assertions over the rendered floats are regenerated accordingly. * fix(tests): reach htslib through variant.h in test_phase.cpp The build copies htslib's headers into src/, so a bare "htslib/vcf.h" only resolves for a file in that directory; from tests/unit/src it built locally only because an unrelated Homebrew prefix happens to sit on the default include path. Include src/variant.h instead, as test_variant.cpp does. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…#240) * fix(bed): read BED files through htslib so gzip and bgzip work (#230) bedData::bedData() read with an ifstream, so a compressed BED could not be parsed. The fopen existence check succeeded on a .bed.gz, the first line then decoded to binary garbage, and parse_coord() reported an invalid coordinate -- naming the wrong cause, since the coordinate was not malformed. ERROR() exits, so the try/catch around the -b constructor never ran. Replace the ifstream with hts_open() plus hts_getline() on a kstring_t, which decodes plain, gzip, and bgzip transparently and detects the encoding from the file's leading bytes rather than its extension. A read failure below EOF is now reported rather than being mistaken for a short file, so a truncated compressed BED no longer parses as a partial region set. The malformed-coordinate ERROR paths are unchanged. htslib was already a link dependency; bcf_open() reads the VCFs. * build(tests): put src/ on the unit-test include path so htslib resolves A quoted include is resolved against the including file's own directory before any -I path, so "htslib/bgzf.h" in tests/unit/src/test_helpers.cpp searched tests/unit/src/htslib/ and stopped. It built locally only because Homebrew's htslib is on the default search path; CI copies the headers into src/ instead of installing them, so the test build failed to find bgzf.h. Add -I$(SRC) to TEST_CXXFLAGS, giving the test objects the same htslib resolution the src objects get for free by living in src/. Verified both ways: with the headers copied to src/htslib (the CI layout) they resolve from there, and with that directory absent they fall back to the system path. CXXFLAGS is untouched, since it must stay identical to src/Makefile. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
parse_variants() rejected a VCF that returned to a contig it had already left, but nothing checked the order of positions within a contig. A record that moved backwards was dropped one at a time by the overlap filter, which warns only above the default verbosity, so the run finished successfully, wrote a full set of outputs, and every denominator silently omitted it. Track the previous record's position alongside prev_end and prev_type, resetting it on a contig change, and ERROR when a record's position decreases. The check runs before FILTER and quality filtering, since being sorted is a property of the file rather than of the records that survive. The ordering that held before was incidental: the overlap filter skips a record unless pos >= prev_end[hap], which is per-haplotype, applies only to accepted records, and runs after the region-membership lookup that #47's monotonic per-region cursor will depend on. Enforcing the precondition here means nothing downstream has to assume it. Equal positions are still accepted -- a site split across rows is ordinary in a sorted VCF -- and so is a contig that starts before the previous one ended, since positions are compared only within a contig. The integration case is worth its own fixture rather than only a unit test precisely because the failure mode it replaces was a successful run. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
BREAKING CHANGE: --max-supercluster-size moves from the short form -s to -sc, matching the existing multi-character convention of -sv, -mq, -ct, -md and -ci. This frees the mnemonic -s for the stratification flag added under #47. No -s alias is retained. A silent alias would let an existing '-s 15000' invocation be reinterpreted as a file path once -s means stratification, so bare -s is rejected outright and an old command line fails loudly. Only the short form changes, so the long form and every prose reference to --max-supercluster-size still read correctly. Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…ts guard (#231) (#241) * refactor(bed): extract classify(), add merge(), and drop the bed_exists guard (#231) Three zero-output-change edits to bed.cpp, all prerequisites of the #47 membership sweep. contains() interleaved locating a variant -- start_idx by upper_bound on starts, stop_idx by lower_bound on stops -- with classifying it from how those indices relate. classify() now takes the decision tree plus the before-all/after-all BED_OUTSIDE early returns, which are load-bearing: a variant left of the first region has start_idx -1 and would otherwise read as BED_BORDER. contains() keeps both binary searches and calls classify(), so it returns exactly what it did before. #47 adds a second way of locating a variant, a monotonic cursor; without the split there would be two copies of one rule, free to drift. contains() opened with `if (!g.bed_exists) return BED_INSIDE;` -- a method on one bedData consulting a global flag about a different one. For a stratification region set that is wrong: membership would depend on whether -b was supplied. The check moves to the sole -b call site in parse_variants(), which also makes contains() testable without global setup. merge() sorts each contig's intervals by start, merges those that overlap or abut, recomputes size, and reports the number coalesced at verbosity >= 2. A new bedData(bed_fn, merge_overlaps) selects it in place of check(). No caller passes true yet: -b keeps the strict check(), which is right for an evaluation region whose malformation silently changes every denominator and wrong for a third-party region set we neither author nor control. Merging is not tidiness -- contains()' binary searches require sorted, disjoint intervals, so an unmerged set returns wrong answers. BedContains.NoBedInside asserted exactly the deleted early return, so it is retargeted to the parse_variants() call site as ParseVariants.NoBedEverythingInside, which also pins that a leftover g.bed cannot filter an unrestricted run. The g.bed_exists = true scaffolding the other BedContains cases carried is now dead setup and is deleted. The existing contains() return-value cases are unchanged, which is the guard on the extraction; reordering classify()'s early returns fails four of them. * refactor(bed): rename merge() to normalize(), and check after it (#231) Review feedback on #241. merge() becomes normalize(), since it sorts as well as merges, and the constructor's merge_overlaps parameter becomes normalize to match. check() now runs unconditionally, after normalize() rather than instead of it. Normalizing repairs disorder and overlap and nothing else, so it does not weaken the validation it precedes: a flipped or zero-length interval is left exactly as it was found and check() still rejects it. normalize() warns when the regions were not already sorted -- sorting is a repair, and a file needing it is not one check() would have accepted -- and both of its messages now name the BED file, since #47 normalizes many region sets and a message about one of them has to say which. bedData gains a filename member, set by the file constructor and empty for an in-memory bedData, to carry that name. BedCtor.MergeOverlapsSkipsCheck is renamed NormalizedOverlapPassesCheck and now asserts that the check still runs and passes, rather than that it was skipped. NormalizeStillRunsCheck pins the other half: a flipped and a zero-length region are still fatal with normalize enabled. The message tests load from a file rather than building in memory, since what they assert on is the file being named. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…234) (#244) * feat(cli): add -st/--stratification manifest loading and validation (#234) Part of #47. The flag and the loader only: region sets are parsed, merged, and resident in Globals, and nothing reads them yet, so no output changes but --help. -st takes one hap.py-compatible manifest TSV whose two tab-separated columns name a stratum and its BED file, with further columns ignored. Relative BED paths resolve against the manifest's own directory rather than the working directory, so the GIAB manifests work unmodified. Blank and '#'-prefixed lines are skipped, and load order is manifest order, which is also the output row order the strata will take. Each region set is normalized at load, since a third-party set we neither author nor control would otherwise have to be sorted and non-overlapping for contains() to answer correctly; normalizing does not exempt it from validation, so a malformation sorting and merging cannot repair is still fatal. Five conditions are fatal: an unreadable manifest, a line naming fewer than two fields, a repeated stratum name (which would make an output row ambiguous), a stratum named '*' (reserved for the all-regions row), and an unreadable stratum BED, which names both the manifest-relative and the resolved absolute path so that a path resolved somewhere unintended is distinguishable from a missing file. A manifest that parses to zero strata warns instead, and the run proceeds as though -st were absent. The realistic failure mode is an assembly or contig-naming mismatch, which would otherwise emit rows of zeroes that read as genuine results, so each region set's contigs are compared against the reference FASTA's once it is loaded. A set sharing no contig is named; every set failing at once is diagnostic rather than incidental, so that case names the likely cause instead. * fix(print): record stratification, supercluster size, and verbosity in parameters.tsv The new -st/--stratification manifest path was not written to parameters.tsv. Auditing every CLI-settable parameter against the writer surfaced two pre-existing omissions with the same shape: -sc/--max-supercluster-size (added in 41cc0a0, never wired in) and -v/--verbosity. No stale rows: all 22 existing keys still map to live Globals members. * refactor(print): write parameters.tsv one row per fprintf, and record nstrata The writer emitted all 26 rows from a single fprintf, whose format string was split across seven fragments and whose arguments ran across five lines, so a key and its value were never adjacent and a row could be added to one without the other. Each row is now its own call, with its key, conversion, and value on one line. Two changes to the rows themselves. The stratification manifest key becomes 'stratification', matching the --stratification flag that sets it. And nstrata is now recorded: the manifest path alone does not say how many region sets were loaded from it, which is precisely what a reader needs to tell a fully-loaded run from one whose manifest parsed to nothing and warned. WriteParams.WritesEveryKeyOnceInOrder pins the whole key list in order, since one call per row makes a dropped or duplicated row a one-line edit away. --------- Co-authored-by: Tim Dunn <timdunn@fulcrumgenomics.com>
…245) * feat(variant): add the match-tier criterion enums and ladder (#193) Add credit_t, allelecount_t, phasematch_t, and matchtier_t, plus the pure match_tier() ladder and the three derivations that feed it. Nothing calls them yet, so vcfdist's output is unchanged. The ladder is a conjunction, so pm implies gm implies am implies lm. MaxAlleleCredit supplies both lower rungs and MinAlleleCredit neither: Min is the stronger predicate, so putting it on the looser rung would rank a genotype error below a weak partial match. ac_errtype_to_allele_count() is total. The 0 -> N pure false positives and N -> 0 pure false negatives map to GAIN and LOSS rather than falling outside the ladder, which changes no tier (only EQUAL reaches gm) but carries the direction of the error. get_phase_match() tests zygosity before phasing, since homozygous and haploid variants keep the PHASE_NONE default and would otherwise report as unphased. PHASEMATCH_NOT_HETEROZYGOUS covers both, as neither occupies a phased pair. * fix(variant): treat homozygous reference as not heterozygous get_phase_match() now asks which genotypes are heterozygous rather than which are not, so the test is total over gt_t: GT_REF_REF joins GT_ALT_ALT and the haploid call in PHASEMATCH_NOT_HETEROZYGOUS, and a genotype added later cannot silently acquire a phase it does not have. parse_variants() stores one variant per non-reference allele, so a 0|0 genotype cannot currently reach the function. It is answered rather than left to fall through to the PHASE_NONE default, which would have reported it as unphased. * fix(variant): name the homozygous genotypes literally in get_phase_match() Test GT_REF_REF and GT_ALT_ALT directly rather than asking which genotypes are heterozygous. New gt_t values are more likely to be heterozygous than homozygous, and those must fall through to the phase comparison rather than be swept into PHASEMATCH_NOT_HETEROZYGOUS by a heterozygous-genotype whitelist. --------- 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.
v3.0.0 Deliverable Tracking
Progress against the v3.0.0 SOW, one issue per deliverable.
Design
vcfdist Extensions
Benchmarking, Release, and Infrastructure
Potential Future Work (not required for v3.0.0)