feature/vcftools-benchmark

[bwalsh]

Add KING-robust kinship matching and VCFtools benchmarks

Summary and use case

Enable a bioinformaticist to compare two observations or search an indexed cohort for potential duplicates and relatives using Manichaikul (KING) kinship estimates. The existing carried-allele index supports allele similarity but cannot distinguish a called reference genotype from missing or filtered data. This branch adds explicit genotype indexing and a king-robust plugin while preserving the existing identity workflow.

It also gives engineers a reproducible way to compare numerical results and computational costs with VCFtools --relatedness2, separating initial indexing from repeated queries. See the KING user story and acceptance criteria and bioinformatics comparison with VCFtools.

Intended base: development (9382c73).

Architecture changes

• Ingestion and storage: add load-samples --index-genotypes --panel PANEL.tsv. Additional SQLite tables store the declared panel, marker identities, observation provenance/completion status, and passing ALT dosages 0, 1, and 2. Validate reference/annotation consistency, observation uniqueness, and QC policy; commit genotype and allele indexing atomically. An absent genotype row means unknown, never reference. The carried-allele table retains its existing identity semantics.
• Plugin API and results: extend PluginContext with panel metadata and explicit genotype access. Add dedicated KinshipResult and KinshipMatches types, CLI rendering, and full-precision JSON. Score jointly called markers, preserve negative estimates, and report zero-denominator cases as unscorable. Existing identity-style plugins retain API version 1 and their result types. See the plugin contract.
• Matching execution: batch SQL retrieval, process all-pairs comparisons in 128-sample tiles, and stream distinct pair results. One-versus-all searches use a bounded heap for top-N selection, with deterministic ties and separate unscorable results. SQLite remains the durable index; packed representations and parallel scoring are deferred. See the performance architecture ADR.
• Benchmarking: add an opt-in runner for prepared data or deterministic synthetic cohorts. Measure index construction, single-query latency, indexed all-pairs matching, and direct VCFtools execution separately. Instrument internal retrieval and scoring phases; record workload dimensions, hashes, timings, peak RSS, storage, and comparator results. Labeled evaluations produce relationship/ancestry/callable-overlap summaries and family-block bootstrap intervals for duplicate retrieval. The benchmark ADR distinguishes the original genotype-discordance proposal from this implemented KING track.
• Runtime: raise the supported Python minimum to 3.13, updating the environment, lockfile, CI workflows, and setup documentation. No new runtime dependencies are added.

Usage and compatibility

vrs-matcher load-samples cohort.vrs.vcf.gz --db cohort.db \
  --index-genotypes --panel king-panel.tsv
vrs-matcher match-samples SAMPLE_A SAMPLE_B --db cohort.db \
  --algorithm king-robust --json
vrs-matcher match-sample SAMPLE_A --db cohort.db \
  --algorithm king-robust --top 10 --json

The KING usage guide documents panel/header requirements, filtering, output fields, and benchmark commands. The initial scope is human autosomal, diploid, biallelic SNVs. Existing allele-only databases still support identity; KING requires re-ingestion from source VCFs because reference calls cannot be reconstructed. shared-variants rejects kinship results.

The default between-family score is distinct from the within-family diagnostic used for VCFtools comparison. On the eight-marker oracle, these are 0.0 and 0.1, respectively; VCFtools reports 0.1. Missing-data tests explicitly verify the difference between the plugin's pairwise callable denominator and VCFtools' individual heterozygote totals.

Validation

Verified at branch head with Python 3.13.5:

• ruff check . and ruff format --check . — passed.
• pytest -q --no-cov — 137 passed, 3 opt-in tests skipped.
• Explicitly enabled VCFtools integration tests — 2 passed, using VCFtools 0.1.16.

Tests cover VCF ingestion, QC boundaries, missing/reference/partial calls, panel validation, rollback, statistical oracles, candidate restrictions, ranking, batching, streaming, CLI/JSON compatibility, synthetic generation, and phase accounting. The existing 1000 Genomes network/SeqRepo integration test was not run. The validation report records the earlier numerical validation and smoke measurements; the performance ADR records subsequent instrumentation experiments.

Limits

Biological acceptance remains pending a genome-wide panel with independently verified replicate and pedigree labels. Synthetic results establish implementation correctness and instrumentation, not biological accuracy or representative throughput. Timings exclude upstream preparation and VRS annotation. VCFtools emits ordered pairs and diagonals, while the plugin all-pairs path emits distinct unordered pairs; reported workload ratios do not establish universal performance superiority. Unscorable search results remain retained in memory.

[bwalsh]

Performance comparison with vcftools 🐌

export GA4GH_VRS_DATAPROXY_URI=seqrepo+file://$HOME/.local/share/seqrepo/2024-12-20
export VCFTOOLS=/opt/homebrew/bin/vcftools

uv run python scripts/benchmark_king.py   --synthetic-samples 128 --synthetic-markers 512 --seed 0   --repeats 3 --vcftools $VCFTOOLS   --output /tmp/king-synthetic-128x512 --smoke

---

KING benchmark report

Mode: smoke. Cohort: 128 observations, 512 markers, 8128 unique pairs.

Input preparation and VRS annotation are excluded; this is not end-to-end timing.
VCFtools computes both pair directions and diagonals; the plugin all-pairs stage
computes each unordered distinct pair once. Query is a separate workload.

Workloads are reported separately: load builds the genotype and allele
SQLite indexes; query matches one sample against the cohort; all matches
all unordered distinct pairs using an existing index; vcftools reads the
prepared VCF and computes --relatedness2. First-run time is repeat 0.
Subsequent median and range exclude repeat 0. Peak RSS is the maximum process
peak across all repeats for that stage, in MiB (1024² bytes). All times
include process startup.

Stage	First run (s)	Subsequent median (s)	Subsequent range (s)	Max peak RSS (MiB)
load	0.830728	0.797891	0.796183–0.799599	45.23
query	0.314482	0.352697	0.326327–0.379067	46.00
all	1.088717	1.143111	1.142282–1.143939	46.44
vcftools	0.093353	0.089595	0.082784–0.096406	2.45

Isolated index, retrieval, and scoring phases

These phase durations are measured inside the worker, excluding subprocess startup. index build times input parsing and database creation. call retrieval times SQLite reads into tiled genotype maps; scoring + iteration covers pair enumeration, score calculation, and result construction but excludes those retrieval calls. Process peak RSS is reported for the worker that measured each phase.

Phase	First run (s)	Subsequent median (s)	Subsequent range (s)	Max worker peak RSS (MiB)
Index build	0.635790	0.631461	0.631310–0.631611	45.23
Call retrieval	0.111596	0.113252	0.112496–0.114009	45.81
Scoring + pair iteration	0.736524	0.747765	0.741924–0.753607	45.81

Workload comparisons

Relative performance is 100 × VCFtools median / VRS-Matcher median; VCFtools is 100%. Each VRS-Matcher median is computed from paired repeats.

VRS-Matcher workload	VRS-Matcher median (s)	VCFtools direct median (s)	Relative performance (% of VCFtools speed)
Indexed all-pairs matching only	1.143111	0.089595	7.8
Index build + all-pairs matching	1.941002	0.089595	4.6

The first comparison measures matching against an already-built VRS-Matcher index; the VCFtools figure still includes reading the VCF, so it is not an algorithm-only comparison. The second includes VRS-Matcher index construction and is the one-shot pipeline view. VCFtools reads the prepared VCF directly. The tools also emit different pair sets as noted above. The query row is a separate one-versus-all workload and has no equivalent single-query VCFtools operation in this benchmark.

Identity-only DB bytes: 9519104.
Genotype + allele DB bytes: 17747968.
Annotated VCF bytes: 701710.

See manifest.json, raw logs, pairs.tsv, and timings.tsv for reproducibility.
No biological accuracy claim is made by a smoke run.

[bwalsh]

Title: feature/vcftools-benchmark

Purpose: Add KING-robust kinship matching and VCFtools benchmarks to enable bioinformaticists to estimate relationships between samples and search indexed cohorts for duplicates and relatives using KING (Manichaikul) kinship estimates. This PR introduces a second matching algorithm alongside the existing identity plugin, an opt-in genotype indexing path, and comprehensive benchmarking capabilities.

Key Changes:

• 3,520 lines added, 128 deleted across 33 files
• 5 commits from author bwalsh (created 41 minutes ago)
• 1 comment on the PR
• Status: Open and mergeable into development branch

Architecture & Implementation

Core additions:

1. Genotype-enabled indexing (genotypes.py, loader.py): New --index-genotypes --panel PANEL.tsv workflow that stores explicit called genotypes (dosage 0, 1, 2) atomically, alongside existing identity-based carried-allele indexing. Panel validation ensures consistent reference/annotation provenance.

2. KING-robust plugin (king.py): Manichaikul equations 11 (between-family) and 9 (within-family) implemented on jointly called markers. Returns dedicated KinshipResult and KinshipMatches types with full-precision JSON. Handles missing calls explicitly and reports unscorable pairs separately.

3. Database schema (db.py): New tables genotype_panel, genotype_marker, genotype_observation, genotype_call to store panel metadata, markers, and per-sample call dosages. Reference calls stored as dosage 0; missing/excluded calls have no row (unknown vs. reference distinction).

4. CLI & matching API extensions (cli.py, matcher.py): Added --index-genotypes, --panel, --json flags. Result rendering for KING includes full precision and explicit status/reason fields. Deprecated the generic JSON output in favor of dedicated kinship result types.

5. Batched retrieval & streaming (king.py, plugins.py): All-pairs matching processes 128-sample tiles to bound memory. Top-N selection uses a bounded heap with deterministic tie-breaking. Streaming pair results to TSV in the benchmark runner.

6. Comprehensive documentation: 7 new docs files covering KING usage (how-to-king-robust.md), benchmark ADR (adr-benchmark-vcftools.md), performance architecture ADR (adr-king-performance.md), index structure (index-description.md), user story with acceptance criteria, validation report, and VCFtools comparison rationale.

7. Benchmark runner (scripts/benchmark_king.py): Opt-in, offline tool separating index build, one-versus-all query, all-pairs matching, and direct VCFtools execution. Supports deterministic synthetic input generation, labeled evaluation, and bootstrap intervals for retrieval accuracy. Measures time, peak RSS, and workload comparisons.

8. Python 3.13 migration: Minimum version raised from 3.12. All CI workflows, pyproject.toml, and documentation updated.

Test Coverage

• 137 passed, 3 opt-in tests skipped in the main suite
• 2 passed in the opt-in VCFtools integration test (tests/integration/test_king_vcftools.py)
• 528 lines of comprehensive acceptance tests (tests/test_king.py) covering:
  • A: Indexing, integrity, QC masking, exclusions, rollback
  • B: Hand-calculated oracle, zero denominators, phase invariance, restrictions
  • C: Plugin integration, ranking, CLI/JSON, batching
  • D: External VCFtools 0.1.16 comparison with exact count validation
  • E: Synthetic benchmark (smoke run recorded in docs/king-validation.md)

Validation & Reproducibility

• Ruff lint and format checks passed
• Smoke benchmark with 8-marker synthetic oracle, 2 observations, 6 repetitions documented
• Validation report records numerical correctness, VCFtools comparison, and synthetic measurements
• Input hashes, manifest, and raw logs included for reproducibility
• Biological validation deferred pending genome-wide panel with verified pedigree/replicate labels

---

Merge Readiness and Risk Assessment

Status: Clean and ready to merge with attention to the notes below.

Strengths:

• Comprehensive scope executed with explicit acceptance criteria and validation
• New KING functionality cleanly separated from existing identity workflow
• Backward compatible: existing allele-only databases still support identity; KING is opt-in
• Detailed documentation and design decisions (ADRs) explain rationale and caveats
• All new runtime dependencies avoided; Python 3.13 is a version bump, not a new dep
• Atomic ingestion with rollback ensures database integrity
• Clear re-ingestion error messages guide users when KING requires new indexing

Observations:

• This is a substantial addition (3.6k LOC); focus review on the plugin contract, batching logic, and database invariants
• Benchmark runner is comprehensive but intentionally defer biological accuracy claims until real-world replication study
• Within-family KING estimate differs from VCFtools due to documented missing-data semantics; this is explicit and tested
• Reference calls now stored in the database; existing allele-only observations cannot be retroactively KING-enabled (rebuild required)

• Verify the PluginContext API: The extension with get_called_genotypes() and get_called_genotypes_many() is new; ensure allele-only samples correctly raise re-ingestion errors and that the batched retrieval SQL is correct for large cohorts.
• Confirm database foreign-key enforcement: The genotype tables use FKs to samples and genotype_marker. Verify that open_db() enables PRAGMA foreign_keys consistently and that the genotype-ingestion transaction properly validates before insertion.
• Review the benchmark's workload isolation: Separate timing for index build, retrieval, and scoring is intentional; ensure the report clearly labels each workload and does not claim universal speed superiority based on the smoke fixture.
• Check CLI output for KinshipMatches when top_n is used: The unscorable list should remain separate even when top-N is requested; spot-check that the CLI renders both lists correctly.

---