v4.5.0: add neat compare-vcfs subcommand (hap.py wrapper with NEAT-aware FN attribution)

[joshfactorial]

Motivation

NEAT 2.1 shipped utilities/vcf_compare_OLD.py — a 721-line Python 2 truth-vs-test VCF comparator. The filename suggests the original authors already considered it deprecated. Its core feature (±50 bp splice-equivalence check for FP/FN pairs) is strictly weaker than what rtg vcfeval and Illumina hap.py provide via haplotype-based equivalence, both of which are free, BSD-licensed, GA4GH-standard, and widely used.

Reimplementing that tool in pysam-based Python 3 would be a rewrite, not a port, and would inherit the entire long tail of VCF edge cases (multi-allelic sites, gVCF, phased GT, SVs) on a tool weaker than the off-the-shelf alternative.

The NEAT-specific value is FN reasoning — explaining misses in terms of the simulation's own parameters: "the FN landed in a region the mutation_bed didn't cover" or "the FN was at a position with delivered coverage below the configured threshold." hap.py cannot do that. A thin wrapper that runs hap.py and then cross-references its FN VCF against the NEAT run config can.

Scope (v4.5.0)

CLI

neat compare-vcfs <golden.vcf> <called.vcf> --neat-run-dir <dir> --output-dir <dir>
                  [--reference <fa>] [--target-bed <bed>] [--happy-bin <path>]

• golden.vcf — the NEAT-simulated truth VCF
• called.vcf — the workflow VCF (output of a variant caller on the simulated reads)
• --neat-run-dir — the directory containing the NEAT simulation output. The wrapper reads simulation_summary.json (new — see below) from here for cross-referencing.
• --happy-bin — optional override; defaults to looking for hap.py on $PATH. If not found, error with installation instructions (conda / docker).

Implementation steps

1. New CLI command at neat/cli/commands/compare_vcfs.py — mirrors the structure of existing subcommands (e.g., gen_mut_model.py).
2. hap.py invocation as a subprocess. Validate it's on $PATH; fail with a clear install hint if not. Pass through --reference and --target-bed directly.
3. Output parser that reads hap.py's TP/FN/FP VCFs (it emits these with INFO tags). Use pysam.
4. simulation_summary.json emission from the read simulator — new minimal artifact. Just the config echo (coverage, read_len, paired_ended, gc_model path, mutation_bed path, target_bed path, threads) plus delivered counts (total reads, total variants, per-contig variant counts). Keep it tight; full per-region telemetry is out of scope for v4.5.0.
5. NEAT-aware FN attribution. For each FN in hap.py's output, tag it with one or more reasons:
   • outside_mutation_bed — position is outside the regions the simulator was told to mutate
   • outside_target_bed — position is outside the target capture regions
   • outside_simulated_contigs — position is on a contig that wasn't in the NEAT run
   • unknown — no NEAT-specific reason found
     The reason set is intentionally narrow for v4.5.0; richer attribution (coverage troughs, GC-bias troughs, hotspot misses) is a follow-up.
6. Report output:
   • comparison_summary.json — precision/recall/F1 from hap.py + per-reason FN counts
   • comparison_summary.txt — human-readable rollup
   • FN_with_reasons.vcf — hap.py's FN VCF with an added NEAT_REASON INFO tag
   • Optional venn plot (--plot flag) of FN reasons; matplotlib is already a dep, so no extra burden

Testing

• Unit tests for the FN-attribution logic (synthetic FN VCF + synthetic simulation_summary.json, assert the right NEAT_REASON tags).
• One end-to-end integration test that runs a small NEAT simulation, runs hap.py against a synthetic "caller" VCF (just the golden VCF with a few variants dropped/added), and checks the report shape. Gated on hap.py being installed; skip cleanly if not.

Documentation

• README section under "Utilities" describing the tool and the install path for hap.py.
• ChangeLog entry under v4.5.0.

Explicitly out of scope

• Reimplementing variant equivalence. Delegate entirely to hap.py/vcfeval.
• Full per-region telemetry instrumentation (per-chunk delivered coverage, GC-bias map per contig, error rates by position). Tracked separately; needed before richer FN attribution.
• SV comparison. hap.py handles small variants; SV truth-vs-test is a different problem and a different tool ecosystem.
• Somatic-specific reporting. Use som.py directly if needed.

Dependencies

• hap.py is a Python 2 / Docker-distributed tool. Modern installations use conda (bioconda::hap.py) or its Docker image. We declare it as an optional runtime dependency, document the install path, and fail cleanly when absent.
• No new Python package dependencies for NEAT itself (pysam and matplotlib are already deps).

Acceptance criteria

• neat compare-vcfs --help documents the subcommand
• Running the subcommand on a NEAT-simulated truth VCF + a synthetic "caller" VCF produces comparison_summary.{json,txt} and FN_with_reasons.vcf
• Each FN has a NEAT_REASON tag matching one of the four reason categories
• Unit tests cover the FN-attribution branches
• Tool degrades gracefully when hap.py isn't installed (clear error, exit 2)
• simulation_summary.json is emitted by every neat read-simulator run