Merge pull request #299 from ncsa/compare-vcfs

Compare vcfs

Author: joshfactorial

ChangeLog.md

@@ -1,3 +1,48 @@
+# NEAT v4.5.0
+
+New `neat compare-vcfs` subcommand: compares a downstream variant caller's VCF
+against the NEAT-simulated truth VCF and attributes each false negative to the
+simulator's own configuration (mutation bed, target bed, simulated contigs).
+Issue #297.
+
+Variant equivalence (multi-allelic normalization, haplotype-level matching) is
+delegated to Illumina's `hap.py`. NEAT adds the false-negative attribution
+layer: each FN is tagged with one or more of `outside_simulated_contigs`,
+`outside_mutation_bed`, `outside_target_bed`, or `unknown`.
+
+**Outputs (in `--output-dir`):**
+
+- `comparison_summary.json` — counts, precision/recall/F1, per-reason FN totals
+- `comparison_summary.txt` — human-readable rollup
+- `FN_with_reasons.vcf` — hap.py's FN records with an added `NEAT_REASON` INFO tag
+- `happy.vcf.gz` and siblings — preserved hap.py output
+- `fn_attribution.png` — optional bar chart, written when `--plot` is set
+
+**Prerequisite artifact:** Every `neat read-simulator` run now writes a small
+`simulation_summary.json` alongside its other outputs, capturing the run's
+config echo (coverage, read length, paired-ended, BED paths, contigs simulated)
+and delivered counts (total reads, total variants, per-contig variants). The
+`compare-vcfs` wrapper reads this file to drive attribution.
+
+**External dependency:** `hap.py` is required at runtime. NEAT does not bundle
+it; install via `conda create -n hap_py_env -c bioconda -c conda-forge hap.py`
+and pass the absolute path via `--happy-bin`, or put it on `$PATH`. Without
+hap.py, the command exits with an install hint.
+
+**Chromosome-name handling:** `compare-vcfs` detects when a BED's chrom names
+don't overlap the reference's (e.g., BED uses `1`/`MT` while reference uses
+`chr1`/`chrM`) and writes a warning into `comparison_summary.json` suggesting
+an alias mapping. Users can apply the mapping via a new `--chrom-aliases TSV`
+flag. NEAT does not auto-normalize — silent renaming would mask real bugs.
+`simulation_summary.json` now also records `delivered.reference_contigs` (the
+full FASTA contig set) alongside `contigs_simulated`.
+
+**Not in this release** (deferred to follow-up issues): full per-region
+simulation telemetry (per-chunk coverage, GC-bias map, error rates by position)
+for richer FN attribution, and SV-comparison support. The simulator's silent
+"skip BED chroms not in reference" warning will be promoted to a fatal error
+in a future release with an opt-in `chrom_aliases` config key.
+
 # NEAT v4.4.4
 Follow-up release on top of v4.4.3 bundling three lines of work: another perf
 pass over the remaining single-thread hot paths (variant overlap checks,

README.md

@@ -47,6 +47,7 @@ To cite this work, please use both of the following:
     * [`neat model-seq-err`](#neat-model-seq-err)
     * [`neat model-qual-score`](#neat-model-qual-score)
     * [`neat model-gc-bias`](#neat-model-gc-bias)
+    * [`neat compare-vcfs`](#neat-compare-vcfs)
   * [Tests](#tests)
     * [Guide to run locally](#guide-to-run-locally)
     * [Note on Sensitive Patient Data](#note-on-sensitive-patient-data)
@@ -598,6 +599,75 @@ neat model-gc-bias    \
 
 and creates `gc_model.pickle.gz` model in the working directory.
 
+### `neat compare-vcfs`
+
+Compares a downstream variant caller's VCF against a NEAT-simulated truth VCF,
+attributing each false negative to the simulator's own configuration. The
+intended workflow is: run NEAT to produce a synthetic FASTQ and golden VCF →
+run your aligner + caller on the FASTQ → compare the caller's VCF against the
+golden with `neat compare-vcfs`.
+
+Variant equivalence (multi-allelic normalization, haplotype-level matching) is
+delegated to `hap.py`. The NEAT-specific value-add is the false-negative
+attribution: each FN is tagged with one or more reasons drawn from the run
+config recorded in `simulation_summary.json`.
+
+```bash
+neat compare-vcfs golden.vcf called.vcf            \
+        --neat-run-dir /path/to/neat/output/dir    \
+        --output-dir /path/to/comparison/output    \
+        --reference reference.fa                   \
+        [--target-bed target.bed]                  \
+        [--happy-bin /abs/path/to/hap.py]          \
+        [--plot]                                   \
+        [--chrom-aliases aliases.tsv]
+```
+
+Outputs (in `--output-dir`):
+
+| File | Purpose |
+|------|---------|
+| `comparison_summary.json`   | Machine-readable: counts (TP/FN/FP), precision/recall/F1, per-reason FN counts |
+| `comparison_summary.txt`    | Human-readable rollup of the same |
+| `FN_with_reasons.vcf`       | hap.py's false-negative records, each annotated with a `NEAT_REASON` INFO tag |
+| `happy.vcf.gz` (+ siblings) | Raw hap.py output preserved for inspection |
+| `fn_attribution.png`        | Optional — only written when `--plot` is set |
+
+**Chromosome-name conventions:** NEAT does NOT silently normalize chrom names
+across the reference, golden VCF, called VCF, and BED files. If your
+`mutation_bed` uses `1`, `2`, … but the reference uses `chr1`, `chr2`, …,
+`compare-vcfs` detects the mismatch up front and writes a warning into
+`comparison_summary.json` suggesting the rename. Pass `--chrom-aliases
+aliases.tsv` (two-column TSV: `bed_name<TAB>canonical_name`) to apply the
+rename at load time. This also handles common mitochondrial variants
+(`M`/`MT`/`chrM`/`chrMT`).
+
+**False-negative reason categories:**
+
+| Tag | Meaning |
+|-----|---------|
+| `outside_simulated_contigs` | FN's chromosome wasn't part of the NEAT run at all |
+| `outside_mutation_bed`      | `mutation_bed` was set in the config and the FN falls outside its regions |
+| `outside_target_bed`        | `target_bed` was set and the FN falls outside its regions |
+| `unknown`                   | None of the above — NEAT has no specific explanation |
+
+**`simulation_summary.json` prerequisite:** Every `neat read-simulator` run now
+emits `simulation_summary.json` into its output dir as part of the standard
+artifacts; `compare-vcfs` reads this file from `--neat-run-dir` to drive
+attribution. No extra step required when running NEAT yourself; if you're
+working with pre-NEAT-4.5 output, re-run the simulator to produce one.
+
+**`hap.py` install:** `hap.py` is an external dependency. The cleanest install
+path is a dedicated conda env:
+
+```bash
+conda create -n hap_py_env -c bioconda -c conda-forge hap.py -y
+```
+
+then point `--happy-bin` at `/path/to/hap_py_env/bin/hap.py` (or put that
+directory on `$PATH`). Without `hap.py` available, `neat compare-vcfs` exits
+with a clear install hint.
+
 ## Tests
 
 We provide unit tests (e.g., mutation and sequencing error models) and basic integration tests for the CLI.

neat/cli/commands/compare_vcfs.py

@@ -0,0 +1,87 @@
+"""
+Command line interface for NEAT's compare-vcfs subcommand (issue #297).
+"""
+import argparse
+
+from ...compare_vcfs import compare_vcfs_runner
+from .base import BaseCommand
+
+
+class Command(BaseCommand):
+    """
+    Compare a downstream variant-caller VCF against a NEAT-simulated truth VCF
+    and attribute false negatives to the simulator's configuration.
+    """
+    name = "compare-vcfs"
+    description = (
+        "Compare a NEAT-simulated truth VCF (golden) against a downstream variant "
+        "caller VCF (called), with NEAT-aware false-negative attribution."
+    )
+
+    def add_arguments(self, parser: argparse.ArgumentParser):
+        parser.add_argument(
+            "golden_vcf",
+            type=str, metavar="golden.vcf",
+            help="NEAT-simulated truth VCF (typically <prefix>_golden.vcf.gz)."
+        )
+        parser.add_argument(
+            "called_vcf",
+            type=str, metavar="called.vcf",
+            help="Downstream variant caller's VCF produced from the simulated reads."
+        )
+        parser.add_argument(
+            "--neat-run-dir",
+            dest="neat_run_dir",
+            type=str, required=True, metavar="DIR",
+            help="Directory containing the NEAT simulator output, including simulation_summary.json."
+        )
+        parser.add_argument(
+            "--output-dir",
+            dest="output_dir",
+            type=str, required=True, metavar="DIR",
+            help="Where to write comparison_summary.{json,txt} and FN_with_reasons.vcf. Created if absent."
+        )
+        parser.add_argument(
+            "--reference",
+            type=str, default=None, metavar="ref.fa",
+            help="Optional reference FASTA forwarded to hap.py."
+        )
+        parser.add_argument(
+            "--target-bed",
+            dest="target_bed",
+            type=str, default=None, metavar="target.bed",
+            help="Optional BED restricting comparison to these regions."
+        )
+        parser.add_argument(
+            "--happy-bin",
+            dest="happy_bin",
+            type=str, default=None, metavar="PATH",
+            help="Explicit path to the hap.py binary. Defaults to looking on $PATH."
+        )
+        parser.add_argument(
+            "--plot",
+            dest="plot",
+            action="store_true",
+            help="Also write fn_attribution.png — a bar chart of FN reason counts."
+        )
+        parser.add_argument(
+            "--chrom-aliases",
+            dest="chrom_aliases",
+            type=str, default=None, metavar="TSV",
+            help="Two-column TSV mapping BED chrom names to reference-canonical names, "
+                 "applied to mutation_bed and target_bed at load time. Used when the BED "
+                 "uses '1' but the reference uses 'chr1', or similar prefix/mt mismatches."
+        )
+
+    def execute(self, arguments: argparse.Namespace):
+        compare_vcfs_runner(
+            golden_vcf=arguments.golden_vcf,
+            called_vcf=arguments.called_vcf,
+            neat_run_dir=arguments.neat_run_dir,
+            output_dir=arguments.output_dir,
+            reference=arguments.reference,
+            target_bed=arguments.target_bed,
+            happy_bin=arguments.happy_bin,
+            plot=arguments.plot,
+            chrom_aliases=arguments.chrom_aliases,
+        )

neat/common/__init__.py

@@ -3,3 +3,4 @@
 from .io import *
 from .constants_and_defaults import *
 from .ploid_functions import *
+from .chrom_names import *

neat/common/chrom_names.py

@@ -0,0 +1,106 @@
+"""
+Chromosome-name normalization utilities for cross-checking reference/VCF/BED
+conventions.
+
+Genomics tooling ships with multiple incompatible naming conventions:
+  - UCSC:    `chr1`, `chr2`, ..., `chrX`, `chrM`
+  - Ensembl: `1`,    `2`,    ..., `X`,    `MT`
+  - NCBI accession: `NC_000001.11`, ...
+  - Custom:  whatever the lab named the scaffolds
+
+A user can run NEAT with a reference using one convention but supply a
+`mutation_bed` / `target_bed` using another. NEAT does not silently normalize;
+instead, callers detect a mismatch up front and either warn (see
+`compare_vcfs`) or apply an explicit user-supplied alias map.
+
+This module is the single source of truth for chrom-name heuristics so the
+simulator side and the compare-vcfs side stay consistent.
+"""
+import logging
+from pathlib import Path
+
+__all__ = [
+    "apply_aliases",
+    "find_aliases",
+    "load_chrom_aliases",
+    "prefix_flip_candidates",
+]
+
+_LOG = logging.getLogger(__name__)
+
+# Common mitochondrial name variants. A reference using any of these
+# represents the same physical chromosome.
+_MITO_VARIANTS = frozenset({"M", "MT", "chrM", "chrMT"})
+
+
+def prefix_flip_candidates(name: str) -> set[str]:
+    """
+    Return the set of alternative names that probably refer to the same
+    contig as `name`. Excludes `name` itself.
+
+    Two heuristics:
+      1. `chr` prefix add/strip — `chr1` ↔ `1`, `chrX` ↔ `X`.
+      2. Mitochondrial variants — names in {M, MT, chrM, chrMT} are all
+         treated as aliases of each other.
+
+    No other heuristic — accession-style names and custom scaffold names
+    require an explicit user-supplied alias map.
+    """
+    candidates: set[str] = set()
+    if name.startswith("chr"):
+        candidates.add(name[3:])
+    else:
+        candidates.add(f"chr{name}")
+    if name in _MITO_VARIANTS:
+        candidates |= _MITO_VARIANTS
+    candidates.discard(name)
+    return candidates
+
+
+def find_aliases(source: set[str] | frozenset[str], target: set[str] | frozenset[str]) -> dict[str, str]:
+    """
+    For each name in `source` that doesn't appear natively in `target`, propose
+    a canonical form from `target` if one of the heuristics matches. Returns a
+    `{source_name: target_name}` dict; names with a native match or no inferrable
+    candidate are omitted.
+    """
+    mapping: dict[str, str] = {}
+    for name in source:
+        if name in target:
+            continue
+        for candidate in prefix_flip_candidates(name):
+            if candidate in target:
+                mapping[name] = candidate
+                break
+    return mapping
+
+
+def load_chrom_aliases(path: Path | str | None) -> dict[str, str]:
+    """
+    Parse a two-column alias TSV: `source_name<TAB>canonical_name`. Lines
+    starting with `#` and empty lines are skipped. Returns `{}` when `path`
+    is None.
+
+    Whitespace tolerated as a fallback separator for human-typed files.
+    """
+    if path is None:
+        return {}
+    aliases: dict[str, str] = {}
+    with open(path) as fh:
+        for lineno, raw in enumerate(fh, start=1):
+            line = raw.strip()
+            if not line or line.startswith("#"):
+                continue
+            parts = line.split("\t")
+            if len(parts) < 2:
+                parts = line.split()
+            if len(parts) < 2:
+                _LOG.warning(f"{path}:{lineno}: skipping malformed alias line: {raw!r}")
+                continue
+            aliases[parts[0]] = parts[1]
+    return aliases
+
+
+def apply_aliases(name: str, aliases: dict[str, str]) -> str:
+    """Return the canonical form of `name`, or `name` unchanged if no alias maps it."""
+    return aliases.get(name, name)

neat/compare_vcfs/__init__.py

@@ -0,0 +1,2 @@
+"""Load the compare-vcfs runner so the CLI can import it from the package root."""
+from .runner import *