From 46dbd49af358f98505b71cf3975ed73d25b4086f Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 18:04:21 -0500 Subject: [PATCH 1/4] Add neat compare-vcfs subcommand; bump to v4.5.0 (#297) Compares a downstream variant caller's VCF against a NEAT-simulated truth VCF and attributes each false negative to the simulator's configuration. Variant equivalence is delegated to hap.py; NEAT contributes the FN attribution against the mutation/target BED and simulated-contig set recorded in a new simulation_summary.json artifact. Outputs in --output-dir: - comparison_summary.json: counts, precision/recall/F1, FN reasons - comparison_summary.txt: human-readable rollup - FN_with_reasons.vcf: hap.py FNs with NEAT_REASON INFO tag - fn_attribution.png: optional bar chart (--plot) The read simulator now emits simulation_summary.json alongside its other outputs (config echo + delivered counts) as a prerequisite for compare-vcfs. hap.py is an optional external dependency (conda bioconda::hap.py); compare-vcfs exits with a clear install hint when not found. Co-Authored-By: Claude Opus 4.7 (1M context) --- ChangeLog.md | 35 ++ README.md | 60 +++ neat/cli/commands/compare_vcfs.py | 78 ++++ neat/compare_vcfs/__init__.py | 2 + neat/compare_vcfs/attribution.py | 131 +++++++ neat/compare_vcfs/happy.py | 138 +++++++ neat/compare_vcfs/reports.py | 244 ++++++++++++ neat/compare_vcfs/runner.py | 203 ++++++++++ neat/read_simulator/runner.py | 10 + .../utils/simulation_summary.py | 144 ++++++++ pyproject.toml | 2 +- tests/test_compare_vcfs/__init__.py | 0 tests/test_compare_vcfs/test_attribution.py | 243 ++++++++++++ tests/test_compare_vcfs/test_happy.py | 255 +++++++++++++ tests/test_compare_vcfs/test_integration.py | 159 ++++++++ tests/test_compare_vcfs/test_reports.py | 300 +++++++++++++++ tests/test_compare_vcfs/test_runner.py | 271 ++++++++++++++ .../test_simulation_summary.py | 348 ++++++++++++++++++ 18 files changed, 2622 insertions(+), 1 deletion(-) create mode 100644 neat/cli/commands/compare_vcfs.py create mode 100644 neat/compare_vcfs/__init__.py create mode 100644 neat/compare_vcfs/attribution.py create mode 100644 neat/compare_vcfs/happy.py create mode 100644 neat/compare_vcfs/reports.py create mode 100644 neat/compare_vcfs/runner.py create mode 100644 neat/read_simulator/utils/simulation_summary.py create mode 100644 tests/test_compare_vcfs/__init__.py create mode 100644 tests/test_compare_vcfs/test_attribution.py create mode 100644 tests/test_compare_vcfs/test_happy.py create mode 100644 tests/test_compare_vcfs/test_integration.py create mode 100644 tests/test_compare_vcfs/test_reports.py create mode 100644 tests/test_compare_vcfs/test_runner.py create mode 100644 tests/test_read_simulator/test_simulation_summary.py diff --git a/ChangeLog.md b/ChangeLog.md index 9e406c1a..12c0c8d2 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,38 @@ +# 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. + +**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. + # 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, diff --git a/README.md b/README.md index 5b568aff..c23354b3 100755 --- a/README.md +++ b/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,65 @@ 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] +``` + +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 | + +**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. diff --git a/neat/cli/commands/compare_vcfs.py b/neat/cli/commands/compare_vcfs.py new file mode 100644 index 00000000..994209a8 --- /dev/null +++ b/neat/cli/commands/compare_vcfs.py @@ -0,0 +1,78 @@ +""" +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 _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." + ) + + 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, + ) diff --git a/neat/compare_vcfs/__init__.py b/neat/compare_vcfs/__init__.py new file mode 100644 index 00000000..e0d6f5be --- /dev/null +++ b/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 * diff --git a/neat/compare_vcfs/attribution.py b/neat/compare_vcfs/attribution.py new file mode 100644 index 00000000..cb4d32f5 --- /dev/null +++ b/neat/compare_vcfs/attribution.py @@ -0,0 +1,131 @@ +""" +NEAT-aware false-negative attribution for `neat compare-vcfs`. + +Each FN from hap.py is tagged with one or more reasons drawn from the +simulator's own configuration: + + - `outside_simulated_contigs` — the FN's chromosome wasn't in the NEAT run. + - `outside_mutation_bed` — `mutation_bed` was set and the FN position + falls outside its regions. + - `outside_target_bed` — `target_bed` was set and the FN position + falls outside its regions. + - `unknown` — none of the above; NEAT can't explain it. + +If the FN's contig wasn't simulated at all, that single tag is reported on its +own (the BED checks would be meaningless). +""" +import logging +from pathlib import Path + +__all__ = [ + "attribute_fn", + "attribute_fns", + "load_bed_intervals", + "position_in_intervals", + "REASON_OUTSIDE_CONTIGS", + "REASON_OUTSIDE_MUTATION_BED", + "REASON_OUTSIDE_TARGET_BED", + "REASON_UNKNOWN", +] + +_LOG = logging.getLogger(__name__) + +REASON_OUTSIDE_CONTIGS = "outside_simulated_contigs" +REASON_OUTSIDE_MUTATION_BED = "outside_mutation_bed" +REASON_OUTSIDE_TARGET_BED = "outside_target_bed" +REASON_UNKNOWN = "unknown" + + +def load_bed_intervals(bed_path: Path | str | None) -> dict[str, list[tuple[int, int]]] | None: + """ + Parse a BED file into a per-contig list of sorted 0-based half-open + intervals. Returns None if `bed_path` is None. + + Comments (`#`), `track`, and `browser` header lines are skipped. Rows with + non-integer start/end are dropped with a debug log entry. + """ + if bed_path is None: + return None + intervals: dict[str, list[tuple[int, int]]] = {} + with open(bed_path) as fh: + for lineno, raw in enumerate(fh, start=1): + line = raw.strip() + if not line or line.startswith("#") or line.startswith(("track", "browser")): + continue + parts = line.split("\t") + if len(parts) < 3: + _LOG.debug(f"{bed_path}:{lineno}: skipping short line") + continue + try: + start, end = int(parts[1]), int(parts[2]) + except ValueError: + _LOG.debug(f"{bed_path}:{lineno}: skipping non-integer interval") + continue + intervals.setdefault(parts[0], []).append((start, end)) + for chrom in intervals: + intervals[chrom].sort() + return intervals + + +def position_in_intervals(pos: int, intervals: list[tuple[int, int]]) -> bool: + """ + Membership test for a 1-based VCF position against 0-based half-open BED + intervals sorted by start. Returns True if any interval contains the + position. Correct for overlapping intervals. + """ + pos_0 = pos - 1 + for start, end in intervals: + if start > pos_0: + return False # sorted: no further interval can match + if pos_0 < end: + return True + return False + + +def attribute_fn( + chrom: str, + pos: int, + contigs_simulated: set[str] | frozenset[str], + mutation_intervals: dict[str, list[tuple[int, int]]] | None, + target_intervals: dict[str, list[tuple[int, int]]] | None, +) -> list[str]: + """ + Return the list of NEAT-specific reasons that explain a single FN. + + If the contig wasn't simulated, that root cause is reported alone — BED + checks are skipped because they presuppose simulation. + """ + if chrom not in contigs_simulated: + return [REASON_OUTSIDE_CONTIGS] + + reasons: list[str] = [] + if mutation_intervals is not None: + chrom_intervals = mutation_intervals.get(chrom, []) + if not position_in_intervals(pos, chrom_intervals): + reasons.append(REASON_OUTSIDE_MUTATION_BED) + if target_intervals is not None: + chrom_intervals = target_intervals.get(chrom, []) + if not position_in_intervals(pos, chrom_intervals): + reasons.append(REASON_OUTSIDE_TARGET_BED) + if not reasons: + reasons.append(REASON_UNKNOWN) + return reasons + + +def attribute_fns(fn_records, summary: dict) -> list[tuple]: + """ + Tag every FN against the run's simulation_summary. + + :param fn_records: iterable of pysam.VariantRecord (FN bucket from hap.py). + :param summary: parsed simulation_summary.json. + :return: list of (record, reasons) tuples; `reasons` is a list[str]. + """ + contigs = frozenset(summary["delivered"].get("contigs_simulated", [])) + cfg = summary.get("config", {}) + mutation_intervals = load_bed_intervals(cfg.get("mutation_bed")) + target_intervals = load_bed_intervals(cfg.get("target_bed")) + + return [ + (rec, attribute_fn(rec.chrom, rec.pos, contigs, mutation_intervals, target_intervals)) + for rec in fn_records + ] diff --git a/neat/compare_vcfs/happy.py b/neat/compare_vcfs/happy.py new file mode 100644 index 00000000..987cfb75 --- /dev/null +++ b/neat/compare_vcfs/happy.py @@ -0,0 +1,138 @@ +""" +hap.py subprocess invocation + output-VCF parsing for `neat compare-vcfs`. + +hap.py emits one VCF per run with per-sample FORMAT annotations describing each +record's classification: + - sample 0 (TRUTH): BD == 'TP' (matched), 'FN' (missed), or '.' (no call) + - sample 1 (QUERY): BD == 'TP' (matched), 'FP' (spurious), or '.' (no call) + +This module wraps the subprocess and classifies each record into one of TP/FN/FP +based on those FORMAT fields. Multi-allelic and gVCF handling are delegated to +hap.py upstream; we read what it emits. +""" +import logging +import subprocess +from pathlib import Path +from typing import Iterable + +import pysam + +__all__ = [ + "HappyExecutionError", + "HappyParseError", + "run_happy", + "parse_happy_output", +] + +_LOG = logging.getLogger(__name__) + + +class HappyExecutionError(RuntimeError): + """Raised when the hap.py subprocess fails.""" + + +class HappyParseError(RuntimeError): + """Raised when the hap.py output VCF is missing expected structure.""" + + +def run_happy( + happy_bin: Path, + golden_vcf: Path, + called_vcf: Path, + output_prefix: Path, + reference: Path | None = None, + target_bed: Path | None = None, + extra_args: Iterable[str] = (), +) -> Path: + """ + Invoke hap.py and return the path to its bgzipped output VCF. + + :param happy_bin: Absolute path to the hap.py executable. + :param golden_vcf: Truth VCF (NEAT golden). + :param called_vcf: Query VCF (downstream variant caller). + :param output_prefix: Path prefix; hap.py appends `.vcf.gz` and writes + sibling artifacts (summary.csv, extended.csv, runinfo.json). + :param reference: Optional reference FASTA forwarded as `-r`. + :param target_bed: Optional restricting BED forwarded as `-T`. + :param extra_args: Optional positional pass-through for future tunables. + :return: Path to `.vcf.gz`. + :raises HappyExecutionError: if hap.py exits non-zero or its output VCF is absent. + """ + cmd = [ + str(happy_bin), str(golden_vcf), str(called_vcf), + "-o", str(output_prefix), + ] + if reference is not None: + cmd += ["-r", str(reference)] + if target_bed is not None: + cmd += ["-T", str(target_bed)] + cmd += list(extra_args) + + _LOG.info(f"Running hap.py: {' '.join(cmd)}") + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + _LOG.error(f"hap.py stderr:\n{result.stderr}") + raise HappyExecutionError( + f"hap.py exited {result.returncode}. See log for stderr. " + f"Command: {' '.join(cmd)}" + ) + if result.stderr: + _LOG.debug(f"hap.py stderr:\n{result.stderr}") + + output_vcf = Path(str(output_prefix) + ".vcf.gz") + if not output_vcf.is_file(): + raise HappyExecutionError( + f"hap.py reported success but its output VCF is missing: {output_vcf}" + ) + return output_vcf + + +def parse_happy_output(vcf_path: Path) -> dict[str, list]: + """ + Group hap.py output records into TP/FN/FP buckets. + + :param vcf_path: hap.py's bgzipped output VCF. + :return: dict with keys 'TP', 'FN', 'FP'; values are lists of pysam.VariantRecord. + :raises HappyParseError: if the VCF lacks the expected TRUTH/QUERY samples + or the BD FORMAT field. + """ + buckets: dict[str, list] = {"TP": [], "FN": [], "FP": []} + with pysam.VariantFile(str(vcf_path)) as vf: + sample_names = list(vf.header.samples) + if len(sample_names) < 2: + raise HappyParseError( + f"{vcf_path} has {len(sample_names)} sample(s); hap.py output must have " + f"TRUTH and QUERY samples." + ) + if "BD" not in vf.header.formats: + raise HappyParseError( + f"{vcf_path} is missing the BD FORMAT field — not a hap.py output VCF?" + ) + truth_name, query_name = sample_names[0], sample_names[1] + + for rec in vf: + truth_bd = rec.samples[truth_name].get("BD") or "." + query_bd = rec.samples[query_name].get("BD") or "." + classification = _classify(truth_bd, query_bd) + if classification is not None: + buckets[classification].append(rec) + return buckets + + +def _classify(truth_bd: str, query_bd: str) -> str | None: + """ + Apply the hap.py decision rules. + + A record is: + - TP if either sample reports TP (matched on at least one side) + - FN if truth reports FN and query did not match + - FP if query reports FP and truth did not match + - Otherwise None (no-call / nocomp / hap.py-internal types) + """ + if truth_bd == "TP" or query_bd == "TP": + return "TP" + if truth_bd == "FN": + return "FN" + if query_bd == "FP": + return "FP" + return None diff --git a/neat/compare_vcfs/reports.py b/neat/compare_vcfs/reports.py new file mode 100644 index 00000000..80ff4491 --- /dev/null +++ b/neat/compare_vcfs/reports.py @@ -0,0 +1,244 @@ +""" +Report generation for `neat compare-vcfs`. + +Three artifacts: + - comparison_summary.json — machine-readable rollup + - comparison_summary.txt — human-readable rollup + - FN_with_reasons.vcf — hap.py's FN records, annotated with NEAT_REASON +""" +import json +import logging +import os +import time +from datetime import datetime, timezone +from pathlib import Path + +import pysam + +from .. import __version__ as NEAT_VERSION + +__all__ = [ + "REPORT_SCHEMA_VERSION", + "build_comparison_summary", + "compute_metrics", + "render_summary_txt", + "summarize_fn_reasons", + "write_comparison_summary_json", + "write_comparison_summary_txt", + "write_fn_attribution_plot", + "write_fn_with_reasons", +] + +_LOG = logging.getLogger(__name__) + +REPORT_SCHEMA_VERSION = "1" + + +def compute_metrics(counts: dict) -> dict: + """ + Precision, recall, F1 from TP/FN/FP. Returns None for any metric that is + undefined (e.g., empty truth set yields undefined recall). + """ + tp = counts.get("TP", 0) + fn = counts.get("FN", 0) + fp = counts.get("FP", 0) + precision = tp / (tp + fp) if (tp + fp) > 0 else None + recall = tp / (tp + fn) if (tp + fn) > 0 else None + if precision is None or recall is None or (precision + recall) == 0: + f1 = None + else: + f1 = 2 * precision * recall / (precision + recall) + return {"precision": precision, "recall": recall, "f1": f1} + + +def summarize_fn_reasons(fn_reasons) -> dict[str, int]: + """Roll up reason tags across all FNs to a {reason: count} dict.""" + counts: dict[str, int] = {} + for _, reasons in fn_reasons: + for r in reasons: + counts[r] = counts.get(r, 0) + 1 + return counts + + +def build_comparison_summary( + *, + golden_vcf: Path, + called_vcf: Path, + neat_run_dir: Path, + simulation_summary_path: Path, + happy_output_vcf: Path, + happy_output_prefix: Path, + counts: dict, + fn_attribution: dict, + fn_with_reasons_vcf: Path, + comparison_summary_json: Path, + comparison_summary_txt: Path, +) -> dict: + """Assemble the comparison_summary dict from the run's artifacts.""" + return { + "schema_version": REPORT_SCHEMA_VERSION, + "neat_version": NEAT_VERSION, + "generated_at": _iso_utc(time.time()), + "inputs": { + "golden_vcf": str(Path(golden_vcf).resolve()), + "called_vcf": str(Path(called_vcf).resolve()), + "neat_run_dir": str(Path(neat_run_dir).resolve()), + "simulation_summary": str(Path(simulation_summary_path).resolve()), + }, + "happy": { + "output_prefix": str(Path(happy_output_prefix).resolve()), + "output_vcf": str(Path(happy_output_vcf).resolve()), + }, + "counts": dict(counts), + "metrics": compute_metrics(counts), + "fn_attribution": dict(fn_attribution), + "outputs": { + "fn_with_reasons_vcf": str(Path(fn_with_reasons_vcf).resolve()), + "comparison_summary_json": str(Path(comparison_summary_json).resolve()), + "comparison_summary_txt": str(Path(comparison_summary_txt).resolve()), + }, + } + + +def write_comparison_summary_json(summary: dict, path: Path) -> Path: + """Atomically write the summary as pretty JSON.""" + path = Path(path) + tmp = path.with_suffix(".json.tmp") + with open(tmp, "w") as fh: + json.dump(summary, fh, indent=2) + os.replace(tmp, path) + return path + + +def write_comparison_summary_txt(summary: dict, path: Path) -> Path: + """Atomically write a human-readable rollup.""" + path = Path(path) + tmp = path.with_suffix(".txt.tmp") + tmp.write_text(render_summary_txt(summary)) + os.replace(tmp, path) + return path + + +def render_summary_txt(summary: dict) -> str: + """Render the summary dict as a fixed-width text report.""" + counts = summary["counts"] + metrics = summary["metrics"] + fn_attr = summary["fn_attribution"] + inputs = summary["inputs"] + outputs = summary["outputs"] + + lines: list[str] = [ + "NEAT compare-vcfs report", + "========================", + "", + f"Generated: {summary['generated_at']}", + f"NEAT version: {summary['neat_version']}", + "", + "Inputs", + "------", + f" Truth (golden) VCF: {inputs['golden_vcf']}", + f" Called VCF: {inputs['called_vcf']}", + f" NEAT run dir: {inputs['neat_run_dir']}", + "", + "Classification (from hap.py)", + "----------------------------", + f" True positives (TP): {counts.get('TP', 0)}", + f" False negatives (FN): {counts.get('FN', 0)}", + f" False positives (FP): {counts.get('FP', 0)}", + "", + "Metrics", + "-------", + f" Precision: {_fmt_metric(metrics['precision'])} (TP / (TP + FP))", + f" Recall: {_fmt_metric(metrics['recall'])} (TP / (TP + FN))", + f" F1: {_fmt_metric(metrics['f1'])}", + "", + "FN attribution", + "--------------", + ] + if fn_attr: + width = max(len(k) for k in fn_attr) + 2 + for reason in sorted(fn_attr): + lines.append(f" {reason:<{width}} {fn_attr[reason]}") + else: + lines.append(" (no false negatives)") + + lines += [ + "", + "Outputs", + "-------", + f" Annotated FN VCF: {outputs['fn_with_reasons_vcf']}", + f" This report (JSON): {outputs['comparison_summary_json']}", + f" This report (text): {outputs['comparison_summary_txt']}", + "", + ] + return "\n".join(lines) + + +def write_fn_with_reasons( + source_happy_vcf: Path, + fn_reasons: list, + output_path: Path, +) -> Path: + """ + Write only FN records to `output_path`, each annotated with a NEAT_REASON + INFO tag. The source VCF is opened only to obtain a compatible header; + records come from the already-classified `fn_reasons` list. + """ + source_happy_vcf = Path(source_happy_vcf) + output_path = Path(output_path) + with pysam.VariantFile(str(source_happy_vcf)) as src: + new_header = src.header.copy() + new_header.info.add( + "NEAT_REASON", ".", "String", + "Comma-separated NEAT-aware false-negative attribution reasons", + ) + with pysam.VariantFile(str(output_path), "w", header=new_header) as dst: + for rec, reasons in fn_reasons: + rec.translate(new_header) + rec.info["NEAT_REASON"] = ",".join(reasons) + dst.write(rec) + return output_path + + +def write_fn_attribution_plot(fn_attribution: dict[str, int], path: Path) -> Path: + """ + Render a horizontal bar chart of FN-reason counts to `path` (PNG). + + A true Venn diagram is awkward with this scheme (outside_simulated_contigs + is mutually exclusive with the BED reasons; unknown can't co-occur with + anything), so the artifact is a clearer bar chart. + """ + # Import here so users running without --plot don't pay the matplotlib import cost + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + path = Path(path) + reasons = sorted(fn_attribution) + counts = [fn_attribution[r] for r in reasons] + + fig, ax = plt.subplots(figsize=(8, max(2.0, 0.5 * len(reasons) + 1.0))) + if reasons: + ax.barh(reasons, counts, color="steelblue") + for i, count in enumerate(counts): + ax.text(count, i, f" {count}", va="center") + ax.set_xlim(0, max(counts) * 1.15 if max(counts) else 1) + else: + ax.text(0.5, 0.5, "no false negatives", ha="center", va="center", + transform=ax.transAxes) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_xlabel("FN count") + ax.set_title("False-negative attribution") + fig.tight_layout() + fig.savefig(path, dpi=120) + plt.close(fig) + return path + + +def _iso_utc(epoch_seconds: float) -> str: + return datetime.fromtimestamp(epoch_seconds, tz=timezone.utc).isoformat(timespec="seconds").replace("+00:00", "Z") + + +def _fmt_metric(v) -> str: + return "N/A" if v is None else f"{v:.4f}" diff --git a/neat/compare_vcfs/runner.py b/neat/compare_vcfs/runner.py new file mode 100644 index 00000000..fdc1777c --- /dev/null +++ b/neat/compare_vcfs/runner.py @@ -0,0 +1,203 @@ +""" +Runner for the `neat compare-vcfs` subcommand (issue #297). + +Compares a downstream variant caller's VCF against a NEAT-simulated truth VCF +and attributes false negatives to the simulator's own configuration (mutation +bed, target bed, simulated contigs) by reading `simulation_summary.json`. + +This module currently implements **only the input-validation scaffold**: +- Validates input file paths. +- Locates `hap.py` on PATH (or via `--happy-bin`). +- Loads and version-checks `simulation_summary.json`. +- Creates the output directory. + +The actual hap.py invocation, output parsing, FN attribution, and report +generation are tracked as follow-up work on the same issue. +""" +import json +import logging +import shutil +from pathlib import Path + +from ..common import validate_input_path +from ..read_simulator.utils.simulation_summary import SCHEMA_VERSION +from .attribution import attribute_fns +from .happy import run_happy, parse_happy_output +from .reports import ( + build_comparison_summary, + summarize_fn_reasons, + write_comparison_summary_json, + write_comparison_summary_txt, + write_fn_attribution_plot, + write_fn_with_reasons, +) + +__all__ = ["compare_vcfs_runner", "load_simulation_summary", "discover_happy"] + +_LOG = logging.getLogger(__name__) + +# An install hint surfaced when hap.py is not on PATH and --happy-bin is not given. +_HAPPY_INSTALL_HINT = ( + "hap.py was not found on $PATH. Install via " + "`conda install -c bioconda hap.py`, " + "or use the official Docker image, or pass --happy-bin /path/to/hap.py." +) + + +class HappyNotFoundError(RuntimeError): + """Raised when hap.py cannot be located.""" + + +class SimulationSummaryError(RuntimeError): + """Raised when simulation_summary.json is missing, malformed, or incompatible.""" + + +def compare_vcfs_runner( + golden_vcf: str, + called_vcf: str, + neat_run_dir: str, + output_dir: str, + reference: str | None = None, + target_bed: str | None = None, + happy_bin: str | None = None, + plot: bool = False, +): + """ + Run the comparison pipeline. + + :param golden_vcf: NEAT-simulated truth VCF. + :param called_vcf: Downstream variant caller's VCF. + :param neat_run_dir: Directory containing the NEAT simulator output, including + `simulation_summary.json`. + :param output_dir: Where to write `comparison_summary.{json,txt}` and + `FN_with_reasons.vcf`. Created if it does not exist. + :param reference: Optional path to the reference FASTA (forwarded to hap.py). + :param target_bed: Optional BED of regions to restrict comparison to. + :param happy_bin: Optional explicit path to the hap.py binary. + """ + golden_path = Path(golden_vcf).resolve() + called_path = Path(called_vcf).resolve() + run_dir_path = Path(neat_run_dir).resolve() + out_dir_path = Path(output_dir).resolve() + + validate_input_path(golden_path) + validate_input_path(called_path) + if not run_dir_path.is_dir(): + raise FileNotFoundError(f"--neat-run-dir does not exist or is not a directory: {run_dir_path}") + + if reference is not None: + validate_input_path(Path(reference).resolve()) + if target_bed is not None: + validate_input_path(Path(target_bed).resolve()) + + happy = discover_happy(happy_bin) + _LOG.info(f"Using hap.py at: {happy}") + + summary = load_simulation_summary(run_dir_path) + _LOG.info( + f"Loaded simulation_summary.json (schema v{summary['schema_version']}, " + f"neat {summary['neat_version']}, " + f"{summary['delivered'].get('total_variants')} simulated variants " + f"across {len(summary['delivered'].get('contigs_simulated', []))} contigs)." + ) + + out_dir_path.mkdir(parents=True, exist_ok=True) + _LOG.info(f"Output directory: {out_dir_path}") + + happy_prefix = out_dir_path / "happy" + happy_vcf = run_happy( + happy_bin=happy, + golden_vcf=golden_path, + called_vcf=called_path, + output_prefix=happy_prefix, + reference=Path(reference).resolve() if reference else None, + target_bed=Path(target_bed).resolve() if target_bed else None, + ) + buckets = parse_happy_output(happy_vcf) + _LOG.info( + f"hap.py classification: TP={len(buckets['TP'])} " + f"FN={len(buckets['FN'])} FP={len(buckets['FP'])}" + ) + + fn_reasons = attribute_fns(buckets["FN"], summary) + reason_counts = summarize_fn_reasons(fn_reasons) + if fn_reasons: + _LOG.info(f"FN attribution: {reason_counts}") + + fn_with_reasons_path = out_dir_path / "FN_with_reasons.vcf" + comparison_json_path = out_dir_path / "comparison_summary.json" + comparison_txt_path = out_dir_path / "comparison_summary.txt" + + write_fn_with_reasons(happy_vcf, fn_reasons, fn_with_reasons_path) + _LOG.info(f"Wrote {fn_with_reasons_path}") + + counts = {k: len(v) for k, v in buckets.items()} + report = build_comparison_summary( + golden_vcf=golden_path, + called_vcf=called_path, + neat_run_dir=run_dir_path, + simulation_summary_path=run_dir_path / "simulation_summary.json", + happy_output_vcf=happy_vcf, + happy_output_prefix=happy_prefix, + counts=counts, + fn_attribution=reason_counts, + fn_with_reasons_vcf=fn_with_reasons_path, + comparison_summary_json=comparison_json_path, + comparison_summary_txt=comparison_txt_path, + ) + write_comparison_summary_json(report, comparison_json_path) + write_comparison_summary_txt(report, comparison_txt_path) + _LOG.info(f"Wrote {comparison_json_path}") + _LOG.info(f"Wrote {comparison_txt_path}") + + if plot: + plot_path = out_dir_path / "fn_attribution.png" + write_fn_attribution_plot(reason_counts, plot_path) + _LOG.info(f"Wrote {plot_path}") + + +def discover_happy(explicit_path: str | None) -> Path: + """ + Resolve the hap.py binary. Explicit path wins; otherwise look on $PATH. + + :param explicit_path: Path passed via --happy-bin, or None. + :return: Absolute path to a hap.py executable. + :raises HappyNotFoundError: if neither location yields an executable. + """ + if explicit_path is not None: + p = Path(explicit_path).resolve() + if not p.is_file(): + raise HappyNotFoundError(f"--happy-bin path does not exist: {p}") + return p + found = shutil.which("hap.py") + if found is None: + raise HappyNotFoundError(_HAPPY_INSTALL_HINT) + return Path(found).resolve() + + +def load_simulation_summary(neat_run_dir: Path) -> dict: + """ + Read and validate `simulation_summary.json` from the NEAT run directory. + + :raises SimulationSummaryError: if the file is missing, malformed, or its + schema version is incompatible with this build. + """ + summary_path = neat_run_dir / "simulation_summary.json" + if not summary_path.is_file(): + raise SimulationSummaryError( + f"{summary_path} not found. Did the simulator run finish? " + "Re-run `neat read-simulator` to produce one." + ) + try: + with open(summary_path) as fh: + data = json.load(fh) + except json.JSONDecodeError as exc: + raise SimulationSummaryError(f"{summary_path} is not valid JSON: {exc}") from exc + + actual = data.get("schema_version") + if actual != SCHEMA_VERSION: + raise SimulationSummaryError( + f"{summary_path} schema_version is {actual!r}, expected {SCHEMA_VERSION!r}. " + "Regenerate by re-running `neat read-simulator` with a current NEAT build." + ) + return data diff --git a/neat/read_simulator/runner.py b/neat/read_simulator/runner.py index 22574ff9..6b3d541d 100644 --- a/neat/read_simulator/runner.py +++ b/neat/read_simulator/runner.py @@ -25,6 +25,7 @@ from ..variants import ContigVariants from .utils.split_inputs import main as split_main from .utils.stitch_outputs import main as stitch_main +from .utils.simulation_summary import write_simulation_summary __all__ = ["read_simulator_runner"] @@ -326,6 +327,15 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): _LOG.info(f"Read simulator complete in {time.time() - analysis_start} s") + write_simulation_summary( + options=options, + output_dir=output_dir, + file_prefix=str(file_prefix), + config_path=config, + analysis_start=analysis_start, + contigs_simulated=list(input_variants_dict.keys()), + ) + def filter_thread_variants(contig_variants: ContigVariants, coords: tuple[int, int]) -> ContigVariants: ret_contig_vars = ContigVariants() for variant_loc in contig_variants.variant_locations: diff --git a/neat/read_simulator/utils/simulation_summary.py b/neat/read_simulator/utils/simulation_summary.py new file mode 100644 index 00000000..e8db4d06 --- /dev/null +++ b/neat/read_simulator/utils/simulation_summary.py @@ -0,0 +1,144 @@ +""" +Emit simulation_summary.json next to the other simulator outputs. + +Consumed by `neat compare-vcfs` (issue #297) to attribute false negatives from a +downstream variant caller against the simulator's configuration (mutation bed, +target bed, simulated contigs). +""" +import json +import logging +import os +import time +from datetime import datetime, timezone +from pathlib import Path +from typing import Iterable + +import pysam + +from ... import __version__ as NEAT_VERSION + +_LOG = logging.getLogger(__name__) + +SCHEMA_VERSION = "1" + + +def write_simulation_summary( + options, + output_dir: Path, + file_prefix: str, + config_path: Path, + analysis_start: float, + contigs_simulated: Iterable[str], +) -> Path: + """ + Write `simulation_summary.json` into `output_dir`. + + Returns the path to the written file. + """ + started_at = _iso_utc(analysis_start) + completed_ts = time.time() + completed_at = _iso_utc(completed_ts) + + vcf_path = _abs_or_none(getattr(options, "vcf", None)) + bam_path = _abs_or_none(getattr(options, "bam", None)) + fq1_path = _abs_or_none(getattr(options, "fq1", None)) + fq2_path = _abs_or_none(getattr(options, "fq2", None)) + + total_variants, variants_by_contig = _count_variants(vcf_path) + total_reads = _count_reads(bam_path, fq1_path, fq2_path, bool(getattr(options, "paired_ended", False))) + + fastq_outputs = [p for p in (fq1_path, fq2_path) if p is not None] + + summary = { + "schema_version": SCHEMA_VERSION, + "neat_version": NEAT_VERSION, + "run": { + "started_at": started_at, + "completed_at": completed_at, + "duration_seconds": round(completed_ts - analysis_start, 3), + "config_file": _abs_or_none(config_path), + "output_dir": str(Path(output_dir).resolve()), + "output_prefix": file_prefix, + }, + "config": { + "reference": _abs_or_none(getattr(options, "reference", None)), + "coverage": getattr(options, "coverage", None), + "read_len": getattr(options, "read_len", None), + "paired_ended": getattr(options, "paired_ended", None), + "fragment_mean": getattr(options, "fragment_mean", None), + "fragment_st_dev": getattr(options, "fragment_st_dev", None), + "ploidy": getattr(options, "ploidy", None), + "rng_seed": getattr(options, "rng_seed", None), + "threads": getattr(options, "threads", None), + "mutation_rate": getattr(options, "mutation_rate", None), + "mutation_bed": _abs_or_none(getattr(options, "mutation_bed", None)), + "target_bed": _abs_or_none(getattr(options, "target_bed", None)), + "discard_bed": _abs_or_none(getattr(options, "discard_bed", None)), + "include_vcf": _abs_or_none(getattr(options, "include_vcf", None)), + "mutation_model": _abs_or_none(getattr(options, "mutation_model", None)), + "gc_model": _abs_or_none(getattr(options, "gc_model", None)), + "error_model": _abs_or_none(getattr(options, "error_model", None)), + "fragment_model": _abs_or_none(getattr(options, "fragment_model", None)), + }, + "outputs": { + "fastq": fastq_outputs if fastq_outputs else None, + "bam": bam_path, + "vcf": vcf_path, + }, + "delivered": { + "total_reads": total_reads, + "total_variants": total_variants, + "variants_by_contig": variants_by_contig, + "contigs_simulated": list(contigs_simulated), + }, + } + + out_path = Path(output_dir) / "simulation_summary.json" + tmp_path = out_path.with_suffix(".json.tmp") + with open(tmp_path, "w") as fh: + json.dump(summary, fh, indent=2) + os.replace(tmp_path, out_path) + _LOG.info(f"Wrote {out_path}") + return out_path + + +def _iso_utc(epoch_seconds: float) -> str: + return datetime.fromtimestamp(epoch_seconds, tz=timezone.utc).isoformat(timespec="seconds").replace("+00:00", "Z") + + +def _abs_or_none(value): + if value is None: + return None + return str(Path(value).resolve()) + + +def _count_variants(vcf_path): + if vcf_path is None or not Path(vcf_path).is_file(): + return None, None + total = 0 + by_contig: dict[str, int] = {} + try: + with pysam.VariantFile(vcf_path) as vf: + for rec in vf: + total += 1 + by_contig[rec.chrom] = by_contig.get(rec.chrom, 0) + 1 + except Exception as exc: + _LOG.warning(f"Could not count variants in {vcf_path}: {exc}") + return None, None + return total, by_contig + + +def _count_reads(bam_path, fq1_path, fq2_path, paired_ended): + if bam_path is not None and Path(bam_path).is_file(): + try: + with pysam.AlignmentFile(bam_path, "rb") as bf: + return bf.count(until_eof=True) + except Exception as exc: + _LOG.warning(f"Could not count reads in {bam_path}: {exc}") + if fq1_path is not None and Path(fq1_path).is_file(): + try: + n = sum(1 for _ in pysam.FastxFile(fq1_path)) + return n * 2 if paired_ended else n + except Exception as exc: + _LOG.warning(f"Could not count reads in {fq1_path}: {exc}") + return None diff --git a/pyproject.toml b/pyproject.toml index fad2097a..6b5c73d7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "neat-genreads" -version = "4.4.4" +version = "4.5.0" description = "NGS Simulation toolkit" readme = "README.md" authors = ["Joshua Allen ", "Keshav Gandhi "] diff --git a/tests/test_compare_vcfs/__init__.py b/tests/test_compare_vcfs/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_compare_vcfs/test_attribution.py b/tests/test_compare_vcfs/test_attribution.py new file mode 100644 index 00000000..47b2a9c3 --- /dev/null +++ b/tests/test_compare_vcfs/test_attribution.py @@ -0,0 +1,243 @@ +""" +Tests for neat/compare_vcfs/attribution.py — NEAT-aware FN classification. +""" +from pathlib import Path +from types import SimpleNamespace + +import pytest + +from neat.compare_vcfs.attribution import ( + REASON_OUTSIDE_CONTIGS, + REASON_OUTSIDE_MUTATION_BED, + REASON_OUTSIDE_TARGET_BED, + REASON_UNKNOWN, + attribute_fn, + attribute_fns, + load_bed_intervals, + position_in_intervals, +) + + +# --------------------------------------------------------------------------- +# load_bed_intervals +# --------------------------------------------------------------------------- + +def test_load_bed_intervals_returns_none_for_none_path(): + assert load_bed_intervals(None) is None + + +def test_load_bed_intervals_parses_basic_bed(tmp_path): + bed = tmp_path / "x.bed" + bed.write_text("chr1\t100\t200\nchr1\t300\t400\nchr2\t50\t75\n") + result = load_bed_intervals(bed) + assert result == {"chr1": [(100, 200), (300, 400)], "chr2": [(50, 75)]} + + +def test_load_bed_intervals_skips_comments_and_headers(tmp_path): + bed = tmp_path / "x.bed" + bed.write_text( + "# header comment\n" + "track name=foo\n" + "browser position chr1\n" + "\n" + "chr1\t100\t200\n" + ) + assert load_bed_intervals(bed) == {"chr1": [(100, 200)]} + + +def test_load_bed_intervals_skips_malformed_rows(tmp_path): + bed = tmp_path / "x.bed" + bed.write_text( + "chr1\t100\t200\n" + "short_row\n" + "chr1\tNaN\t300\n" + "chr2\t10\t20\n" + ) + assert load_bed_intervals(bed) == {"chr1": [(100, 200)], "chr2": [(10, 20)]} + + +def test_load_bed_intervals_sorts_intervals_per_contig(tmp_path): + bed = tmp_path / "x.bed" + bed.write_text("chr1\t300\t400\nchr1\t100\t200\nchr1\t500\t600\n") + result = load_bed_intervals(bed) + assert result["chr1"] == [(100, 200), (300, 400), (500, 600)] + + +def test_load_bed_intervals_ignores_extra_columns(tmp_path): + bed = tmp_path / "x.bed" + bed.write_text("chr1\t100\t200\tname1\t1000\t+\n") + assert load_bed_intervals(bed) == {"chr1": [(100, 200)]} + + +def test_load_bed_intervals_returns_empty_dict_for_empty_file(tmp_path): + bed = tmp_path / "empty.bed" + bed.write_text("") + assert load_bed_intervals(bed) == {} + + +# --------------------------------------------------------------------------- +# position_in_intervals — coordinate boundaries are the tricky part +# --------------------------------------------------------------------------- + +def test_position_in_intervals_inside_interval(): + """VCF pos 150 → 0-based 149; falls in [100, 200).""" + assert position_in_intervals(150, [(100, 200)]) is True + + +def test_position_in_intervals_at_start_inclusive(): + """BED start is inclusive. VCF pos 101 → 0-based 100; [100, 200) starts at 100.""" + assert position_in_intervals(101, [(100, 200)]) is True + + +def test_position_in_intervals_at_end_inclusive_for_1based(): + """BED end is exclusive. VCF pos 200 → 0-based 199; in [100, 200).""" + assert position_in_intervals(200, [(100, 200)]) is True + + +def test_position_in_intervals_just_past_end(): + """VCF pos 201 → 0-based 200; [100, 200) does not include 200.""" + assert position_in_intervals(201, [(100, 200)]) is False + + +def test_position_in_intervals_before_first(): + assert position_in_intervals(50, [(100, 200)]) is False + + +def test_position_in_intervals_between_intervals(): + assert position_in_intervals(250, [(100, 200), (300, 400)]) is False + + +def test_position_in_intervals_picks_right_interval_in_sorted_list(): + intervals = [(100, 200), (300, 400), (500, 600)] + assert position_in_intervals(350, intervals) is True + assert position_in_intervals(550, intervals) is True + assert position_in_intervals(450, intervals) is False + + +def test_position_in_intervals_handles_empty_list(): + assert position_in_intervals(100, []) is False + + +def test_position_in_intervals_finds_match_in_overlapping_intervals(): + """Overlapping intervals: ensure we don't miss containment in an earlier interval.""" + intervals = [(100, 500), (200, 300)] # overlapping; sorted by start + assert position_in_intervals(450, intervals) is True + + +# --------------------------------------------------------------------------- +# attribute_fn — single-record classification +# --------------------------------------------------------------------------- + +def test_attribute_fn_returns_outside_contigs_when_chrom_not_simulated(): + reasons = attribute_fn("chrZ", 100, frozenset({"chr1"}), None, None) + assert reasons == [REASON_OUTSIDE_CONTIGS] + + +def test_attribute_fn_outside_contigs_does_not_combine_with_bed_reasons(): + """If the contig isn't simulated, BED checks are skipped entirely.""" + reasons = attribute_fn("chrZ", 100, frozenset({"chr1"}), + {"chr1": [(0, 1000)]}, {"chr1": [(0, 1000)]}) + assert reasons == [REASON_OUTSIDE_CONTIGS] + + +def test_attribute_fn_unknown_when_no_beds_configured(): + reasons = attribute_fn("chr1", 100, frozenset({"chr1"}), None, None) + assert reasons == [REASON_UNKNOWN] + + +def test_attribute_fn_outside_mutation_bed_only(): + reasons = attribute_fn( + "chr1", 100, frozenset({"chr1"}), + mutation_intervals={"chr1": [(500, 600)]}, # 100 outside + target_intervals=None, + ) + assert reasons == [REASON_OUTSIDE_MUTATION_BED] + + +def test_attribute_fn_outside_target_bed_only(): + reasons = attribute_fn( + "chr1", 100, frozenset({"chr1"}), + mutation_intervals=None, + target_intervals={"chr1": [(500, 600)]}, + ) + assert reasons == [REASON_OUTSIDE_TARGET_BED] + + +def test_attribute_fn_multiple_reasons_combined(): + """A FN outside both beds gets both tags, in canonical order.""" + reasons = attribute_fn( + "chr1", 100, frozenset({"chr1"}), + mutation_intervals={"chr1": [(500, 600)]}, + target_intervals={"chr1": [(700, 800)]}, + ) + assert reasons == [REASON_OUTSIDE_MUTATION_BED, REASON_OUTSIDE_TARGET_BED] + + +def test_attribute_fn_unknown_when_inside_all_configured_beds(): + """If the FN is inside every configured bed, NEAT has no explanation.""" + reasons = attribute_fn( + "chr1", 150, frozenset({"chr1"}), + mutation_intervals={"chr1": [(100, 200)]}, + target_intervals={"chr1": [(100, 200)]}, + ) + assert reasons == [REASON_UNKNOWN] + + +def test_attribute_fn_bed_missing_chrom_is_outside(): + """A bed that doesn't mention this chrom counts as 'outside' for it.""" + reasons = attribute_fn( + "chr1", 100, frozenset({"chr1"}), + mutation_intervals={"chr2": [(0, 1000)]}, # chr1 absent + target_intervals=None, + ) + assert reasons == [REASON_OUTSIDE_MUTATION_BED] + + +# --------------------------------------------------------------------------- +# attribute_fns — the integration entry point used by the runner +# --------------------------------------------------------------------------- + +def _fake_record(chrom: str, pos: int): + return SimpleNamespace(chrom=chrom, pos=pos) + + +def test_attribute_fns_returns_one_tag_set_per_record(tmp_path): + mut = tmp_path / "mut.bed" + mut.write_text("chr1\t0\t1000\n") + summary = { + "delivered": {"contigs_simulated": ["chr1", "chr2"]}, + "config": {"mutation_bed": str(mut), "target_bed": None}, + } + fns = [_fake_record("chr1", 500), _fake_record("chr2", 100), _fake_record("chrZ", 100)] + result = attribute_fns(fns, summary) + assert [r for _, r in result] == [ + [REASON_UNKNOWN], # chr1:500 inside mut bed + [REASON_OUTSIDE_MUTATION_BED], # chr2:100 — chr2 not in mut bed + [REASON_OUTSIDE_CONTIGS], # chrZ wasn't simulated + ] + + +def test_attribute_fns_no_beds_configured(tmp_path): + summary = { + "delivered": {"contigs_simulated": ["chr1"]}, + "config": {"mutation_bed": None, "target_bed": None}, + } + fns = [_fake_record("chr1", 100)] + result = attribute_fns(fns, summary) + assert result[0][1] == [REASON_UNKNOWN] + + +def test_attribute_fns_returns_empty_list_for_empty_input(): + summary = {"delivered": {"contigs_simulated": []}, "config": {}} + assert attribute_fns([], summary) == [] + + +def test_attribute_fns_pairs_each_record_with_reasons(tmp_path): + summary = { + "delivered": {"contigs_simulated": ["chr1"]}, + "config": {"mutation_bed": None, "target_bed": None}, + } + rec = _fake_record("chr1", 42) + [(returned_rec, reasons)] = attribute_fns([rec], summary) + assert returned_rec is rec + assert reasons == [REASON_UNKNOWN] diff --git a/tests/test_compare_vcfs/test_happy.py b/tests/test_compare_vcfs/test_happy.py new file mode 100644 index 00000000..558c15d4 --- /dev/null +++ b/tests/test_compare_vcfs/test_happy.py @@ -0,0 +1,255 @@ +""" +Tests for neat/compare_vcfs/happy.py — hap.py subprocess invocation and +output-VCF parsing into TP/FN/FP buckets. +""" +import subprocess +from pathlib import Path +from types import SimpleNamespace + +import pysam +import pytest + +from neat.compare_vcfs.happy import ( + HappyExecutionError, + HappyParseError, + _classify, + parse_happy_output, + run_happy, +) + + +# --------------------------------------------------------------------------- +# run_happy +# --------------------------------------------------------------------------- + +def test_run_happy_builds_expected_command(tmp_path, monkeypatch): + """The subprocess invocation must include all forwarded options in order.""" + captured = {} + + def fake_run(cmd, capture_output, text): + captured["cmd"] = cmd + # Pretend hap.py wrote its output + Path(str(tmp_path / "happy") + ".vcf.gz").touch() + return SimpleNamespace(returncode=0, stdout="", stderr="") + + monkeypatch.setattr(subprocess, "run", fake_run) + + run_happy( + happy_bin=Path("/bin/hap.py"), + golden_vcf=Path("/tmp/g.vcf"), + called_vcf=Path("/tmp/c.vcf"), + output_prefix=tmp_path / "happy", + reference=Path("/tmp/ref.fa"), + target_bed=Path("/tmp/t.bed"), + ) + + assert captured["cmd"][:5] == [ + "/bin/hap.py", "/tmp/g.vcf", "/tmp/c.vcf", "-o", str(tmp_path / "happy") + ] + assert "-r" in captured["cmd"] and "/tmp/ref.fa" in captured["cmd"] + assert "-T" in captured["cmd"] and "/tmp/t.bed" in captured["cmd"] + + +def test_run_happy_omits_optional_flags_when_not_given(tmp_path, monkeypatch): + captured = {} + + def fake_run(cmd, **kw): + captured["cmd"] = cmd + Path(str(tmp_path / "happy") + ".vcf.gz").touch() + return SimpleNamespace(returncode=0, stdout="", stderr="") + + monkeypatch.setattr(subprocess, "run", fake_run) + run_happy( + happy_bin=Path("/bin/hap.py"), + golden_vcf=Path("/tmp/g.vcf"), + called_vcf=Path("/tmp/c.vcf"), + output_prefix=tmp_path / "happy", + ) + assert "-r" not in captured["cmd"] + assert "-T" not in captured["cmd"] + + +def test_run_happy_returns_output_vcf_path(tmp_path, monkeypatch): + expected = Path(str(tmp_path / "happy") + ".vcf.gz") + + def fake_run(cmd, **kw): + expected.touch() + return SimpleNamespace(returncode=0, stdout="", stderr="") + + monkeypatch.setattr(subprocess, "run", fake_run) + got = run_happy( + happy_bin=Path("/bin/hap.py"), + golden_vcf=Path("/tmp/g.vcf"), + called_vcf=Path("/tmp/c.vcf"), + output_prefix=tmp_path / "happy", + ) + assert got == expected + + +def test_run_happy_raises_when_returncode_nonzero(tmp_path, monkeypatch): + def fake_run(cmd, **kw): + return SimpleNamespace(returncode=1, stdout="", stderr="boom") + + monkeypatch.setattr(subprocess, "run", fake_run) + with pytest.raises(HappyExecutionError, match="exited 1"): + run_happy( + happy_bin=Path("/bin/hap.py"), + golden_vcf=Path("/tmp/g.vcf"), + called_vcf=Path("/tmp/c.vcf"), + output_prefix=tmp_path / "happy", + ) + + +def test_run_happy_raises_when_output_vcf_missing(tmp_path, monkeypatch): + """Even with returncode=0, missing the expected output VCF is a fatal contract violation.""" + def fake_run(cmd, **kw): + return SimpleNamespace(returncode=0, stdout="", stderr="") + + monkeypatch.setattr(subprocess, "run", fake_run) + with pytest.raises(HappyExecutionError, match="output VCF is missing"): + run_happy( + happy_bin=Path("/bin/hap.py"), + golden_vcf=Path("/tmp/g.vcf"), + called_vcf=Path("/tmp/c.vcf"), + output_prefix=tmp_path / "happy", + ) + + +def test_run_happy_passes_extra_args(tmp_path, monkeypatch): + captured = {} + + def fake_run(cmd, **kw): + captured["cmd"] = cmd + Path(str(tmp_path / "happy") + ".vcf.gz").touch() + return SimpleNamespace(returncode=0, stdout="", stderr="") + + monkeypatch.setattr(subprocess, "run", fake_run) + run_happy( + happy_bin=Path("/bin/hap.py"), + golden_vcf=Path("/tmp/g.vcf"), + called_vcf=Path("/tmp/c.vcf"), + output_prefix=tmp_path / "happy", + extra_args=["--engine", "vcfeval"], + ) + assert "--engine" in captured["cmd"] and "vcfeval" in captured["cmd"] + + +# --------------------------------------------------------------------------- +# _classify +# --------------------------------------------------------------------------- + +def test_classify_tp_when_truth_tp(): + assert _classify("TP", ".") == "TP" + + +def test_classify_tp_when_query_tp(): + assert _classify(".", "TP") == "TP" + + +def test_classify_tp_when_both_tp(): + assert _classify("TP", "TP") == "TP" + + +def test_classify_fn(): + assert _classify("FN", ".") == "FN" + + +def test_classify_fp(): + assert _classify(".", "FP") == "FP" + + +def test_classify_returns_none_for_nocall(): + assert _classify(".", ".") is None + + +def test_classify_returns_none_for_unknown_codes(): + assert _classify("N", "N") is None + + +def test_classify_tp_beats_other_signals(): + """If a record has TP on one side and FN on the other, hap.py considers it matched.""" + assert _classify("TP", "FN") == "TP" + assert _classify("FN", "TP") == "TP" + + +# --------------------------------------------------------------------------- +# parse_happy_output — synthetic VCFs +# --------------------------------------------------------------------------- + +_HAPPY_HEADER = """##fileformat=VCFv4.2 +##contig= +##FORMAT= +##FORMAT= +##FORMAT= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH\tQUERY +""" + + +def _make_happy_vcf(tmp_path: Path, lines: list[str]) -> Path: + """Write a tiny hap.py-shaped VCF and bgzip+index it for pysam.""" + raw = tmp_path / "happy.vcf" + raw.write_text(_HAPPY_HEADER + "\n".join(lines) + ("\n" if lines else "")) + bgz_path = pysam.tabix_index(str(raw), preset="vcf", force=True) + return Path(bgz_path) + + +def test_parse_happy_output_buckets_tp_fn_fp(tmp_path): + vcf = _make_happy_vcf(tmp_path, [ + "chr1\t100\t.\tA\tT\t.\t.\t.\tGT:BD:BVT\t1|0:TP:SNP\t1|0:TP:SNP", + "chr1\t200\t.\tA\tT\t.\t.\t.\tGT:BD:BVT\t1|0:FN:SNP\t.:.:.", + "chr1\t300\t.\tA\tT\t.\t.\t.\tGT:BD:BVT\t.:.:.\t1|0:FP:SNP", + "chr1\t400\t.\tA\tT\t.\t.\t.\tGT:BD:BVT\t.:.:.\t.:.:.", # no-call, dropped + ]) + buckets = parse_happy_output(vcf) + assert len(buckets["TP"]) == 1 + assert len(buckets["FN"]) == 1 + assert len(buckets["FP"]) == 1 + assert buckets["TP"][0].pos == 100 + assert buckets["FN"][0].pos == 200 + assert buckets["FP"][0].pos == 300 + + +def test_parse_happy_output_empty_file_returns_empty_buckets(tmp_path): + vcf = _make_happy_vcf(tmp_path, []) + buckets = parse_happy_output(vcf) + assert buckets == {"TP": [], "FN": [], "FP": []} + + +def test_parse_happy_output_raises_on_single_sample(tmp_path): + """A VCF with only TRUTH is not a hap.py output.""" + header = """##fileformat=VCFv4.2 +##contig= +##FORMAT= +##FORMAT= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH +chr1\t100\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:TP +""" + raw = tmp_path / "single.vcf" + raw.write_text(header) + bgz = Path(pysam.tabix_index(str(raw), preset="vcf", force=True)) + with pytest.raises(HappyParseError, match="TRUTH and QUERY"): + parse_happy_output(bgz) + + +def test_parse_happy_output_raises_when_BD_format_absent(tmp_path): + """A VCF lacking the BD FORMAT is not a hap.py output.""" + header = """##fileformat=VCFv4.2 +##contig= +##FORMAT= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH\tQUERY +chr1\t100\t.\tA\tT\t.\t.\t.\tGT\t1|0\t1|0 +""" + raw = tmp_path / "noBD.vcf" + raw.write_text(header) + bgz = Path(pysam.tabix_index(str(raw), preset="vcf", force=True)) + with pytest.raises(HappyParseError, match="BD FORMAT"): + parse_happy_output(bgz) + + +def test_parse_happy_output_preserves_record_chrom_pos(tmp_path): + vcf = _make_happy_vcf(tmp_path, [ + "chr1\t150\t.\tA\tT\t.\t.\t.\tGT:BD:BVT\t1|0:FN:SNP\t.:.:.", + "chr1\t250\t.\tA\tT\t.\t.\t.\tGT:BD:BVT\t1|0:FN:INDEL\t.:.:.", + ]) + fns = parse_happy_output(vcf)["FN"] + assert [(r.chrom, r.pos) for r in fns] == [("chr1", 150), ("chr1", 250)] diff --git a/tests/test_compare_vcfs/test_integration.py b/tests/test_compare_vcfs/test_integration.py new file mode 100644 index 00000000..a3d82112 --- /dev/null +++ b/tests/test_compare_vcfs/test_integration.py @@ -0,0 +1,159 @@ +""" +Real-hap.py integration test for `neat compare-vcfs` (issue #297). + +Runs an actual NEAT simulation, then invokes the real hap.py binary. Skipped +cleanly when hap.py isn't available — set NEAT_HAPPY_BIN to the absolute path +of a working hap.py to enable. + +Notes on hap.py packaging: the conda `bioconda::hap.py` package is Python-2-based. +Its shebang resolves `python` via PATH, so the test prepends the env's bin +directory so the child process picks up python2.7 from the same env. +""" +import gzip +import json +import os +import shutil +from pathlib import Path + +import pysam +import pytest + +from neat.compare_vcfs.runner import compare_vcfs_runner +from neat.read_simulator.runner import read_simulator_runner + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def happy_bin(): + """Skip the test if a working hap.py is unavailable.""" + explicit = os.environ.get("NEAT_HAPPY_BIN") + if explicit and Path(explicit).is_file(): + return Path(explicit) + found = shutil.which("hap.py") + if found: + return Path(found) + pytest.skip( + "hap.py is not available. Install via `conda create -n hap_py_env " + "-c bioconda -c conda-forge hap.py -y` and set NEAT_HAPPY_BIN to its " + "absolute path to enable this test." + ) + + +@pytest.fixture +def happy_env_path(happy_bin): + """Prepend hap.py's env bin to PATH so its #!/usr/bin/env python shebang + resolves to the python interpreter shipped alongside it.""" + return str(happy_bin.parent) + + +def _write_ref(path: Path) -> Path: + """A small but realistic reference: ~2kb with mixed GC content.""" + # Two contigs so we can also exercise per-contig stats. + seq_a = "ACGT" * 250 # 1000 bp + seq_b = ("AAAGGCCC" * 125) # 1000 bp, higher GC + path.write_text(f">chr1\n{seq_a}\n>chr2\n{seq_b}\n", encoding="utf-8") + # Index for hap.py (it needs a .fai) + pysam.FastaFile(str(path)) + return path + + +def _write_config(path: Path, ref_path: Path) -> Path: + cfg = ( + f"reference: {ref_path}\n" + "produce_fastq: false\n" + "produce_bam: false\n" + "produce_vcf: true\n" + "read_len: 100\n" + "coverage: 5\n" + "rng_seed: 42\n" + "mutation_rate: 0.01\n" + "overwrite_output: true\n" + "cleanup_splits: true\n" + ) + path.write_text(cfg, encoding="utf-8") + return path + + +def _called_vcf_dropping_first_variant(golden_vcf: Path, output_vcf: Path) -> Path: + """Read golden.vcf.gz, drop the first variant, write a new bgzipped+indexed VCF. + + Dropping a variant creates exactly one false negative; everything else is + a true positive (the rest of the truth's variants are also in the caller VCF). + """ + written = 0 + skipped = False + with pysam.VariantFile(str(golden_vcf)) as src: + # Keep the same samples / header + with pysam.VariantFile(str(output_vcf), "wz", header=src.header) as dst: + for rec in src: + if not skipped: + skipped = True + continue + dst.write(rec) + written += 1 + pysam.tabix_index(str(output_vcf), preset="vcf", force=True) + assert skipped, "golden VCF had no variants to drop; test setup wrong" + return output_vcf + + +# =========================================================================== +# Integration test — real hap.py end-to-end +# =========================================================================== + +def test_compare_vcfs_real_happy_end_to_end(tmp_path, happy_bin, happy_env_path, monkeypatch): + """ + Run NEAT, drop one variant from the golden VCF to create a caller VCF with + exactly one FN, then run compare-vcfs against the real hap.py binary. + """ + monkeypatch.setenv("PATH", happy_env_path + os.pathsep + os.environ["PATH"]) + + # 1. Simulate + sim_out = tmp_path / "sim_out" + sim_out.mkdir() + ref = _write_ref(tmp_path / "ref.fa") + cfg = _write_config(tmp_path / "conf.yml", ref) + read_simulator_runner(str(cfg), str(sim_out), "run") + + golden = sim_out / "run_golden.vcf.gz" + assert golden.is_file(), "NEAT did not produce a golden VCF" + assert (sim_out / "simulation_summary.json").is_file() + + # 2. Build a "called" VCF missing one variant + called = tmp_path / "called.vcf.gz" + _called_vcf_dropping_first_variant(golden, called) + + # 3. Run compare-vcfs against real hap.py + cmp_out = tmp_path / "cmp_out" + compare_vcfs_runner( + golden_vcf=str(golden), + called_vcf=str(called), + neat_run_dir=str(sim_out), + output_dir=str(cmp_out), + reference=str(ref), + happy_bin=str(happy_bin), + ) + + # 4. Assert the three reports exist and look sensible + assert (cmp_out / "comparison_summary.json").is_file() + assert (cmp_out / "comparison_summary.txt").is_file() + assert (cmp_out / "FN_with_reasons.vcf").is_file() + assert (cmp_out / "happy.vcf.gz").is_file() + + report = json.loads((cmp_out / "comparison_summary.json").read_text()) + assert report["schema_version"] == "1" + # We dropped exactly one variant — expect at least one FN + assert report["counts"]["FN"] >= 1 + assert report["counts"]["FP"] == 0 + assert report["metrics"]["precision"] == pytest.approx(1.0) + # Every FN must be attributed somehow + total_attribution = sum(report["fn_attribution"].values()) + assert total_attribution >= report["counts"]["FN"] + + # And the annotated VCF must carry NEAT_REASON + with pysam.VariantFile(str(cmp_out / "FN_with_reasons.vcf")) as vf: + assert "NEAT_REASON" in vf.header.info + for rec in vf: + assert "NEAT_REASON" in rec.info diff --git a/tests/test_compare_vcfs/test_reports.py b/tests/test_compare_vcfs/test_reports.py new file mode 100644 index 00000000..f466452b --- /dev/null +++ b/tests/test_compare_vcfs/test_reports.py @@ -0,0 +1,300 @@ +""" +Tests for neat/compare_vcfs/reports.py — metric computation, summary builders, +text rendering, and the annotated-FN VCF writer. +""" +import json +from pathlib import Path +from types import SimpleNamespace + +import pysam +import pytest + +from neat.compare_vcfs.attribution import ( + REASON_OUTSIDE_MUTATION_BED, + REASON_OUTSIDE_TARGET_BED, + REASON_UNKNOWN, +) +from neat.compare_vcfs.reports import ( + REPORT_SCHEMA_VERSION, + build_comparison_summary, + compute_metrics, + render_summary_txt, + summarize_fn_reasons, + write_comparison_summary_json, + write_comparison_summary_txt, + write_fn_attribution_plot, + write_fn_with_reasons, +) + + +# =========================================================================== +# compute_metrics +# =========================================================================== + +def test_compute_metrics_typical_counts(): + m = compute_metrics({"TP": 90, "FN": 10, "FP": 5}) + assert m["precision"] == pytest.approx(90 / 95) + assert m["recall"] == pytest.approx(0.9) + p, r = m["precision"], m["recall"] + assert m["f1"] == pytest.approx(2 * p * r / (p + r)) + + +def test_compute_metrics_perfect_call(): + m = compute_metrics({"TP": 100, "FN": 0, "FP": 0}) + assert m["precision"] == 1.0 + assert m["recall"] == 1.0 + assert m["f1"] == 1.0 + + +def test_compute_metrics_no_truth_no_calls_returns_none_metrics(): + """Empty buckets → precision/recall/f1 all undefined.""" + m = compute_metrics({"TP": 0, "FN": 0, "FP": 0}) + assert m == {"precision": None, "recall": None, "f1": None} + + +def test_compute_metrics_no_calls_means_precision_undefined(): + m = compute_metrics({"TP": 0, "FN": 5, "FP": 0}) + assert m["precision"] is None + assert m["recall"] == 0.0 # tp / (tp+fn) is defined (0 / 5) + assert m["f1"] is None + + +def test_compute_metrics_no_truth_means_recall_undefined(): + m = compute_metrics({"TP": 0, "FN": 0, "FP": 5}) + assert m["precision"] == 0.0 + assert m["recall"] is None + assert m["f1"] is None + + +# =========================================================================== +# summarize_fn_reasons +# =========================================================================== + +def test_summarize_fn_reasons_counts_each_tag(): + fn_reasons = [ + (SimpleNamespace(), [REASON_UNKNOWN]), + (SimpleNamespace(), [REASON_OUTSIDE_TARGET_BED]), + (SimpleNamespace(), [REASON_OUTSIDE_TARGET_BED, REASON_OUTSIDE_MUTATION_BED]), + ] + assert summarize_fn_reasons(fn_reasons) == { + REASON_UNKNOWN: 1, + REASON_OUTSIDE_TARGET_BED: 2, + REASON_OUTSIDE_MUTATION_BED: 1, + } + + +def test_summarize_fn_reasons_empty_input(): + assert summarize_fn_reasons([]) == {} + + +# =========================================================================== +# build_comparison_summary +# =========================================================================== + +_MISSING = object() + + +def _build_minimal_summary(tmp_path, counts=_MISSING, fn_attr=_MISSING): + if counts is _MISSING: + counts = {"TP": 1, "FN": 1, "FP": 1} + if fn_attr is _MISSING: + fn_attr = {REASON_UNKNOWN: 1} + return build_comparison_summary( + golden_vcf=tmp_path / "g.vcf", + called_vcf=tmp_path / "c.vcf", + neat_run_dir=tmp_path / "run", + simulation_summary_path=tmp_path / "run" / "simulation_summary.json", + happy_output_vcf=tmp_path / "out" / "happy.vcf.gz", + happy_output_prefix=tmp_path / "out" / "happy", + counts=counts, + fn_attribution=fn_attr, + fn_with_reasons_vcf=tmp_path / "out" / "FN_with_reasons.vcf", + comparison_summary_json=tmp_path / "out" / "comparison_summary.json", + comparison_summary_txt=tmp_path / "out" / "comparison_summary.txt", + ) + + +def test_build_comparison_summary_has_required_top_level_keys(tmp_path): + s = _build_minimal_summary(tmp_path) + assert set(s) == { + "schema_version", "neat_version", "generated_at", + "inputs", "happy", "counts", "metrics", "fn_attribution", "outputs", + } + + +def test_build_comparison_summary_schema_version_matches_constant(tmp_path): + s = _build_minimal_summary(tmp_path) + assert s["schema_version"] == REPORT_SCHEMA_VERSION + + +def test_build_comparison_summary_resolves_paths_absolute(tmp_path, monkeypatch): + monkeypatch.chdir(tmp_path) + s = build_comparison_summary( + golden_vcf=Path("g.vcf"), + called_vcf=Path("c.vcf"), + neat_run_dir=Path("run"), + simulation_summary_path=Path("run/simulation_summary.json"), + happy_output_vcf=Path("out/happy.vcf.gz"), + happy_output_prefix=Path("out/happy"), + counts={"TP": 0, "FN": 0, "FP": 0}, + fn_attribution={}, + fn_with_reasons_vcf=Path("out/FN_with_reasons.vcf"), + comparison_summary_json=Path("out/comparison_summary.json"), + comparison_summary_txt=Path("out/comparison_summary.txt"), + ) + for v in s["inputs"].values(): + assert Path(v).is_absolute() + for v in s["outputs"].values(): + assert Path(v).is_absolute() + + +def test_build_comparison_summary_carries_counts_and_metrics(tmp_path): + s = _build_minimal_summary(tmp_path, counts={"TP": 8, "FN": 2, "FP": 1}) + assert s["counts"] == {"TP": 8, "FN": 2, "FP": 1} + assert s["metrics"]["precision"] == pytest.approx(8 / 9) + assert s["metrics"]["recall"] == pytest.approx(0.8) + + +# =========================================================================== +# write_comparison_summary_json / write_comparison_summary_txt +# =========================================================================== + +def test_write_comparison_summary_json_round_trips(tmp_path): + s = _build_minimal_summary(tmp_path) + path = write_comparison_summary_json(s, tmp_path / "comparison_summary.json") + assert path.is_file() + parsed = json.loads(path.read_text()) + assert parsed["schema_version"] == REPORT_SCHEMA_VERSION + assert parsed["counts"] == {"TP": 1, "FN": 1, "FP": 1} + + +def test_write_comparison_summary_json_atomic(tmp_path): + """The temp file should not linger after success.""" + s = _build_minimal_summary(tmp_path) + write_comparison_summary_json(s, tmp_path / "report.json") + assert not (tmp_path / "report.json.tmp").exists() + + +def test_write_comparison_summary_txt_contains_expected_sections(tmp_path): + s = _build_minimal_summary(tmp_path) + path = write_comparison_summary_txt(s, tmp_path / "report.txt") + content = path.read_text() + for section in ["NEAT compare-vcfs report", "Inputs", "Classification", "Metrics", + "FN attribution", "Outputs"]: + assert section in content + + +def test_write_comparison_summary_txt_includes_counts_and_metrics(tmp_path): + s = _build_minimal_summary(tmp_path, counts={"TP": 90, "FN": 10, "FP": 5}) + txt = write_comparison_summary_txt(s, tmp_path / "report.txt").read_text() + assert "TP): 90" in txt + assert "FN): 10" in txt + assert "FP): 5" in txt + assert "0.9474" in txt # precision + assert "0.9000" in txt # recall + + +def test_render_summary_txt_handles_no_fns(tmp_path): + s = _build_minimal_summary(tmp_path, counts={"TP": 10, "FN": 0, "FP": 0}, fn_attr={}) + txt = render_summary_txt(s) + assert "(no false negatives)" in txt + + +def test_render_summary_txt_renders_NA_for_undefined_metrics(tmp_path): + s = _build_minimal_summary(tmp_path, counts={"TP": 0, "FN": 0, "FP": 0}, fn_attr={}) + txt = render_summary_txt(s) + assert "N/A" in txt + + +# =========================================================================== +# write_fn_with_reasons +# =========================================================================== + +_HAPPY_HEADER = """##fileformat=VCFv4.2 +##contig= +##FORMAT= +##FORMAT= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH\tQUERY +""" + + +def _make_source_vcf(tmp_path: Path, lines: list[str]) -> Path: + raw = tmp_path / "happy.vcf" + raw.write_text(_HAPPY_HEADER + "\n".join(lines) + ("\n" if lines else "")) + return Path(pysam.tabix_index(str(raw), preset="vcf", force=True)) + + +def test_write_fn_with_reasons_adds_neat_reason_info(tmp_path): + source = _make_source_vcf(tmp_path, [ + "chr1\t100\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:FN\t.:.", + "chr1\t200\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:FN\t.:.", + ]) + from neat.compare_vcfs.happy import parse_happy_output + fn_records = parse_happy_output(source)["FN"] + fn_reasons = [ + (fn_records[0], [REASON_UNKNOWN]), + (fn_records[1], [REASON_OUTSIDE_MUTATION_BED, REASON_OUTSIDE_TARGET_BED]), + ] + out = write_fn_with_reasons(source, fn_reasons, tmp_path / "FN_with_reasons.vcf") + assert out.is_file() + + # Re-open and verify the NEAT_REASON tag is present on each record + with pysam.VariantFile(str(out)) as vf: + assert "NEAT_REASON" in vf.header.info + records = list(vf) + assert len(records) == 2 + assert records[0].info["NEAT_REASON"] == (REASON_UNKNOWN,) + assert tuple(records[1].info["NEAT_REASON"]) == ( + REASON_OUTSIDE_MUTATION_BED, REASON_OUTSIDE_TARGET_BED, + ) + + +def test_write_fn_with_reasons_writes_only_fn_records(tmp_path): + """Source has TP/FN/FP; output should only contain the FN we passed in.""" + source = _make_source_vcf(tmp_path, [ + "chr1\t100\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:TP\t1|0:TP", + "chr1\t200\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:FN\t.:.", + "chr1\t300\t.\tA\tT\t.\t.\t.\tGT:BD\t.:.\t1|0:FP", + ]) + from neat.compare_vcfs.happy import parse_happy_output + fn_records = parse_happy_output(source)["FN"] + out = write_fn_with_reasons( + source, + [(fn_records[0], [REASON_UNKNOWN])], + tmp_path / "FN_with_reasons.vcf", + ) + with pysam.VariantFile(str(out)) as vf: + records = list(vf) + assert len(records) == 1 + assert records[0].pos == 200 + + +def test_write_fn_with_reasons_handles_empty_fn_list(tmp_path): + source = _make_source_vcf(tmp_path, [ + "chr1\t100\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:TP\t1|0:TP", + ]) + out = write_fn_with_reasons(source, [], tmp_path / "FN_with_reasons.vcf") + with pysam.VariantFile(str(out)) as vf: + assert "NEAT_REASON" in vf.header.info + assert list(vf) == [] + + +# =========================================================================== +# write_fn_attribution_plot +# =========================================================================== + +def test_write_fn_attribution_plot_writes_png(tmp_path): + path = write_fn_attribution_plot( + {REASON_OUTSIDE_MUTATION_BED: 12, REASON_OUTSIDE_TARGET_BED: 5, REASON_UNKNOWN: 3}, + tmp_path / "fn_attribution.png", + ) + assert path.is_file() + # PNG magic bytes + assert path.read_bytes()[:8] == b"\x89PNG\r\n\x1a\n" + + +def test_write_fn_attribution_plot_handles_no_fns(tmp_path): + """An empty attribution dict still produces a (placeholder) PNG.""" + path = write_fn_attribution_plot({}, tmp_path / "fn_attribution.png") + assert path.is_file() + assert path.read_bytes()[:8] == b"\x89PNG\r\n\x1a\n" diff --git a/tests/test_compare_vcfs/test_runner.py b/tests/test_compare_vcfs/test_runner.py new file mode 100644 index 00000000..a5ad9df4 --- /dev/null +++ b/tests/test_compare_vcfs/test_runner.py @@ -0,0 +1,271 @@ +""" +Tests for neat/compare_vcfs/runner.py — the scaffold for `neat compare-vcfs`. + +Covers input validation, hap.py discovery, simulation_summary.json loading, +and the scaffold's NotImplementedError contract. +""" +import json +import os +import stat +from pathlib import Path + +import pytest + +from neat.compare_vcfs.runner import ( + HappyNotFoundError, + SimulationSummaryError, + compare_vcfs_runner, + discover_happy, + load_simulation_summary, +) +from neat.read_simulator.utils.simulation_summary import SCHEMA_VERSION + + +# --------------------------------------------------------------------------- +# Synthetic-file builders +# --------------------------------------------------------------------------- + +def _touch(path: Path, content: str = "stub") -> Path: + """Write non-empty content; validate_input_path sys.exits on empty files.""" + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(content) + return path + + +def _make_summary_dir(tmp_path: Path, schema_version: str = SCHEMA_VERSION, + extra: dict | None = None) -> Path: + """Create a NEAT-run-dir with a valid (or version-bumped) simulation_summary.json.""" + run_dir = tmp_path / "neat_run" + run_dir.mkdir() + summary = { + "schema_version": schema_version, + "neat_version": "4.4.4", + "run": {}, + "config": {}, + "outputs": {}, + "delivered": { + "total_variants": 0, + "contigs_simulated": ["chr1"], + }, + } + if extra: + summary.update(extra) + (run_dir / "simulation_summary.json").write_text(json.dumps(summary)) + return run_dir + + +def _make_executable(path: Path) -> Path: + path.touch() + path.chmod(path.stat().st_mode | stat.S_IXUSR) + return path + + +# =========================================================================== +# discover_happy +# =========================================================================== + +def test_discover_happy_returns_path_from_PATH(tmp_path, monkeypatch): + fake = _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + assert discover_happy(None) == fake.resolve() + + +def test_discover_happy_raises_when_not_on_PATH(monkeypatch): + monkeypatch.setenv("PATH", "") + with pytest.raises(HappyNotFoundError, match="not found on \\$PATH"): + discover_happy(None) + + +def test_discover_happy_returns_explicit_path_when_file_exists(tmp_path): + fake = _make_executable(tmp_path / "my_happy") + assert discover_happy(str(fake)) == fake.resolve() + + +def test_discover_happy_raises_when_explicit_path_missing(tmp_path): + with pytest.raises(HappyNotFoundError, match="path does not exist"): + discover_happy(str(tmp_path / "nope")) + + +def test_discover_happy_explicit_path_wins_over_PATH(tmp_path, monkeypatch): + """An explicit --happy-bin should be used even if hap.py is also on PATH.""" + on_path = _make_executable(tmp_path / "hap.py") + explicit = _make_executable(tmp_path / "alt_hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + assert discover_happy(str(explicit)) == explicit.resolve() + assert discover_happy(str(explicit)) != on_path.resolve() + + +# =========================================================================== +# load_simulation_summary +# =========================================================================== + +def test_load_simulation_summary_returns_dict_for_valid_file(tmp_path): + run_dir = _make_summary_dir(tmp_path) + data = load_simulation_summary(run_dir) + assert data["schema_version"] == SCHEMA_VERSION + assert data["neat_version"] == "4.4.4" + + +def test_load_simulation_summary_raises_when_missing(tmp_path): + run_dir = tmp_path / "empty_run" + run_dir.mkdir() + with pytest.raises(SimulationSummaryError, match="not found"): + load_simulation_summary(run_dir) + + +def test_load_simulation_summary_raises_on_malformed_json(tmp_path): + run_dir = tmp_path / "broken" + run_dir.mkdir() + (run_dir / "simulation_summary.json").write_text("{not valid json") + with pytest.raises(SimulationSummaryError, match="not valid JSON"): + load_simulation_summary(run_dir) + + +def test_load_simulation_summary_raises_on_version_mismatch(tmp_path): + run_dir = _make_summary_dir(tmp_path, schema_version="99") + with pytest.raises(SimulationSummaryError, match="schema_version"): + load_simulation_summary(run_dir) + + +def test_load_simulation_summary_raises_on_missing_schema_version(tmp_path): + """A summary without schema_version is treated as a version mismatch (None ≠ '1').""" + run_dir = tmp_path / "no_version" + run_dir.mkdir() + (run_dir / "simulation_summary.json").write_text(json.dumps({"foo": "bar"})) + with pytest.raises(SimulationSummaryError, match="schema_version"): + load_simulation_summary(run_dir) + + +# =========================================================================== +# compare_vcfs_runner — input validation +# =========================================================================== + +def test_runner_exits_when_golden_missing(tmp_path): + """Missing truth VCF: validate_input_path sys.exit(5).""" + run_dir = _make_summary_dir(tmp_path) + called = _touch(tmp_path / "called.vcf") + with pytest.raises(SystemExit) as excinfo: + compare_vcfs_runner( + golden_vcf=str(tmp_path / "nope.vcf"), + called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(tmp_path / "out"), + ) + assert excinfo.value.code == 5 + + +def test_runner_exits_when_called_missing(tmp_path): + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + with pytest.raises(SystemExit) as excinfo: + compare_vcfs_runner( + golden_vcf=str(golden), + called_vcf=str(tmp_path / "nope.vcf"), + neat_run_dir=str(run_dir), + output_dir=str(tmp_path / "out"), + ) + assert excinfo.value.code == 5 + + +def test_runner_raises_when_run_dir_missing(tmp_path): + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + with pytest.raises(FileNotFoundError, match="--neat-run-dir"): + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(tmp_path / "nope_dir"), + output_dir=str(tmp_path / "out"), + ) + + +def test_runner_raises_when_happy_missing(tmp_path, monkeypatch): + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + monkeypatch.setenv("PATH", "") # nothing on PATH + with pytest.raises(HappyNotFoundError): + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(tmp_path / "out"), + ) + + +def test_runner_completes_end_to_end_writing_all_reports(tmp_path, monkeypatch): + """End-to-end happy path: mocked hap.py writes a tiny output VCF; runner + produces all three reports and does not raise.""" + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + out_dir = tmp_path / "out" + + from neat.compare_vcfs import runner as runner_mod + + happy_header = """##fileformat=VCFv4.2 +##contig= +##FORMAT= +##FORMAT= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH\tQUERY +chr1\t100\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:TP\t1|0:TP +chr1\t200\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:FN\t.:. +chr1\t300\t.\tA\tT\t.\t.\t.\tGT:BD\t.:.\t1|0:FP +""" + + def fake_run_happy(happy_bin, golden_vcf, called_vcf, output_prefix, **kw): + import pysam + raw = Path(str(output_prefix) + ".vcf") + raw.write_text(happy_header) + bgz = Path(pysam.tabix_index(str(raw), preset="vcf", force=True)) + return bgz + + monkeypatch.setattr(runner_mod, "run_happy", fake_run_happy) + + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(out_dir), + ) + + assert (out_dir / "comparison_summary.json").is_file() + assert (out_dir / "comparison_summary.txt").is_file() + assert (out_dir / "FN_with_reasons.vcf").is_file() + + import json as _json + report = _json.loads((out_dir / "comparison_summary.json").read_text()) + assert report["counts"] == {"TP": 1, "FN": 1, "FP": 1} + assert report["metrics"]["precision"] == pytest.approx(0.5) + assert report["metrics"]["recall"] == pytest.approx(0.5) + assert report["metrics"]["f1"] == pytest.approx(0.5) + # chr1 IS simulated and no beds configured → FN attributed as 'unknown' + assert report["fn_attribution"] == {"unknown": 1} + + +def test_runner_validates_optional_reference_when_provided(tmp_path): + """A passed --reference that doesn't exist should fail fast via validate_input_path.""" + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + with pytest.raises(SystemExit) as excinfo: + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(tmp_path / "out"), + reference=str(tmp_path / "missing_ref.fa"), + ) + assert excinfo.value.code == 5 + + +def test_runner_validates_optional_target_bed_when_provided(tmp_path): + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + with pytest.raises(SystemExit) as excinfo: + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(tmp_path / "out"), + target_bed=str(tmp_path / "missing_target.bed"), + ) + assert excinfo.value.code == 5 diff --git a/tests/test_read_simulator/test_simulation_summary.py b/tests/test_read_simulator/test_simulation_summary.py new file mode 100644 index 00000000..38cc216f --- /dev/null +++ b/tests/test_read_simulator/test_simulation_summary.py @@ -0,0 +1,348 @@ +""" +Tests for neat/read_simulator/utils/simulation_summary.py + +Covers the write_simulation_summary() helper that emits the per-run JSON manifest +consumed by `neat compare-vcfs`, plus the private count helpers. +""" +import json +import re +from datetime import datetime, timezone +from pathlib import Path +from types import SimpleNamespace + +import pysam +import pytest + +from neat.read_simulator.utils.simulation_summary import ( + SCHEMA_VERSION, + _abs_or_none, + _count_reads, + _count_variants, + _iso_utc, + write_simulation_summary, +) + + +# --------------------------------------------------------------------------- +# Synthetic-file builders +# --------------------------------------------------------------------------- + +_VCF_HEADER = ( + "##fileformat=VCFv4.2\n" + "##contig=\n" + "##contig=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" +) + + +def _make_vcf(tmp_path: Path, records: list[tuple[str, int]]) -> Path: + """records is a list of (chrom, pos) tuples; alt is always 'T' over ref 'A'.""" + path = tmp_path / "variants.vcf" + with open(path, "w") as fh: + fh.write(_VCF_HEADER) + for chrom, pos in records: + fh.write(f"{chrom}\t{pos}\t.\tA\tT\t.\t.\t.\n") + return path + + +def _make_bam(tmp_path: Path, n_reads: int) -> Path: + """Write a tiny BAM with n_reads identical aligned reads on chr1.""" + path = tmp_path / "reads.bam" + header = {"HD": {"VN": "1.6"}, "SQ": [{"SN": "chr1", "LN": 1000}]} + with pysam.AlignmentFile(str(path), "wb", header=header) as bf: + for i in range(n_reads): + a = pysam.AlignedSegment() + a.query_name = f"r{i}" + a.query_sequence = "ACGT" + a.flag = 0 + a.reference_id = 0 + a.reference_start = 100 + i + a.mapping_quality = 60 + a.cigar = ((0, 4),) + a.query_qualities = pysam.qualitystring_to_array("IIII") + bf.write(a) + return path + + +def _make_fastq(tmp_path: Path, n_records: int, name: str = "r1.fastq") -> Path: + path = tmp_path / name + with open(path, "w") as fh: + for i in range(n_records): + fh.write(f"@read{i}\nACGT\n+\nIIII\n") + return path + + +def _make_options(**overrides) -> SimpleNamespace: + """Build a stub Options that exposes the attrs write_simulation_summary reads.""" + defaults = dict( + reference=None, + coverage=30, + read_len=150, + paired_ended=False, + fragment_mean=None, fragment_st_dev=None, + ploidy=2, + rng_seed=42, + threads=4, + mutation_rate=None, + mutation_bed=None, target_bed=None, discard_bed=None, include_vcf=None, + mutation_model=None, gc_model=None, error_model=None, fragment_model=None, + fq1=None, fq2=None, bam=None, vcf=None, + ) + defaults.update(overrides) + return SimpleNamespace(**defaults) + + +# =========================================================================== +# _abs_or_none +# =========================================================================== + +def test_abs_or_none_returns_none_for_none(): + assert _abs_or_none(None) is None + + +def test_abs_or_none_resolves_relative_path(tmp_path, monkeypatch): + monkeypatch.chdir(tmp_path) + result = _abs_or_none("foo/bar.txt") + assert Path(result).is_absolute() + assert result.endswith("foo/bar.txt") + + +def test_abs_or_none_preserves_absolute_path(tmp_path): + target = tmp_path / "x.txt" + assert _abs_or_none(target) == str(target.resolve()) + + +# =========================================================================== +# _iso_utc +# =========================================================================== + +def test_iso_utc_format_ends_in_Z(): + """Schema requires an ISO-8601 UTC timestamp with the 'Z' suffix.""" + s = _iso_utc(0) + assert s.endswith("Z") + assert re.match(r"^\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}Z$", s) + + +def test_iso_utc_epoch_zero_is_1970(): + assert _iso_utc(0) == "1970-01-01T00:00:00Z" + + +def test_iso_utc_round_trips_through_datetime(): + epoch = 1_700_000_000.0 + parsed = datetime.fromisoformat(_iso_utc(epoch).replace("Z", "+00:00")) + assert parsed == datetime.fromtimestamp(epoch, tz=timezone.utc).replace(microsecond=0) + + +# =========================================================================== +# _count_variants +# =========================================================================== + +def test_count_variants_returns_none_pair_for_none_path(): + assert _count_variants(None) == (None, None) + + +def test_count_variants_returns_none_pair_for_missing_file(tmp_path): + assert _count_variants(str(tmp_path / "does_not_exist.vcf")) == (None, None) + + +def test_count_variants_counts_and_groups_by_contig(tmp_path): + vcf = _make_vcf(tmp_path, [("chr1", 100), ("chr1", 200), ("chr2", 50)]) + total, by_contig = _count_variants(str(vcf)) + assert total == 3 + assert by_contig == {"chr1": 2, "chr2": 1} + + +def test_count_variants_empty_vcf_returns_zero(tmp_path): + vcf = _make_vcf(tmp_path, []) + total, by_contig = _count_variants(str(vcf)) + assert total == 0 + assert by_contig == {} + + +def test_count_variants_malformed_vcf_returns_none_pair(tmp_path, caplog): + bad = tmp_path / "bad.vcf" + bad.write_text("not a real vcf at all\n") + result = _count_variants(str(bad)) + assert result == (None, None) + + +# =========================================================================== +# _count_reads +# =========================================================================== + +def test_count_reads_returns_none_when_no_inputs(): + assert _count_reads(None, None, None, paired_ended=False) is None + + +def test_count_reads_prefers_bam_over_fastq(tmp_path): + bam = _make_bam(tmp_path, 7) + fq = _make_fastq(tmp_path, 999) + assert _count_reads(str(bam), str(fq), None, paired_ended=False) == 7 + + +def test_count_reads_falls_back_to_fastq_when_no_bam(tmp_path): + fq = _make_fastq(tmp_path, 12) + assert _count_reads(None, str(fq), None, paired_ended=False) == 12 + + +def test_count_reads_doubles_for_paired_ended_fastq(tmp_path): + fq = _make_fastq(tmp_path, 10, name="r1.fastq") + assert _count_reads(None, str(fq), None, paired_ended=True) == 20 + + +def test_count_reads_returns_none_for_missing_bam_and_missing_fq(tmp_path): + assert _count_reads( + str(tmp_path / "nope.bam"), str(tmp_path / "nope.fq"), None, paired_ended=False + ) is None + + +# =========================================================================== +# write_simulation_summary — schema shape and content +# =========================================================================== + +def test_write_simulation_summary_writes_file_in_output_dir(tmp_path): + options = _make_options() + path = write_simulation_summary( + options=options, output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, + contigs_simulated=[], + ) + assert path == tmp_path / "simulation_summary.json" + assert path.is_file() + + +def test_write_simulation_summary_has_expected_top_level_keys(tmp_path): + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + summary = json.loads((tmp_path / "simulation_summary.json").read_text()) + assert set(summary) == {"schema_version", "neat_version", "run", "config", "outputs", "delivered"} + + +def test_write_simulation_summary_schema_version_matches_constant(tmp_path): + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + summary = json.loads((tmp_path / "simulation_summary.json").read_text()) + assert summary["schema_version"] == SCHEMA_VERSION + + +def test_write_simulation_summary_neat_version_is_populated(tmp_path): + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + summary = json.loads((tmp_path / "simulation_summary.json").read_text()) + assert isinstance(summary["neat_version"], str) and summary["neat_version"] + + +def test_write_simulation_summary_echoes_config_fields(tmp_path): + options = _make_options(coverage=42, read_len=200, paired_ended=True, ploidy=4, rng_seed=99, threads=8) + write_simulation_summary( + options=options, output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + cfg = json.loads((tmp_path / "simulation_summary.json").read_text())["config"] + assert cfg["coverage"] == 42 + assert cfg["read_len"] == 200 + assert cfg["paired_ended"] is True + assert cfg["ploidy"] == 4 + assert cfg["rng_seed"] == 99 + assert cfg["threads"] == 8 + + +def test_write_simulation_summary_optional_config_fields_default_to_null(tmp_path): + """All Path-typed optional fields and rate must serialize as JSON null when absent.""" + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + cfg = json.loads((tmp_path / "simulation_summary.json").read_text())["config"] + optional_keys = [ + "reference", "mutation_rate", "mutation_bed", "target_bed", "discard_bed", + "include_vcf", "mutation_model", "gc_model", "error_model", "fragment_model", + ] + for k in optional_keys: + assert cfg[k] is None, f"expected {k} to be null" + + +def test_write_simulation_summary_resolves_paths_to_absolute(tmp_path, monkeypatch): + """Relative inputs (cwd, reference, config_path) must serialize as absolute paths.""" + monkeypatch.chdir(tmp_path) + ref = tmp_path / "ref.fa" + ref.touch() + out_rel = tmp_path / "out_rel" + out_rel.mkdir() # runner creates output_dir before calling the helper + options = _make_options(reference="ref.fa") # relative + write_simulation_summary( + options=options, output_dir="out_rel", + file_prefix="run", config_path="cfg.yml", analysis_start=0.0, + contigs_simulated=[], + ) + summary = json.loads((out_rel / "simulation_summary.json").read_text()) + assert Path(summary["run"]["output_dir"]).is_absolute() + assert Path(summary["run"]["config_file"]).is_absolute() + assert Path(summary["config"]["reference"]).is_absolute() + + +def test_write_simulation_summary_contigs_simulated_preserves_order(tmp_path): + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, + contigs_simulated=["chr3", "chr1", "chr2"], + ) + summary = json.loads((tmp_path / "simulation_summary.json").read_text()) + assert summary["delivered"]["contigs_simulated"] == ["chr3", "chr1", "chr2"] + + +def test_write_simulation_summary_outputs_none_when_nothing_produced(tmp_path): + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + outputs = json.loads((tmp_path / "simulation_summary.json").read_text())["outputs"] + assert outputs == {"fastq": None, "bam": None, "vcf": None} + + +def test_write_simulation_summary_outputs_populated_for_multi_output_run(tmp_path): + vcf = _make_vcf(tmp_path, [("chr1", 100), ("chr1", 200), ("chr2", 50)]) + bam = _make_bam(tmp_path, 5) + fq1 = _make_fastq(tmp_path, 3, "r1.fastq") + fq2 = _make_fastq(tmp_path, 3, "r2.fastq") + options = _make_options(paired_ended=True, vcf=vcf, bam=bam, fq1=fq1, fq2=fq2) + write_simulation_summary( + options=options, output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, + contigs_simulated=["chr1", "chr2"], + ) + summary = json.loads((tmp_path / "simulation_summary.json").read_text()) + assert summary["outputs"]["vcf"] == str(vcf.resolve()) + assert summary["outputs"]["bam"] == str(bam.resolve()) + assert summary["outputs"]["fastq"] == [str(fq1.resolve()), str(fq2.resolve())] + assert summary["delivered"]["total_variants"] == 3 + assert summary["delivered"]["variants_by_contig"] == {"chr1": 2, "chr2": 1} + assert summary["delivered"]["total_reads"] == 5 # BAM preferred over FASTQ + + +def test_write_simulation_summary_duration_matches_elapsed(tmp_path): + """duration_seconds should be the difference between now and analysis_start.""" + import time + start = time.time() - 10 # ten seconds ago + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=start, contigs_simulated=[], + ) + summary = json.loads((tmp_path / "simulation_summary.json").read_text()) + assert 9.5 <= summary["run"]["duration_seconds"] <= 11.0 + + +def test_write_simulation_summary_no_tmp_file_remains(tmp_path): + """Atomic write contract: .json.tmp must not linger after success.""" + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, contigs_simulated=[], + ) + assert not (tmp_path / "simulation_summary.json.tmp").exists() + assert (tmp_path / "simulation_summary.json").exists() From 6961af3768cfe53f1a74a692adf2118437d13883 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 20:44:41 -0500 Subject: [PATCH 2/4] Cover compare-vcfs CLI, --plot, stale-bed, multi-FN gaps Audit found five untested paths in the compare-vcfs work: - CLI layer: subcommand registration + argument routing - --plot flag: runner integration (on/off) - Stale BED path in simulation_summary: undefined failure mode - Chrom-prefix mismatch ('1' vs 'chr1'): documents a user gotcha - Multi-FN scaling through real hap.py Adds 9 tests across 4 files (+1 new test_cli.py). Full suite: 727 passing. Co-Authored-By: Claude Opus 4.7 (1M context) --- tests/test_compare_vcfs/test_attribution.py | 15 ++++ tests/test_compare_vcfs/test_cli.py | 97 ++++++++++++++++++++ tests/test_compare_vcfs/test_integration.py | 60 +++++++++++++ tests/test_compare_vcfs/test_runner.py | 98 +++++++++++++++++++++ 4 files changed, 270 insertions(+) create mode 100644 tests/test_compare_vcfs/test_cli.py diff --git a/tests/test_compare_vcfs/test_attribution.py b/tests/test_compare_vcfs/test_attribution.py index 47b2a9c3..589102ad 100644 --- a/tests/test_compare_vcfs/test_attribution.py +++ b/tests/test_compare_vcfs/test_attribution.py @@ -183,6 +183,21 @@ def test_attribute_fn_unknown_when_inside_all_configured_beds(): assert reasons == [REASON_UNKNOWN] +def test_attribute_fn_chrom_prefix_mismatch_treats_as_outside(): + """ + Known gotcha: if a BED uses '1' but the run records 'chr1' as simulated + (or vice versa), the unmatched chrom name means the FN tagged 'outside' + the BED. Documents that users must align chrom conventions across + reference, BED, golden VCF, and caller VCF — NEAT does not normalize. + """ + reasons = attribute_fn( + "chr1", 100, frozenset({"chr1"}), + mutation_intervals={"1": [(0, 1000)]}, # bare-numeric chrom in BED + target_intervals=None, + ) + assert reasons == [REASON_OUTSIDE_MUTATION_BED] + + def test_attribute_fn_bed_missing_chrom_is_outside(): """A bed that doesn't mention this chrom counts as 'outside' for it.""" reasons = attribute_fn( diff --git a/tests/test_compare_vcfs/test_cli.py b/tests/test_compare_vcfs/test_cli.py new file mode 100644 index 00000000..341d4baf --- /dev/null +++ b/tests/test_compare_vcfs/test_cli.py @@ -0,0 +1,97 @@ +""" +CLI-layer tests for `neat compare-vcfs` — subcommand registration, argument +parsing, and routing into the runner. +""" +import argparse + +import pytest + +from neat.cli.cli import Cli +from neat.cli.commands.compare_vcfs import Command + + +def _compare_vcfs_subparser(): + cli = Cli() + subparsers_action = next( + a for a in cli.parser._actions if isinstance(a, argparse._SubParsersAction) + ) + return subparsers_action._name_parser_map.get("compare-vcfs") + + +# =========================================================================== +# Subcommand registration +# =========================================================================== + +def test_compare_vcfs_subcommand_is_registered(): + """`neat compare-vcfs` must be discovered by Cli's pkgutil-based scan.""" + sub = _compare_vcfs_subparser() + assert sub is not None, "compare-vcfs subcommand is not registered" + + +def test_compare_vcfs_help_text_lists_all_documented_flags(): + """Regression guard: every flag described in the docs must appear in --help.""" + sub = _compare_vcfs_subparser() + help_text = sub.format_help() + for flag in [ + "--neat-run-dir", + "--output-dir", + "--reference", + "--target-bed", + "--happy-bin", + "--plot", + ]: + assert flag in help_text, f"--help missing {flag}" + + +# =========================================================================== +# Argument routing — Command.execute → compare_vcfs_runner +# =========================================================================== + +def test_command_routes_all_args_to_runner(monkeypatch): + """Parsed args must reach compare_vcfs_runner with the right kwargs.""" + captured = {} + + def fake_runner(**kwargs): + captured.update(kwargs) + + from neat.cli.commands import compare_vcfs as cmd_mod + monkeypatch.setattr(cmd_mod, "compare_vcfs_runner", fake_runner) + + parser = argparse.ArgumentParser() + cmd = Command(parser) + args = parser.parse_args([ + "/g.vcf", "/c.vcf", + "--neat-run-dir", "/run", + "--output-dir", "/out", + "--reference", "/ref.fa", + "--target-bed", "/t.bed", + "--happy-bin", "/bin/hap.py", + "--plot", + ]) + cmd.execute(args) + + assert captured == { + "golden_vcf": "/g.vcf", + "called_vcf": "/c.vcf", + "neat_run_dir": "/run", + "output_dir": "/out", + "reference": "/ref.fa", + "target_bed": "/t.bed", + "happy_bin": "/bin/hap.py", + "plot": True, + } + + +def test_command_optional_flags_default_to_none_or_false(): + """When optional flags aren't supplied, parsed values match runner defaults.""" + parser = argparse.ArgumentParser() + Command(parser) + args = parser.parse_args([ + "/g.vcf", "/c.vcf", + "--neat-run-dir", "/run", + "--output-dir", "/out", + ]) + assert args.reference is None + assert args.target_bed is None + assert args.happy_bin is None + assert args.plot is False diff --git a/tests/test_compare_vcfs/test_integration.py b/tests/test_compare_vcfs/test_integration.py index a3d82112..08c1f555 100644 --- a/tests/test_compare_vcfs/test_integration.py +++ b/tests/test_compare_vcfs/test_integration.py @@ -99,6 +99,21 @@ def _called_vcf_dropping_first_variant(golden_vcf: Path, output_vcf: Path) -> Pa return output_vcf +def _called_vcf_dropping_first_n(golden_vcf: Path, output_vcf: Path, n: int) -> int: + """Drop the first n variants from golden_vcf into output_vcf; returns the + number actually dropped (capped at the total available).""" + actually_dropped = 0 + with pysam.VariantFile(str(golden_vcf)) as src: + with pysam.VariantFile(str(output_vcf), "wz", header=src.header) as dst: + for rec in src: + if actually_dropped < n: + actually_dropped += 1 + continue + dst.write(rec) + pysam.tabix_index(str(output_vcf), preset="vcf", force=True) + return actually_dropped + + # =========================================================================== # Integration test — real hap.py end-to-end # =========================================================================== @@ -157,3 +172,48 @@ def test_compare_vcfs_real_happy_end_to_end(tmp_path, happy_bin, happy_env_path, assert "NEAT_REASON" in vf.header.info for rec in vf: assert "NEAT_REASON" in rec.info + + +def test_compare_vcfs_real_happy_multi_fn(tmp_path, happy_bin, happy_env_path, monkeypatch): + """ + Drop multiple variants to verify FN counts scale linearly and that each + FN gets a NEAT_REASON tag (locks in attribution behavior at scale). + """ + monkeypatch.setenv("PATH", happy_env_path + os.pathsep + os.environ["PATH"]) + + sim_out = tmp_path / "sim_out" + sim_out.mkdir() + ref = _write_ref(tmp_path / "ref.fa") + cfg = _write_config(tmp_path / "conf.yml", ref) + read_simulator_runner(str(cfg), str(sim_out), "run") + golden = sim_out / "run_golden.vcf.gz" + + called = tmp_path / "called.vcf.gz" + n_dropped = _called_vcf_dropping_first_n(golden, called, n=3) + assert n_dropped == 3, "golden produced fewer than 3 variants — config too small" + + cmp_out = tmp_path / "cmp_out" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(sim_out), + output_dir=str(cmp_out), + reference=str(ref), + happy_bin=str(happy_bin), + plot=True, + ) + + report = json.loads((cmp_out / "comparison_summary.json").read_text()) + assert report["counts"]["FN"] >= 3, "expected ≥3 FNs after dropping 3 variants" + assert report["counts"]["FP"] == 0 + assert sum(report["fn_attribution"].values()) >= report["counts"]["FN"] + + # Every FN record must be annotated; this is the contract step 6 promises + with pysam.VariantFile(str(cmp_out / "FN_with_reasons.vcf")) as vf: + fn_records = list(vf) + assert len(fn_records) == report["counts"]["FN"] + for rec in fn_records: + tags = rec.info["NEAT_REASON"] + assert tags, f"FN at {rec.chrom}:{rec.pos} has empty NEAT_REASON" + + # --plot was on, so the bar chart should be present + assert (cmp_out / "fn_attribution.png").is_file() diff --git a/tests/test_compare_vcfs/test_runner.py b/tests/test_compare_vcfs/test_runner.py index a5ad9df4..823ee924 100644 --- a/tests/test_compare_vcfs/test_runner.py +++ b/tests/test_compare_vcfs/test_runner.py @@ -257,6 +257,104 @@ def test_runner_validates_optional_reference_when_provided(tmp_path): assert excinfo.value.code == 5 +# =========================================================================== +# Shared fake-hap.py helpers for plot / stale-bed tests +# =========================================================================== + +_HAPPY_HEADER_ONE_FN = """##fileformat=VCFv4.2 +##contig= +##FORMAT= +##FORMAT= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH\tQUERY +chr1\t200\t.\tA\tT\t.\t.\t.\tGT:BD\t1|0:FN\t.:. +""" + + +def _install_fake_happy(monkeypatch, body: str): + """Patch run_happy to write `body` as the hap.py output VCF.""" + from neat.compare_vcfs import runner as runner_mod + import pysam + + def fake_run_happy(happy_bin, golden_vcf, called_vcf, output_prefix, **kw): + raw = Path(str(output_prefix) + ".vcf") + raw.write_text(body) + return Path(pysam.tabix_index(str(raw), preset="vcf", force=True)) + + monkeypatch.setattr(runner_mod, "run_happy", fake_run_happy) + + +# =========================================================================== +# --plot wiring +# =========================================================================== + +def test_runner_writes_plot_when_enabled(tmp_path, monkeypatch): + """plot=True must produce fn_attribution.png alongside the standard reports.""" + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + _install_fake_happy(monkeypatch, _HAPPY_HEADER_ONE_FN) + + out_dir = tmp_path / "out" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(out_dir), + plot=True, + ) + plot_path = out_dir / "fn_attribution.png" + assert plot_path.is_file() + assert plot_path.read_bytes()[:8] == b"\x89PNG\r\n\x1a\n" + + +def test_runner_omits_plot_by_default(tmp_path, monkeypatch): + """plot defaults to False; fn_attribution.png must NOT be written.""" + run_dir = _make_summary_dir(tmp_path) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + _install_fake_happy(monkeypatch, _HAPPY_HEADER_ONE_FN) + + out_dir = tmp_path / "out" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(out_dir), + ) + assert not (out_dir / "fn_attribution.png").exists() + + +# =========================================================================== +# Stale BED path in simulation_summary +# =========================================================================== + +def test_runner_fails_clearly_when_summary_bed_path_no_longer_exists(tmp_path, monkeypatch): + """If simulation_summary references a BED that's been moved/deleted, the + failure must surface (currently as FileNotFoundError from load_bed_intervals). + Locking in current behavior so a regression to a silent skip is caught.""" + bed = tmp_path / "stale.bed" + bed.write_text("chr1\t0\t1000\n") + run_dir = _make_summary_dir(tmp_path, extra={ + "config": {"mutation_bed": str(bed), "target_bed": None}, + }) + bed.unlink() + + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + _install_fake_happy(monkeypatch, _HAPPY_HEADER_ONE_FN) + + with pytest.raises(FileNotFoundError): + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(tmp_path / "out"), + ) + + def test_runner_validates_optional_target_bed_when_provided(tmp_path): run_dir = _make_summary_dir(tmp_path) golden = _touch(tmp_path / "golden.vcf") From 73cd132853b1bdb5db41b2a51c331196a9dcec56 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 22:33:58 -0500 Subject: [PATCH 3/4] Detect chrom-naming mismatches in compare-vcfs; add --chrom-aliases MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When a mutation_bed or target_bed uses different chrom names than the reference (e.g., '1' vs 'chr1', 'MT' vs 'chrM'), NEAT previously produced silent garbage — every FN tagged 'outside_mutation_bed' even though they were actually inside. Changes: - new neat/common/chrom_names.py with prefix/mitochondrial heuristics and a user-supplied alias-file parser - simulation_summary.json now records delivered.reference_contigs (full FASTA contig set, separate from contigs_simulated) - compare_vcfs_runner runs detect_chrom_naming_mismatches before attribution; emits WARNING logs AND a 'warnings' array in comparison_summary.json with suggested aliases - new --chrom-aliases TSV flag applies user-supplied mappings to BED chrom names at load time - README + ChangeLog updated; mismatch behavior documented as non- normalizing by default NEAT does not auto-normalize: silent renaming masks real bugs. Detection + warnings + opt-in normalization is the policy. 34 new tests (chrom_names utility, schema, attribution detection, runner integration). Full suite: 761 passing. Co-Authored-By: Claude Opus 4.7 (1M context) --- ChangeLog.md | 12 +- README.md | 12 +- neat/cli/commands/compare_vcfs.py | 9 ++ neat/common/__init__.py | 1 + neat/common/chrom_names.py | 106 +++++++++++++ neat/compare_vcfs/attribution.py | 89 ++++++++++- neat/compare_vcfs/reports.py | 8 + neat/compare_vcfs/runner.py | 16 +- neat/read_simulator/runner.py | 1 + .../utils/simulation_summary.py | 2 + tests/test_common/__init__.py | 0 tests/test_common/test_chrom_names.py | 144 ++++++++++++++++++ tests/test_compare_vcfs/test_attribution.py | 112 ++++++++++++++ tests/test_compare_vcfs/test_cli.py | 4 + tests/test_compare_vcfs/test_reports.py | 3 +- tests/test_compare_vcfs/test_runner.py | 62 ++++++++ .../test_simulation_summary.py | 26 ++++ 17 files changed, 597 insertions(+), 10 deletions(-) create mode 100644 neat/common/chrom_names.py create mode 100644 tests/test_common/__init__.py create mode 100644 tests/test_common/test_chrom_names.py diff --git a/ChangeLog.md b/ChangeLog.md index 12c0c8d2..5cc4f683 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -29,9 +29,19 @@ 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. +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 diff --git a/README.md b/README.md index c23354b3..f67d670d 100755 --- a/README.md +++ b/README.md @@ -619,7 +619,8 @@ neat compare-vcfs golden.vcf called.vcf \ --reference reference.fa \ [--target-bed target.bed] \ [--happy-bin /abs/path/to/hap.py] \ - [--plot] + [--plot] \ + [--chrom-aliases aliases.tsv] ``` Outputs (in `--output-dir`): @@ -632,6 +633,15 @@ Outputs (in `--output-dir`): | `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_namecanonical_name`) to apply the +rename at load time. This also handles common mitochondrial variants +(`M`/`MT`/`chrM`/`chrMT`). + **False-negative reason categories:** | Tag | Meaning | diff --git a/neat/cli/commands/compare_vcfs.py b/neat/cli/commands/compare_vcfs.py index 994209a8..90376206 100644 --- a/neat/cli/commands/compare_vcfs.py +++ b/neat/cli/commands/compare_vcfs.py @@ -64,6 +64,14 @@ def add_arguments(self, parser: argparse.ArgumentParser): 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( @@ -75,4 +83,5 @@ def execute(self, arguments: argparse.Namespace): target_bed=arguments.target_bed, happy_bin=arguments.happy_bin, plot=arguments.plot, + chrom_aliases=arguments.chrom_aliases, ) diff --git a/neat/common/__init__.py b/neat/common/__init__.py index ab96bc3f..93cd432f 100644 --- a/neat/common/__init__.py +++ b/neat/common/__init__.py @@ -3,3 +3,4 @@ from .io import * from .constants_and_defaults import * from .ploid_functions import * +from .chrom_names import * diff --git a/neat/common/chrom_names.py b/neat/common/chrom_names.py new file mode 100644 index 00000000..46f1d8d2 --- /dev/null +++ b/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_namecanonical_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) diff --git a/neat/compare_vcfs/attribution.py b/neat/compare_vcfs/attribution.py index cb4d32f5..88ab3f86 100644 --- a/neat/compare_vcfs/attribution.py +++ b/neat/compare_vcfs/attribution.py @@ -17,9 +17,12 @@ import logging from pathlib import Path +from ..common.chrom_names import find_aliases + __all__ = [ "attribute_fn", "attribute_fns", + "detect_chrom_naming_mismatches", "load_bed_intervals", "position_in_intervals", "REASON_OUTSIDE_CONTIGS", @@ -36,13 +39,20 @@ REASON_UNKNOWN = "unknown" -def load_bed_intervals(bed_path: Path | str | None) -> dict[str, list[tuple[int, int]]] | None: +def load_bed_intervals( + bed_path: Path | str | None, + aliases: dict[str, str] | None = None, +) -> dict[str, list[tuple[int, int]]] | None: """ Parse a BED file into a per-contig list of sorted 0-based half-open intervals. Returns None if `bed_path` is None. Comments (`#`), `track`, and `browser` header lines are skipped. Rows with non-integer start/end are dropped with a debug log entry. + + If `aliases` is provided, each BED row's chrom field is remapped through + the dict (`aliases.get(chrom, chrom)`) before being stored, so downstream + lookups can use reference-canonical names. """ if bed_path is None: return None @@ -61,7 +71,10 @@ def load_bed_intervals(bed_path: Path | str | None) -> dict[str, list[tuple[int, except ValueError: _LOG.debug(f"{bed_path}:{lineno}: skipping non-integer interval") continue - intervals.setdefault(parts[0], []).append((start, end)) + chrom = parts[0] + if aliases: + chrom = aliases.get(chrom, chrom) + intervals.setdefault(chrom, []).append((start, end)) for chrom in intervals: intervals[chrom].sort() return intervals @@ -112,20 +125,86 @@ def attribute_fn( return reasons -def attribute_fns(fn_records, summary: dict) -> list[tuple]: +def attribute_fns( + fn_records, + summary: dict, + aliases: dict[str, str] | None = None, +) -> list[tuple]: """ Tag every FN against the run's simulation_summary. :param fn_records: iterable of pysam.VariantRecord (FN bucket from hap.py). :param summary: parsed simulation_summary.json. + :param aliases: optional user-supplied {bed_name: canonical_name} map applied + to BED chrom names at load time. :return: list of (record, reasons) tuples; `reasons` is a list[str]. """ contigs = frozenset(summary["delivered"].get("contigs_simulated", [])) cfg = summary.get("config", {}) - mutation_intervals = load_bed_intervals(cfg.get("mutation_bed")) - target_intervals = load_bed_intervals(cfg.get("target_bed")) + mutation_intervals = load_bed_intervals(cfg.get("mutation_bed"), aliases=aliases) + target_intervals = load_bed_intervals(cfg.get("target_bed"), aliases=aliases) return [ (rec, attribute_fn(rec.chrom, rec.pos, contigs, mutation_intervals, target_intervals)) for rec in fn_records ] + + +def detect_chrom_naming_mismatches( + summary: dict, + aliases: dict[str, str] | None = None, +) -> list[dict]: + """ + Inspect each configured BED's chrom set against the reference's; return a + list of warning records (one per mismatched BED) or empty if all BEDs + overlap (after any user-supplied aliases are applied). + + The "reference" set comes from `summary['delivered']['reference_contigs']` + (the full FASTA contig list), with a fallback to `contigs_simulated` for + backward compatibility with summaries written before that field existed. + """ + aliases = aliases or {} + reference_chroms = frozenset( + summary["delivered"].get("reference_contigs") + or summary["delivered"].get("contigs_simulated", []) + ) + if not reference_chroms: + return [] + + warnings: list[dict] = [] + for bed_label in ("mutation_bed", "target_bed"): + bed_path = summary.get("config", {}).get(bed_label) + if not bed_path: + continue + try: + intervals = load_bed_intervals(bed_path, aliases=None) + except FileNotFoundError: + continue # surface this separately; not a naming issue + if not intervals: + continue + + raw_chroms = frozenset(intervals.keys()) + mapped_chroms = frozenset(aliases.get(c, c) for c in raw_chroms) + if mapped_chroms & reference_chroms: + continue # at least one chrom matches; not a mismatch worth flagging + + suggested = find_aliases(raw_chroms, reference_chroms) + message = ( + f"{bed_label} chrom names don't overlap with reference contigs" + + ( + f"; suggested --chrom-aliases mappings: {suggested}" + if suggested + else " and no naming convention could be inferred. Attribution against " + "this BED will report every FN as 'outside'." + ) + ) + warnings.append({ + "type": "chrom_naming_mismatch", + "bed": bed_label, + "bed_path": str(bed_path), + "bed_chroms_sample": sorted(raw_chroms)[:5], + "reference_chroms_sample": sorted(reference_chroms)[:5], + "suggested_aliases": suggested, + "message": message, + }) + return warnings diff --git a/neat/compare_vcfs/reports.py b/neat/compare_vcfs/reports.py index 80ff4491..39980db6 100644 --- a/neat/compare_vcfs/reports.py +++ b/neat/compare_vcfs/reports.py @@ -73,6 +73,7 @@ def build_comparison_summary( fn_with_reasons_vcf: Path, comparison_summary_json: Path, comparison_summary_txt: Path, + warnings: list[dict] | None = None, ) -> dict: """Assemble the comparison_summary dict from the run's artifacts.""" return { @@ -92,6 +93,7 @@ def build_comparison_summary( "counts": dict(counts), "metrics": compute_metrics(counts), "fn_attribution": dict(fn_attribution), + "warnings": list(warnings) if warnings else [], "outputs": { "fn_with_reasons_vcf": str(Path(fn_with_reasons_vcf).resolve()), "comparison_summary_json": str(Path(comparison_summary_json).resolve()), @@ -162,6 +164,12 @@ def render_summary_txt(summary: dict) -> str: else: lines.append(" (no false negatives)") + warnings = summary.get("warnings") or [] + if warnings: + lines += ["", "Warnings", "--------"] + for w in warnings: + lines.append(f" - {w.get('message', w)}") + lines += [ "", "Outputs", diff --git a/neat/compare_vcfs/runner.py b/neat/compare_vcfs/runner.py index fdc1777c..bf103ca3 100644 --- a/neat/compare_vcfs/runner.py +++ b/neat/compare_vcfs/runner.py @@ -20,8 +20,9 @@ from pathlib import Path from ..common import validate_input_path +from ..common.chrom_names import load_chrom_aliases from ..read_simulator.utils.simulation_summary import SCHEMA_VERSION -from .attribution import attribute_fns +from .attribution import attribute_fns, detect_chrom_naming_mismatches from .happy import run_happy, parse_happy_output from .reports import ( build_comparison_summary, @@ -61,6 +62,7 @@ def compare_vcfs_runner( target_bed: str | None = None, happy_bin: str | None = None, plot: bool = False, + chrom_aliases: str | None = None, ): """ Run the comparison pipeline. @@ -89,6 +91,11 @@ def compare_vcfs_runner( validate_input_path(Path(reference).resolve()) if target_bed is not None: validate_input_path(Path(target_bed).resolve()) + if chrom_aliases is not None: + validate_input_path(Path(chrom_aliases).resolve()) + aliases = load_chrom_aliases(chrom_aliases) + if aliases: + _LOG.info(f"Loaded {len(aliases)} chrom alias mapping(s) from {chrom_aliases}") happy = discover_happy(happy_bin) _LOG.info(f"Using hap.py at: {happy}") @@ -119,7 +126,11 @@ def compare_vcfs_runner( f"FN={len(buckets['FN'])} FP={len(buckets['FP'])}" ) - fn_reasons = attribute_fns(buckets["FN"], summary) + chrom_warnings = detect_chrom_naming_mismatches(summary, aliases=aliases) + for w in chrom_warnings: + _LOG.warning(w["message"]) + + fn_reasons = attribute_fns(buckets["FN"], summary, aliases=aliases) reason_counts = summarize_fn_reasons(fn_reasons) if fn_reasons: _LOG.info(f"FN attribution: {reason_counts}") @@ -144,6 +155,7 @@ def compare_vcfs_runner( fn_with_reasons_vcf=fn_with_reasons_path, comparison_summary_json=comparison_json_path, comparison_summary_txt=comparison_txt_path, + warnings=chrom_warnings, ) write_comparison_summary_json(report, comparison_json_path) write_comparison_summary_txt(report, comparison_txt_path) diff --git a/neat/read_simulator/runner.py b/neat/read_simulator/runner.py index 6b3d541d..cde187a1 100644 --- a/neat/read_simulator/runner.py +++ b/neat/read_simulator/runner.py @@ -334,6 +334,7 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): config_path=config, analysis_start=analysis_start, contigs_simulated=list(input_variants_dict.keys()), + reference_contigs=list(reference_keys_with_lens.keys()), ) def filter_thread_variants(contig_variants: ContigVariants, coords: tuple[int, int]) -> ContigVariants: diff --git a/neat/read_simulator/utils/simulation_summary.py b/neat/read_simulator/utils/simulation_summary.py index e8db4d06..5f8f3bad 100644 --- a/neat/read_simulator/utils/simulation_summary.py +++ b/neat/read_simulator/utils/simulation_summary.py @@ -29,6 +29,7 @@ def write_simulation_summary( config_path: Path, analysis_start: float, contigs_simulated: Iterable[str], + reference_contigs: Iterable[str] | None = None, ) -> Path: """ Write `simulation_summary.json` into `output_dir`. @@ -90,6 +91,7 @@ def write_simulation_summary( "total_variants": total_variants, "variants_by_contig": variants_by_contig, "contigs_simulated": list(contigs_simulated), + "reference_contigs": list(reference_contigs) if reference_contigs is not None else list(contigs_simulated), }, } diff --git a/tests/test_common/__init__.py b/tests/test_common/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_common/test_chrom_names.py b/tests/test_common/test_chrom_names.py new file mode 100644 index 00000000..90a948a9 --- /dev/null +++ b/tests/test_common/test_chrom_names.py @@ -0,0 +1,144 @@ +""" +Tests for neat/common/chrom_names.py — chromosome-name normalization helpers. +""" +from pathlib import Path + +from neat.common.chrom_names import ( + apply_aliases, + find_aliases, + load_chrom_aliases, + prefix_flip_candidates, +) + + +# =========================================================================== +# prefix_flip_candidates +# =========================================================================== + +def test_prefix_flip_strips_chr_prefix(): + assert prefix_flip_candidates("chr1") == {"1"} + + +def test_prefix_flip_adds_chr_prefix(): + assert prefix_flip_candidates("1") == {"chr1"} + + +def test_prefix_flip_handles_sex_chroms(): + assert prefix_flip_candidates("chrX") == {"X"} + assert prefix_flip_candidates("Y") == {"chrY"} + + +def test_prefix_flip_mitochondrial_chrM_yields_all_three_alternatives(): + """`chrM` should yield {M, MT, chrMT} — all known mt-name conventions.""" + assert prefix_flip_candidates("chrM") == {"M", "MT", "chrMT"} + + +def test_prefix_flip_mitochondrial_MT_yields_all_three_alternatives(): + assert prefix_flip_candidates("MT") == {"M", "chrM", "chrMT"} + + +def test_prefix_flip_excludes_self(): + """The input name is never present in its own candidate set.""" + for name in ("chr1", "1", "chrM", "MT", "scaffold_42"): + assert name not in prefix_flip_candidates(name) + + +def test_prefix_flip_custom_scaffold_name_only_does_prefix(): + """Non-standard names get the chr-prefix heuristic but no mitochondrial mapping.""" + assert prefix_flip_candidates("scaffold_42") == {"chrscaffold_42"} + + +# =========================================================================== +# find_aliases +# =========================================================================== + +def test_find_aliases_native_match_omitted(): + """Names that already appear in target are skipped.""" + assert find_aliases({"chr1"}, {"chr1", "chr2"}) == {} + + +def test_find_aliases_prefix_flip_proposed(): + assert find_aliases({"1", "2"}, {"chr1", "chr2"}) == {"1": "chr1", "2": "chr2"} + + +def test_find_aliases_mitochondrial_proposed(): + assert find_aliases({"MT"}, {"chrM"}) == {"MT": "chrM"} + + +def test_find_aliases_unmatchable_omitted(): + """A name with no candidate in target is not in the result — caller knows the + BED is genuinely non-overlapping rather than just misnamed.""" + assert find_aliases({"weird_contig"}, {"chr1"}) == {} + + +def test_find_aliases_mixed_some_aliased_some_native(): + """A BED can have both naming conventions; aliases only cover mismatches.""" + result = find_aliases({"chr1", "2", "MT"}, {"chr1", "chr2", "chrM"}) + assert result == {"2": "chr2", "MT": "chrM"} + + +def test_find_aliases_empty_inputs(): + assert find_aliases(set(), {"chr1"}) == {} + assert find_aliases({"chr1"}, set()) == {} + + +# =========================================================================== +# load_chrom_aliases +# =========================================================================== + +def test_load_chrom_aliases_returns_empty_for_none(): + assert load_chrom_aliases(None) == {} + + +def test_load_chrom_aliases_parses_tab_separated(tmp_path): + f = tmp_path / "aliases.tsv" + f.write_text("1\tchr1\n2\tchr2\nMT\tchrM\n") + assert load_chrom_aliases(f) == {"1": "chr1", "2": "chr2", "MT": "chrM"} + + +def test_load_chrom_aliases_skips_comments_and_blanks(tmp_path): + f = tmp_path / "aliases.tsv" + f.write_text( + "# This is a comment\n" + "\n" + "1\tchr1\n" + "# another comment\n" + "2\tchr2\n" + ) + assert load_chrom_aliases(f) == {"1": "chr1", "2": "chr2"} + + +def test_load_chrom_aliases_tolerates_whitespace_separator(tmp_path): + """Human-typed files often use spaces instead of literal tabs.""" + f = tmp_path / "aliases.tsv" + f.write_text("1 chr1\n2 chr2\n") + assert load_chrom_aliases(f) == {"1": "chr1", "2": "chr2"} + + +def test_load_chrom_aliases_skips_malformed_lines(tmp_path): + f = tmp_path / "aliases.tsv" + f.write_text("just_one_token\n1\tchr1\n") + assert load_chrom_aliases(f) == {"1": "chr1"} + + +def test_load_chrom_aliases_later_entries_override_earlier(tmp_path): + """If a source name is listed twice, the last entry wins (documents behavior).""" + f = tmp_path / "aliases.tsv" + f.write_text("1\tchr1\n1\tchrONE\n") + assert load_chrom_aliases(f) == {"1": "chrONE"} + + +# =========================================================================== +# apply_aliases +# =========================================================================== + +def test_apply_aliases_returns_canonical_when_mapped(): + assert apply_aliases("1", {"1": "chr1"}) == "chr1" + + +def test_apply_aliases_returns_input_when_unmapped(): + assert apply_aliases("chr1", {"1": "chr1"}) == "chr1" + + +def test_apply_aliases_empty_map_passthrough(): + assert apply_aliases("chr1", {}) == "chr1" diff --git a/tests/test_compare_vcfs/test_attribution.py b/tests/test_compare_vcfs/test_attribution.py index 589102ad..bc70e4f6 100644 --- a/tests/test_compare_vcfs/test_attribution.py +++ b/tests/test_compare_vcfs/test_attribution.py @@ -13,6 +13,7 @@ REASON_UNKNOWN, attribute_fn, attribute_fns, + detect_chrom_naming_mismatches, load_bed_intervals, position_in_intervals, ) @@ -256,3 +257,114 @@ def test_attribute_fns_pairs_each_record_with_reasons(tmp_path): [(returned_rec, reasons)] = attribute_fns([rec], summary) assert returned_rec is rec assert reasons == [REASON_UNKNOWN] + + +def test_attribute_fns_applies_aliases_to_bed_chrom_keys(tmp_path): + """User aliases must rewrite BED chrom names so the FN check finds them.""" + bed = tmp_path / "mut.bed" + bed.write_text("1\t0\t1000\n") # BED uses '1', reference uses 'chr1' + summary = { + "delivered": {"contigs_simulated": ["chr1"], "reference_contigs": ["chr1"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + fns = [_fake_record("chr1", 500)] + # Without aliases — FN is outside (BED chrom doesn't match ref) + [(_, no_alias_reasons)] = attribute_fns(fns, summary) + assert no_alias_reasons == [REASON_OUTSIDE_MUTATION_BED] + # With aliases — FN is inside (BED's '1' is normalized to 'chr1') + [(_, aliased_reasons)] = attribute_fns(fns, summary, aliases={"1": "chr1"}) + assert aliased_reasons == [REASON_UNKNOWN] + + +# =========================================================================== +# detect_chrom_naming_mismatches +# =========================================================================== + +def test_detect_chrom_mismatch_empty_when_no_beds_configured(): + summary = { + "delivered": {"reference_contigs": ["chr1"], "contigs_simulated": ["chr1"]}, + "config": {"mutation_bed": None, "target_bed": None}, + } + assert detect_chrom_naming_mismatches(summary) == [] + + +def test_detect_chrom_mismatch_empty_when_overlap_exists(tmp_path): + """If even one chrom overlaps, no warning — partial match is acceptable.""" + bed = tmp_path / "mut.bed" + bed.write_text("chr1\t0\t1000\nweird\t0\t100\n") + summary = { + "delivered": {"reference_contigs": ["chr1", "chr2"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + assert detect_chrom_naming_mismatches(summary) == [] + + +def test_detect_chrom_mismatch_suggests_prefix_aliases(tmp_path): + """When the BED uses '1'/'2' and reference uses 'chr1'/'chr2', the warning + must include a prefix-flip mapping in suggested_aliases.""" + bed = tmp_path / "mut.bed" + bed.write_text("1\t0\t1000\n2\t0\t1000\n") + summary = { + "delivered": {"reference_contigs": ["chr1", "chr2"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + warnings = detect_chrom_naming_mismatches(summary) + assert len(warnings) == 1 + w = warnings[0] + assert w["type"] == "chrom_naming_mismatch" + assert w["bed"] == "mutation_bed" + assert w["suggested_aliases"] == {"1": "chr1", "2": "chr2"} + assert "--chrom-aliases" in w["message"] + + +def test_detect_chrom_mismatch_no_suggestion_when_inscrutable(tmp_path): + """Names with no prefix-flip and no mt match → warning with empty suggested_aliases.""" + bed = tmp_path / "mut.bed" + bed.write_text("weird_contig\t0\t1000\n") + summary = { + "delivered": {"reference_contigs": ["chr1"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + warnings = detect_chrom_naming_mismatches(summary) + assert len(warnings) == 1 + assert warnings[0]["suggested_aliases"] == {} + assert "no naming convention" in warnings[0]["message"] + + +def test_detect_chrom_mismatch_silenced_by_user_aliases(tmp_path): + """If user supplies aliases that resolve the mismatch, no warning is emitted.""" + bed = tmp_path / "mut.bed" + bed.write_text("1\t0\t1000\n") + summary = { + "delivered": {"reference_contigs": ["chr1"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + assert detect_chrom_naming_mismatches(summary, aliases={"1": "chr1"}) == [] + + +def test_detect_chrom_mismatch_falls_back_to_contigs_simulated(tmp_path): + """When reference_contigs is absent, fall back to contigs_simulated (back-compat).""" + bed = tmp_path / "mut.bed" + bed.write_text("1\t0\t1000\n") + summary = { + "delivered": {"contigs_simulated": ["chr1"]}, # no reference_contigs + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + warnings = detect_chrom_naming_mismatches(summary) + assert len(warnings) == 1 + assert warnings[0]["suggested_aliases"] == {"1": "chr1"} + + +def test_detect_chrom_mismatch_reports_one_warning_per_bed(tmp_path): + """Each mismatched BED gets its own warning entry.""" + mut = tmp_path / "mut.bed" + mut.write_text("1\t0\t1000\n") + tgt = tmp_path / "tgt.bed" + tgt.write_text("2\t0\t1000\n") + summary = { + "delivered": {"reference_contigs": ["chr1", "chr2"]}, + "config": {"mutation_bed": str(mut), "target_bed": str(tgt)}, + } + warnings = detect_chrom_naming_mismatches(summary) + assert len(warnings) == 2 + assert {w["bed"] for w in warnings} == {"mutation_bed", "target_bed"} diff --git a/tests/test_compare_vcfs/test_cli.py b/tests/test_compare_vcfs/test_cli.py index 341d4baf..31f97bfd 100644 --- a/tests/test_compare_vcfs/test_cli.py +++ b/tests/test_compare_vcfs/test_cli.py @@ -39,6 +39,7 @@ def test_compare_vcfs_help_text_lists_all_documented_flags(): "--target-bed", "--happy-bin", "--plot", + "--chrom-aliases", ]: assert flag in help_text, f"--help missing {flag}" @@ -67,6 +68,7 @@ def fake_runner(**kwargs): "--target-bed", "/t.bed", "--happy-bin", "/bin/hap.py", "--plot", + "--chrom-aliases", "/aliases.tsv", ]) cmd.execute(args) @@ -79,6 +81,7 @@ def fake_runner(**kwargs): "target_bed": "/t.bed", "happy_bin": "/bin/hap.py", "plot": True, + "chrom_aliases": "/aliases.tsv", } @@ -95,3 +98,4 @@ def test_command_optional_flags_default_to_none_or_false(): assert args.target_bed is None assert args.happy_bin is None assert args.plot is False + assert args.chrom_aliases is None diff --git a/tests/test_compare_vcfs/test_reports.py b/tests/test_compare_vcfs/test_reports.py index f466452b..062a84dd 100644 --- a/tests/test_compare_vcfs/test_reports.py +++ b/tests/test_compare_vcfs/test_reports.py @@ -118,7 +118,8 @@ def test_build_comparison_summary_has_required_top_level_keys(tmp_path): s = _build_minimal_summary(tmp_path) assert set(s) == { "schema_version", "neat_version", "generated_at", - "inputs", "happy", "counts", "metrics", "fn_attribution", "outputs", + "inputs", "happy", "counts", "metrics", "fn_attribution", + "warnings", "outputs", } diff --git a/tests/test_compare_vcfs/test_runner.py b/tests/test_compare_vcfs/test_runner.py index 823ee924..979e5f3d 100644 --- a/tests/test_compare_vcfs/test_runner.py +++ b/tests/test_compare_vcfs/test_runner.py @@ -367,3 +367,65 @@ def test_runner_validates_optional_target_bed_when_provided(tmp_path): target_bed=str(tmp_path / "missing_target.bed"), ) assert excinfo.value.code == 5 + + +# =========================================================================== +# Chrom-naming mismatch end-to-end through the runner +# =========================================================================== + +def test_runner_emits_chrom_mismatch_warning_into_json_report(tmp_path, monkeypatch): + """A BED with mismatched chrom names must surface as a warning in the JSON.""" + mut = tmp_path / "mut.bed" + mut.write_text("1\t0\t1000\n") + run_dir = _make_summary_dir(tmp_path, extra={ + "delivered": {"reference_contigs": ["chr1"], "contigs_simulated": ["chr1"]}, + "config": {"mutation_bed": str(mut), "target_bed": None}, + }) + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + _install_fake_happy(monkeypatch, _HAPPY_HEADER_ONE_FN) + + out_dir = tmp_path / "out" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(out_dir), + ) + report = json.loads((out_dir / "comparison_summary.json").read_text()) + assert len(report["warnings"]) == 1 + w = report["warnings"][0] + assert w["type"] == "chrom_naming_mismatch" + assert w["suggested_aliases"] == {"1": "chr1"} + + +def test_runner_chrom_aliases_silences_warning_and_fixes_attribution(tmp_path, monkeypatch): + """Passing --chrom-aliases must both remove the warning AND make the FN + correctly attributed (no longer 'outside_mutation_bed').""" + mut = tmp_path / "mut.bed" + mut.write_text("1\t0\t1000\n") # FN at chr1:200 falls inside if alias applied + run_dir = _make_summary_dir(tmp_path, extra={ + "delivered": {"reference_contigs": ["chr1"], "contigs_simulated": ["chr1"]}, + "config": {"mutation_bed": str(mut), "target_bed": None}, + }) + aliases = tmp_path / "aliases.tsv" + aliases.write_text("1\tchr1\n") + + golden = _touch(tmp_path / "golden.vcf") + called = _touch(tmp_path / "called.vcf") + _make_executable(tmp_path / "hap.py") + monkeypatch.setenv("PATH", str(tmp_path), prepend=os.pathsep) + _install_fake_happy(monkeypatch, _HAPPY_HEADER_ONE_FN) + + out_dir = tmp_path / "out" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(run_dir), + output_dir=str(out_dir), + chrom_aliases=str(aliases), + ) + report = json.loads((out_dir / "comparison_summary.json").read_text()) + assert report["warnings"] == [] + # FN at chr1:200, BED's '1' mapped to 'chr1' → position 200 IS inside [0, 1000) + assert report["fn_attribution"] == {"unknown": 1} diff --git a/tests/test_read_simulator/test_simulation_summary.py b/tests/test_read_simulator/test_simulation_summary.py index 38cc216f..59d84650 100644 --- a/tests/test_read_simulator/test_simulation_summary.py +++ b/tests/test_read_simulator/test_simulation_summary.py @@ -297,6 +297,32 @@ def test_write_simulation_summary_contigs_simulated_preserves_order(tmp_path): assert summary["delivered"]["contigs_simulated"] == ["chr3", "chr1", "chr2"] +def test_write_simulation_summary_records_reference_contigs(tmp_path): + """`reference_contigs` captures the full FASTA contig set, separate from + `contigs_simulated` (which is what the simulator iterated over).""" + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, + contigs_simulated=["chr1"], + reference_contigs=["chr1", "chr2", "chrX", "chrM"], + ) + delivered = json.loads((tmp_path / "simulation_summary.json").read_text())["delivered"] + assert delivered["reference_contigs"] == ["chr1", "chr2", "chrX", "chrM"] + assert delivered["contigs_simulated"] == ["chr1"] + + +def test_write_simulation_summary_reference_contigs_defaults_to_contigs_simulated(tmp_path): + """When reference_contigs is omitted, fall back to contigs_simulated so old + callers that don't pass it still get a populated field.""" + write_simulation_summary( + options=_make_options(), output_dir=tmp_path, file_prefix="run", + config_path=tmp_path / "cfg.yml", analysis_start=0.0, + contigs_simulated=["chr1", "chr2"], + ) + delivered = json.loads((tmp_path / "simulation_summary.json").read_text())["delivered"] + assert delivered["reference_contigs"] == ["chr1", "chr2"] + + def test_write_simulation_summary_outputs_none_when_nothing_produced(tmp_path): write_simulation_summary( options=_make_options(), output_dir=tmp_path, file_prefix="run", From bf5ab16a92efa46e4b41fe6dbdff17a44037dcf1 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 22:44:30 -0500 Subject: [PATCH 4/4] Close outside_* semantic gap + fill pre-merge test gaps Semantic fix: when chrom-name mismatches make a BED unusable, the runner now skips attribution against that BED entirely (via new skip_beds param on attribute_fns) instead of tagging every FN as 'outside_'. The warning in comparison_summary.json explains the cause; the fn_attribution counts now reflect what NEAT actually checked. Test additions (7): - render_summary_txt with/without warnings section (rendering was untested) - detect_chrom_naming_mismatches skips empty BEDs (no spurious warning) - detect_chrom_naming_mismatches still warns when user aliases don't resolve the mismatch (suggested mapping comes from raw BED chroms, not user input) - attribute_fns skip_beds suppresses outside_ tags - attribute_fns skip_beds doesn't affect outside_simulated_contigs - Real hap.py end-to-end with mismatched mutation_bed (warning fires, outside_mutation_bed absent, --chrom-aliases silences warning) Existing test updated: - test_runner_emits_chrom_mismatch_warning_into_json_report now also asserts the semantic guarantee that outside_mutation_bed is absent. Full suite: 768 passing. Co-Authored-By: Claude Opus 4.7 (1M context) --- neat/compare_vcfs/attribution.py | 16 ++++- neat/compare_vcfs/runner.py | 7 +- tests/test_compare_vcfs/test_attribution.py | 61 ++++++++++++++++ tests/test_compare_vcfs/test_integration.py | 77 +++++++++++++++++++++ tests/test_compare_vcfs/test_reports.py | 20 ++++++ tests/test_compare_vcfs/test_runner.py | 4 ++ 6 files changed, 182 insertions(+), 3 deletions(-) diff --git a/neat/compare_vcfs/attribution.py b/neat/compare_vcfs/attribution.py index 88ab3f86..dc641451 100644 --- a/neat/compare_vcfs/attribution.py +++ b/neat/compare_vcfs/attribution.py @@ -129,6 +129,7 @@ def attribute_fns( fn_records, summary: dict, aliases: dict[str, str] | None = None, + skip_beds: set[str] | None = None, ) -> list[tuple]: """ Tag every FN against the run's simulation_summary. @@ -137,12 +138,23 @@ def attribute_fns( :param summary: parsed simulation_summary.json. :param aliases: optional user-supplied {bed_name: canonical_name} map applied to BED chrom names at load time. + :param skip_beds: optional set of BED labels (e.g., {"mutation_bed"}) to + treat as if not configured. Used by the runner to skip BEDs whose chrom + names are entirely mismatched against the reference — attributing FNs + against an unusable BED would produce misleading `outside_*` counts. :return: list of (record, reasons) tuples; `reasons` is a list[str]. """ + skip_beds = skip_beds or set() contigs = frozenset(summary["delivered"].get("contigs_simulated", [])) cfg = summary.get("config", {}) - mutation_intervals = load_bed_intervals(cfg.get("mutation_bed"), aliases=aliases) - target_intervals = load_bed_intervals(cfg.get("target_bed"), aliases=aliases) + mutation_intervals = ( + None if "mutation_bed" in skip_beds + else load_bed_intervals(cfg.get("mutation_bed"), aliases=aliases) + ) + target_intervals = ( + None if "target_bed" in skip_beds + else load_bed_intervals(cfg.get("target_bed"), aliases=aliases) + ) return [ (rec, attribute_fn(rec.chrom, rec.pos, contigs, mutation_intervals, target_intervals)) diff --git a/neat/compare_vcfs/runner.py b/neat/compare_vcfs/runner.py index bf103ca3..d75c3395 100644 --- a/neat/compare_vcfs/runner.py +++ b/neat/compare_vcfs/runner.py @@ -130,7 +130,12 @@ def compare_vcfs_runner( for w in chrom_warnings: _LOG.warning(w["message"]) - fn_reasons = attribute_fns(buckets["FN"], summary, aliases=aliases) + # A fully mismatched BED is unusable for attribution; skip it so we don't + # report misleading 'outside_' counts for chroms NEAT never actually + # checked against the BED. The warning surfaces the underlying cause. + unusable_beds = {w["bed"] for w in chrom_warnings if w.get("type") == "chrom_naming_mismatch"} + + fn_reasons = attribute_fns(buckets["FN"], summary, aliases=aliases, skip_beds=unusable_beds) reason_counts = summarize_fn_reasons(fn_reasons) if fn_reasons: _LOG.info(f"FN attribution: {reason_counts}") diff --git a/tests/test_compare_vcfs/test_attribution.py b/tests/test_compare_vcfs/test_attribution.py index bc70e4f6..6bff8f61 100644 --- a/tests/test_compare_vcfs/test_attribution.py +++ b/tests/test_compare_vcfs/test_attribution.py @@ -355,6 +355,67 @@ def test_detect_chrom_mismatch_falls_back_to_contigs_simulated(tmp_path): assert warnings[0]["suggested_aliases"] == {"1": "chr1"} +def test_detect_chrom_mismatch_skips_empty_bed(tmp_path): + """A BED with only comments/blanks has no chroms — no warning should fire.""" + bed = tmp_path / "empty.bed" + bed.write_text("# only a comment\n\n") + summary = { + "delivered": {"reference_contigs": ["chr1"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + assert detect_chrom_naming_mismatches(summary) == [] + + +def test_detect_chrom_mismatch_warns_when_user_aliases_dont_resolve(tmp_path): + """Wrong user aliases (mapping to a name not in the reference) should NOT + silence the warning — and the suggested mapping must come from the raw BED + chroms, not the user's incorrect alias output.""" + bed = tmp_path / "mut.bed" + bed.write_text("1\t0\t1000\n") + summary = { + "delivered": {"reference_contigs": ["chr1"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + warnings = detect_chrom_naming_mismatches(summary, aliases={"1": "definitely_not_a_chrom"}) + assert len(warnings) == 1 + assert warnings[0]["suggested_aliases"] == {"1": "chr1"} + + +# =========================================================================== +# attribute_fns with skip_beds (semantic gap: unusable BED → no outside_* tag) +# =========================================================================== + +def test_attribute_fns_skip_beds_suppresses_outside_mutation_bed(tmp_path): + """When the runner marks mutation_bed as unusable, attribute_fns must not + report 'outside_mutation_bed' for any FN — even when chrom is in the BED + on paper. Locks in the semantic fix for the chrom-mismatch case.""" + bed = tmp_path / "mut.bed" + bed.write_text("chrZ\t0\t1000\n") # BED has chroms unrelated to FN locations + summary = { + "delivered": {"contigs_simulated": ["chr1"], "reference_contigs": ["chr1"]}, + "config": {"mutation_bed": str(bed), "target_bed": None}, + } + fns = [_fake_record("chr1", 500)] + # Without skip_beds: outside_mutation_bed (chr1 not in BED keys) + [(_, default_reasons)] = attribute_fns(fns, summary) + assert default_reasons == [REASON_OUTSIDE_MUTATION_BED] + # With skip_beds: BED is ignored entirely → no other reason applies → unknown + [(_, skipped_reasons)] = attribute_fns(fns, summary, skip_beds={"mutation_bed"}) + assert skipped_reasons == [REASON_UNKNOWN] + + +def test_attribute_fns_skip_beds_does_not_affect_outside_simulated_contigs(): + """skip_beds is about BED checks; the simulated-contigs check still fires + for chroms NEAT never ran on.""" + summary = { + "delivered": {"contigs_simulated": ["chr1"], "reference_contigs": ["chr1"]}, + "config": {"mutation_bed": None, "target_bed": None}, + } + fns = [_fake_record("chrZ", 100)] + [(_, reasons)] = attribute_fns(fns, summary, skip_beds={"mutation_bed", "target_bed"}) + assert reasons == [REASON_OUTSIDE_CONTIGS] + + def test_detect_chrom_mismatch_reports_one_warning_per_bed(tmp_path): """Each mismatched BED gets its own warning entry.""" mut = tmp_path / "mut.bed" diff --git a/tests/test_compare_vcfs/test_integration.py b/tests/test_compare_vcfs/test_integration.py index 08c1f555..cc903a16 100644 --- a/tests/test_compare_vcfs/test_integration.py +++ b/tests/test_compare_vcfs/test_integration.py @@ -217,3 +217,80 @@ def test_compare_vcfs_real_happy_multi_fn(tmp_path, happy_bin, happy_env_path, m # --plot was on, so the bar chart should be present assert (cmp_out / "fn_attribution.png").is_file() + + +def test_compare_vcfs_real_happy_with_chrom_mismatched_bed(tmp_path, happy_bin, happy_env_path, monkeypatch): + """ + End-to-end against real hap.py with a mutation_bed that uses '1'/'2' while + the reference uses 'chr1'/'chr2'. Verifies that: + - the warning surfaces in the JSON report + - FNs are NOT mislabeled as 'outside_mutation_bed' (the semantic fix) + - --chrom-aliases silences the warning AND restores correct attribution + """ + monkeypatch.setenv("PATH", happy_env_path + os.pathsep + os.environ["PATH"]) + + # Hand-craft a mutation_bed that won't match the reference's chr-prefixed names. + # The simulator will log warnings and effectively ignore it (NEAT's pre-existing + # "skip BED chroms not in reference" behavior). simulation_summary still records + # the BED path; compare-vcfs reads it and detects the mismatch. + bed = tmp_path / "mismatched_mut.bed" + bed.write_text("1\t0\t1000\n2\t0\t1000\n") + + sim_out = tmp_path / "sim_out" + sim_out.mkdir() + ref = _write_ref(tmp_path / "ref.fa") + cfg_path = tmp_path / "conf.yml" + cfg_path.write_text( + f"reference: {ref}\n" + "produce_fastq: false\n" + "produce_bam: false\n" + "produce_vcf: true\n" + "read_len: 100\n" + "coverage: 5\n" + "rng_seed: 42\n" + "mutation_rate: 0.01\n" + f"mutation_bed: {bed}\n" + "overwrite_output: true\n" + "cleanup_splits: true\n", + encoding="utf-8", + ) + read_simulator_runner(str(cfg_path), str(sim_out), "run") + golden = sim_out / "run_golden.vcf.gz" + + called = tmp_path / "called.vcf.gz" + _called_vcf_dropping_first_variant(golden, called) + + # ------------------------------------------------------------------ + # Run 1: no --chrom-aliases → mismatch warning, no misleading outside_* + # ------------------------------------------------------------------ + cmp_out_no_aliases = tmp_path / "cmp_no_aliases" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(sim_out), + output_dir=str(cmp_out_no_aliases), + reference=str(ref), + happy_bin=str(happy_bin), + ) + report = json.loads((cmp_out_no_aliases / "comparison_summary.json").read_text()) + chrom_warnings = [w for w in report["warnings"] if w["type"] == "chrom_naming_mismatch"] + assert chrom_warnings, "expected chrom-mismatch warning to fire" + assert chrom_warnings[0]["suggested_aliases"] == {"1": "chr1", "2": "chr2"} + # Semantic guarantee: mismatched BED must NOT produce outside_mutation_bed + assert "outside_mutation_bed" not in report["fn_attribution"] + + # ------------------------------------------------------------------ + # Run 2: --chrom-aliases supplied → no warning, BED becomes usable + # ------------------------------------------------------------------ + aliases = tmp_path / "aliases.tsv" + aliases.write_text("1\tchr1\n2\tchr2\n") + cmp_out_aliased = tmp_path / "cmp_aliased" + compare_vcfs_runner( + golden_vcf=str(golden), called_vcf=str(called), + neat_run_dir=str(sim_out), + output_dir=str(cmp_out_aliased), + reference=str(ref), + happy_bin=str(happy_bin), + chrom_aliases=str(aliases), + ) + report_aliased = json.loads((cmp_out_aliased / "comparison_summary.json").read_text()) + assert not [w for w in report_aliased["warnings"] if w["type"] == "chrom_naming_mismatch"] diff --git a/tests/test_compare_vcfs/test_reports.py b/tests/test_compare_vcfs/test_reports.py index 062a84dd..4bf4828c 100644 --- a/tests/test_compare_vcfs/test_reports.py +++ b/tests/test_compare_vcfs/test_reports.py @@ -207,6 +207,26 @@ def test_render_summary_txt_renders_NA_for_undefined_metrics(tmp_path): assert "N/A" in txt +def test_render_summary_txt_includes_warnings_section(tmp_path): + """When warnings are present, the txt report must surface them in a dedicated section.""" + s = _build_minimal_summary(tmp_path) + s["warnings"] = [ + {"type": "chrom_naming_mismatch", "message": "mutation_bed names don't overlap reference"}, + {"type": "chrom_naming_mismatch", "message": "target_bed names don't overlap reference"}, + ] + txt = render_summary_txt(s) + assert "Warnings" in txt + assert "mutation_bed names don't overlap reference" in txt + assert "target_bed names don't overlap reference" in txt + + +def test_render_summary_txt_omits_warnings_section_when_empty(tmp_path): + s = _build_minimal_summary(tmp_path) + s["warnings"] = [] + txt = render_summary_txt(s) + assert "\nWarnings\n" not in txt + + # =========================================================================== # write_fn_with_reasons # =========================================================================== diff --git a/tests/test_compare_vcfs/test_runner.py b/tests/test_compare_vcfs/test_runner.py index 979e5f3d..f1203f0f 100644 --- a/tests/test_compare_vcfs/test_runner.py +++ b/tests/test_compare_vcfs/test_runner.py @@ -398,6 +398,10 @@ def test_runner_emits_chrom_mismatch_warning_into_json_report(tmp_path, monkeypa w = report["warnings"][0] assert w["type"] == "chrom_naming_mismatch" assert w["suggested_aliases"] == {"1": "chr1"} + # Semantic guarantee: the FN must NOT be tagged 'outside_mutation_bed' + # because the BED was unusable. The warning surfaces the cause instead. + assert "outside_mutation_bed" not in report["fn_attribution"] + assert report["fn_attribution"] == {"unknown": 1} def test_runner_chrom_aliases_silences_warning_and_fixes_attribution(tmp_path, monkeypatch):