refactor(phase): emit one summary VCF record per variant (#103) by TimD1-bot · Pull Request #208 · TimD1/vcfdist

TimD1-bot · 2026-08-05T21:34:56Z

Note

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

Fixes #103. Step 3 of 8 toward #48. Was stacked on #204 (#102), whose per-variant ploidy this
change renders each record's GT from; #204 has since merged, so this is rebased onto dev and
targets it directly.

Problem

write_summary_vcf emitted one record per haplotype per variant, through four
for (int qhi = 0; qhi < HAPS; qhi++) loops. By the time the writer runs, a homozygous
variant is already a single entry — load_and_merge_callset_vars_across_haps collapsed
its two parse-time copies by exact (pos, ref, alt) equality and recorded GT_ALT1_ALT1
(src/cluster.cpp:137-150). The writer then re-split it. So this is not a merge, it is
declining to re-split.

Change

One record per variant. The four qhi loops are gone, and the three query branches
collapse into one: whether the truth sample is populated or empty is the only thing that
differed between them.
print_var_sample (src/variant.cpp:396) renders GT from orig_gts[vi] and
ploidies[vi]. Deliberately not calc_gts: that is vcfdist's inferred genotype, and
putting it in the GT column would present the inference as the caller's claim. The truth
sample likewise reports its own orig_gt rather than the phase-adjusted haplotype the old
per-haplotype record carried; BS and FE still report the phase relationship between the
two samples, so nothing is lost.
The per-haplotype fields — BD, BC, RD, QD, BK, SG — become comma-separated
lists with one value per GT allele. QQ, SC, PS, PB, BS, VP, FE, GE are
already per-variant and stay Number=1.
Value order is 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, so allele
i reads lane i ^ swap.

Reference alleles report `.`

A GT allele that is reference was never evaluated, so every per-haplotype field reports
. for it rather than the untouched lane's zeroed-out contents. Reporting the lane would
invent a second decision per heterozygous call: a het truth 0|1 that is a TP would render
BD=FN,TP, with the FN describing its reference allele. With ., the number of
TP/FP/FN values in the file is exactly what it was before — on chr20, 177,276 / 3,768 /
7,648 both ways, across 100,207 records before and 72,584 after.

`Number=.`, not `Number=P`

Number=P — one value per GT allele — is what these fields actually carry, but it is a VCF
4.4 addition and BCF_VL_P first appears in htslib 1.23 (absent in 1.19 through 1.22.1), so
consumers on older bcftools or pysam would report a cardinality error. These fields are
informational and downstream pipelines routinely pin older htslib, so compatibility wins.
Number=. produces byte-identical records and gives up only the declared cardinality, which
is why the tests assert the value count directly. Since nothing in the format then conveys
either the count or the order, both are stated in the field descriptions. Number=P is the
right end state once htslib 1.23 is widespread — a header-only change.

Het-alt records stay split

A 1|2 record is parsed into two entries with different ALTs (src/variant.cpp:903-908) and
the cross-haplotype merge cannot rejoin them (src/cluster.cpp:137-150). They remain two
co-located biallelic records, for two reasons:

The two alleles need not share a position after normalization. vcfdist trims and
re-anchors per allele (src/variant.cpp:947-956), so REF=T ALT=G,TT can yield a SNP and
an INS at different normalized positions. One record would mean falling back to the source
record's un-normalized POS/REF/ALT while every other record in the file carries
normalized alleles.
Credit would have to be aggregated across two independently evaluated entries, which may
sit in different superclusters. Unlike the homozygous case, this needs a new cross-entry
rule.

The consequence is a permanent count residual at het-alt sites relative to tools that count a
1|2 record as one location. The new integration test pins this asymmetry explicitly, so it
is deliberate rather than incidental.

`BD`/`BK` are lists here, and only here

Making them multi-valued is the mechanical consequence of un-splitting and combines no
information, so it needs no aggregation policy. #49 then replaces the list with a single
per-site value derived from its match-tier ladder, which is where "which haplotype's decision
wins" belongs.

Verification

pytest: 109 passed, including the 714-test unit suite.
Guardrail: on the chr20 fixture, against dev's tip (3469842), precision-recall-summary.tsv
is byte-identical before and after, as are precision-recall.tsv, phasing-summary.tsv, genotype-errors.tsv,
switchflips.tsv, phase-blocks.tsv, query.tsv, and truth.tsv. Un-splitting the output
records does not change the counts; the counting convention is D4: GA4GH Benchmarking VCF Compatibility #49's change, not this one.
The TP/FP/FN value tallies across summary.vcf are unchanged (above), so no evaluation
value is invented or lost.
WriteSummaryVcf.PerAlleleValuesFollowTheGenotypeSwap was checked against a mutant that
drops the calcgt_is_swapped lookup: it reports FP,. / 0.000000,. / 7,. and fails.
doxygen: no warnings. Build at -Wall -Wextra: no new warnings.

Tests

record_shapes, a new integration fixture whose query and truth are identical so only the
record shape is under test: a hom SNP, a hom CPX (split into an INS and a DEL at parse
time), a hom deletion, and a het-alt 1|2. Pins one record with GT=1|1 and two values per
per-haplotype field for each homozygous call, two co-located records for the het-alt, and
the absence of the pre-change per-haplotype shape.
Unit tests in tests/unit/src/test_phase.cpp for the hom SNP, hom indel, het (reference
allele dotted), haploid (one value), allele order on a homozygous record, the swap, the
het-alt split, and the header cardinality declarations.
Existing summary.vcf assertions updated to the un-split shape; the must_not_contain
guards were widened from ":FP:" to ":FP" so they still bite against a list.

TimD1 · 2026-08-07T03:05:48Z

+ * @return The credit with six digits after the decimal point
+ */
+static std::string credit_str(float credit) {
+    char buf[32];


no static length arrays. Use the C++ stdlib

avoid defining a custom function for this if at all possible

TimD1 · 2026-08-07T03:07:37Z

+void ctgVariants::print_var_sample(FILE* out_fp, int vi, int sc_idx, int phase_block,
+        bool phase_switch, bool phase_flip, bool query /* = false */) {
+
+    // a haploid record carries one bare allele; an unknown ploidy (0) is rendered as diploid


there should be no records with "unknown ploidy". They should all be 1 or 2.

TimD1 · 2026-08-07T03:08:20Z

+
+    bool swap = this->calcgt_is_swapped(vi);
+    std::string errtypes, credits, ref_eds, query_eds, match_types, sync_groups;
+    for (int ai = 0; ai < alleles; ai++) {


This isn't iterating over alleles, it's iterating over haplotypes. I use hi for this elsewhere.

TimD1 · 2026-08-08T02:16:12Z

+    // 4.4's Number=P declares. BCF_VL_P only reaches htslib in 1.23, so a consumer on any older
+    // bcftools or pysam would report a cardinality error; Number=. produces byte-identical records
+    // and merely gives up the declared cardinality, so the count and order are stated here instead.
+    const char* per_allele = " One value per allele of this sample's GT, in GT allele order, "


use a string, not char*

TimD1 · 2026-08-08T02:19:54Z

+        bool phase_switch, bool phase_flip, bool query /* = false */) {
+
+    // a haploid record carries one bare allele; an unknown ploidy (0) is rendered as diploid
+    int alleles = this->ploidies[vi] == 1 ? 1 : HAPS;


please create a new enum for ploidies that is either HAPLOID=1 or DIPLOID=2, following the existing style

TimD1 · 2026-08-08T02:27:28Z

+
+    // a haploid record carries one bare allele; an unknown ploidy (0) is rendered as diploid
+    int alleles = this->ploidies[vi] == 1 ? 1 : HAPS;
+    const std::string gt = alleles == 1 ? "1" : gt_strs[this->orig_gts[vi]];


I think the best thing here is to define a helper function that maps orig_gt, phasing, and ploidy to the display GT string.

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.

#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.

TimD1 · 2026-08-08T03:04:24Z

    int rec_idx = -1;             ///< source VCF record ordinal (0-based, -1 = unknown)
    int alt_idx = -1;             ///< original ALT ordinal (1-based, -1 = unknown)
-    uint8_t ploidy = 0;           ///< variant ploidy from std::abs(ngt) (0 = unknown)
+    ploidy_t ploidy = PLOIDY_DIPLOID; ///< haplotypes the variant was called on, from std::abs(ngt)


this is variant ploidy, not haplotype it was called on

TimD1 · 2026-08-08T03:04:46Z

+        }
+
+        // the evaluation lanes are keyed by matched_gt's haplotypes, not orig_gt's
+        hap_t hi_resolved = swap ? other_hap(hi) : hi;


I think hi_matched is better

TimD1 · 2026-08-08T03:04:58Z

    std::vector<int> rec_idxs;      ///< source VCF record ordinal (0-based, -1 = unknown)
    std::vector<int> alt_idxs;      ///< original ALT ordinal (1-based, -1 = unknown)
-    std::vector<uint8_t> ploidies;  ///< variant ploidy from std::abs(ngt) (0 = unknown)
+    std::vector<ploidy_t> ploidies; ///< haplotypes each variant was called on, from std::abs(ngt)


TimD1 · 2026-08-08T03:06:41Z

 /** @brief Returns the other haplotype of the pair. */
 constexpr hap_t other_hap(hap_t h) { return h == HAP1 ? HAP2 : HAP1; }

+/** @brief How many haplotypes a variant's genotype was called on. */


I'm technically encoding a haploid variant as 0/1, but that's just an implementation detail. It's really 1.
So the explanatory comment should just say "variant ploidy", not "how many haplotypes it's called on", because 0/1 is called on one haplotype but is diploid.

…#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.

Base automatically changed from 102_td_remove-per-contig-ploidy to dev August 6, 2026 20:23

TimD1-bot force-pushed the 103_td_one-entry-per-variant branch from f29bd04 to ac6a9e4 Compare August 6, 2026 20:33

TimD1 force-pushed the 103_td_one-entry-per-variant branch from ac6a9e4 to 222a57c Compare August 8, 2026 02:01

TimD1 requested changes Aug 8, 2026

View reviewed changes

TimD1 added 2 commits August 7, 2026 22:39

TimD1-bot force-pushed the 103_td_one-entry-per-variant branch from 222a57c to a88d34e Compare August 8, 2026 02:51

TimD1 requested changes Aug 8, 2026

View reviewed changes

TimD1 approved these changes Aug 8, 2026

View reviewed changes

TimD1 merged commit d752879 into dev Aug 8, 2026
1 check passed

TimD1 deleted the 103_td_one-entry-per-variant branch August 8, 2026 03:29

Navigation Menu

refactor(phase): emit one summary VCF record per variant (#103) - #208

refactor(phase): emit one summary VCF record per variant (#103)#208
TimD1 merged 3 commits into
devfrom
103_td_one-entry-per-variant

TimD1-bot commented Aug 5, 2026 •

edited

Loading

Uh oh!

TimD1 Aug 7, 2026

Uh oh!

TimD1 Aug 7, 2026

Uh oh!

TimD1 Aug 7, 2026

Uh oh!

TimD1 Aug 7, 2026

Uh oh!

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

TimD1 Aug 8, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

TimD1-bot commented Aug 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Problem

Change

Reference alleles report .

Number=., not Number=P

Het-alt records stay split

BD/BK are lists here, and only here

Verification

Tests

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

TimD1-bot commented Aug 5, 2026 •

edited

Loading

Reference alleles report `.`

`Number=.`, not `Number=P`

`BD`/`BK` are lists here, and only here