← Reproduction studies

Do variant-caller validations survive a benchmark update?

A reproduction study of NIST Genome in a Bottle HG001 releases v3.3.2 and v4.2.1, run as a Facetwork FFL workflow. September 2026.

1. Purpose of the analysis

Genome in a Bottle (GIAB) is the reference truth set used worldwide to validate germline variant callers. A laboratory measures its caller's sensitivity and precision by comparing calls against a GIAB benchmark VCF restricted to GIAB's high-confidence regions.

NIST periodically publishes new benchmark releases. That creates a practical question with real consequences for any lab that has already validated:

If a caller was validated against HG001 v3.3.2, does that validation still hold on v4.2.1 — and if not, exactly where does it stop holding?

Answering it requires knowing not just how many variants changed, but which kind of claim changed: whether the sites themselves moved, or only the genotype evidence attached to them. Those imply different remediation.

2. What the original work does

NIST's Genome in a Bottle consortium does not sequence a genome and publish the calls. It integrates many independent callsets — multiple sequencing technologies, multiple variant callers — and arbitrates between them to produce a benchmark that is trustworthy enough to grade other people's callers against.

That integration is visible in the released VCF's own annotations, which is how a consumer can tell how much evidence stands behind any single call:

The benchmark ships as two artifacts that must be used together: the VCF of high-confidence calls, and a BED of the regions in which the benchmark makes any claim at all. A caller is graded only inside that BED — outside it, the benchmark is silent rather than negative. Method background: Zook et al., Nature Biotechnology 2019.

NIST describe the v4 series as substantially expanding the high-confidence regions over v3 by incorporating additional sequencing technologies and improved integration, and as revising genotype representation. A ~48% larger BED is the published evidence for that expansion.

What we do differently

These are complementary, not competing. We do not re-derive the benchmark and could not: the reads, the callers and the curation are not reproduced here. What this workflow answers is the question a lab actually faces when a new release lands — does my previous validation still hold? — which the release artifacts do not state and which requires comparing them.

3. What our FFL workflow found

The workflow fetches both published releases, compares them keyed on variant identity (CHROM, POS, REF, ALT), and reports the highest level at which they agree. Results below are over the complete GRCh38 benchmarks — 3.6 M and 3.9 M variants — compared in under a minute on a single commodity node.

Two things had to be right before the numbers meant anything, and neither is obvious. A VCF line can carry several ALT alleles at once, so keying on the whole ALT field asks "were these alleles grouped identically?" when the question is "is this allele present in both?"; the comparison therefore splits ALT and compares alleles. And the two releases do not span the same genome — see below.

The releases do not cover the same genome

Before any content can be compared, the two artifacts have to be describing the same key space — and here they are not. v3.3.2 covers chromosomes 1–X; v4.2.1 covers 1–22 and contains no chrX at all. The release filenames say so outright: …CHROM1-X_v.3.3.2… versus …1_22_v4.2.1….

That is 101,253 alleles — 41% of everything a naive comparison calls "retired". They were not retired. They are out of scope, and the distinction is the whole answer for a lab: "your chrX validation is void" and "the new release makes no claim about chrX" demand completely different responses.

Set membership, over the in-scope alleles only:

Without separating scope from content the same comparison reports 247,380 retired variants and 93.2% carryover. Both are wrong, in the same direction, and neither looks wrong: the run is green, the arithmetic is sound, and the error is entirely in what the numbers were taken to mean.

Why every genotype differs

The genotype column differs for every shared variant, in three distinct ways at once. A representative record:

chr1:817186 G>A
  v3.3.2  FORMAT=GT:DP:ADALL:AD:GQ:IGT:IPS:PS   1|1:823:0,381:78,454:283:1/1:.:PATMAT
  v4.2.1  FORMAT=GT:PS:DP:ADALL:AD:GQ           1/1:.:1041:78,459:0,88:348

1. FORMAT restructured — IGT and IPS dropped, PS moved forward.
2. Phasing dropped — 1|1 became 1/1.
3. Evidence values genuinely increased — depth 823 → 1041, genotype quality 283 → 348. This is not reformatting; v4.2.1 rests on more data.

The operational answer

Site-level validation carries over. Genotype-level validation does not. A lab measuring sensitivity/precision on variant sites can rely on its v3.3.2 work for the shared set and needs only to evaluate the 523,745 added alleles — and must note that chrX is no longer covered at all. A lab whose acceptance criteria touch genotype quality, depth, or phasing must re-run in full: every one of those records changed.

4. How our workflow compares to the published work

These are complementary rather than competing artifacts, and the distinction matters:

The algorithmic difference that matters

The conventional comparison for this task is a checksum or a diff, and both give the wrong answer here. A checksum reports failed; a naive diff reports every record as changed. Neither distinguishes "the sites are identical and only the genotype encoding moved" from "the calls disagree" — yet those demand opposite responses.

Our comparison instead reports the highest level at which two outputs agree:

exists → schema → cardinality → keys → values → bytes

and, critically, reports every normalisation required to reach that level. Sorting keys or applying a tolerance until two files match does not demonstrate a reproduction — it manufactures one — so each is declared in the result rather than hidden inside it.

This study also found the rung the ladder was missing. Every level above assumes the two artifacts describe the same key space, and these two did not. A comparison that reports content differences between artifacts covering different chromosomes answers a question nobody asked, in a form indistinguishable from the answer to the one they did — so scope is now checked beneath the ladder rather than assumed by it.

Scale

The comparison indexes one file at a time into a compact digest map — 16 bytes per key — rather than materialising records, so peak memory tracks the key count rather than the file size. Measured:

The naive version needed more memory for a sixth of the data than the streaming one needs for all of it. That is not an optimisation footnote. The first whole-genome attempt exhausted a 23 GiB VM, the runner stopped answering its heartbeat, the reaper reclaimed the task as a dead host, and it dead-lettered after five retries and 24 minutes — a failure that reads as an infrastructure fault and was arithmetic. A task-level heartbeat cannot rescue a process that is thrashing; the memory had to come down.

5. Reproducing this

The workflow is 31 lines of FFL (plus 70 of commentary) composed entirely from built-in facets — no domain-specific handlers were written for this. fw.http.Fetch acquires each published release and records its digest and provenance; fw.compare.Tabular compares and writes the report.

An earlier version also called fw.compare.Summarise, which aggregates many file-pair verdicts into the weakest one. Over a single pair it was handed an empty verdict list and correctly answered "nothing was compared" — the guard working, not a result. It was removed rather than left in to look thorough.

fw ffl run --primary examples/repro-giab/giab_benchmark_diff.ffl \
           --workflow repro.giab.CompareBenchmarkReleases

A page about reproducibility should say whether it reproduces. It was run twice, an hour apart, on different runner containers of a distributed fleet:

Byte-identical. That is a weaker claim than it sounds — the inputs are fixed artifacts and the comparison is deterministic, so anything else would have been a bug. It is worth stating precisely because the interesting failures in this project have been the ones that looked fine: an orchestrator that retries a task can leave a duplicate output record and make a perfectly reproducible pipeline read as irreproducible, which is why the comparison keys on record identity rather than file position.

6. Exactly what was compared

A reproduction study that cannot be checked is an assertion. fw.http.Fetch writes a provenance sidecar beside every download, so these are the precise bytes behind every number on this page — not "the v3.3.2 release" but one identified artifact. Anyone can fetch the same URLs and verify the digests before deciding whether our figures are worth arguing with.

Both are served by NIST with a Last-Modified older than four years, so the sources are stable and a later rerun compares the same artifacts rather than a moving target — which is itself the failure mode that froze an earlier pipeline here on a -latest URL.

7. Limitations

• This page has been wrong twice, in the same way both times — a number that was arithmetically correct and meant something other than what it was labelled. First it reported figures from a partial slice as though they were whole-genome (the slice was chromosome-1 heavy and under-counted added sites threefold). Then it reported 242,958 variants as retired when 41% of them were simply outside v4.2.1's scope. Both drafts ran green. The corrections are recorded here rather than quietly swapped in, because anyone who acted on the earlier figures needs to know they moved.
• Multi-allelic records are split so alleles are compared individually. Trimming to parsimony and MNP decomposition are not applied, and left-alignment through repeat regions is not either — that one needs the 3 GB GRCh38 reference for random access. This is a measured decision, not an oversight: ablated over the full pair, splitting ALT rescues 396 records from "wholly unmatched", adding trimming rescues 24 more, and MNP decomposition rescues none. Zero records in either release are even trimmable, because both come from pipelines that already normalise. Against an arbitrary variant caller's output the calculation would come out differently, and these are the numbers to re-measure before assuming otherwise.
• Scope detection uses the first key column (here CHROM) and only while its cardinality stays modest. It catches a whole partition missing from one side; it does not catch a partially-covered one. The high-confidence BED files define coverage precisely and are not consulted — a variant inside both VCFs but outside the newer BED region is still counted as shared here.
• GIAB's methodology is not re-executed, so this cannot detect an error in the benchmark construction itself — only differences between two published outputs.