From 46dbd49af358f98505b71cf3975ed73d25b4086f Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 18:04:21 -0500 Subject: [PATCH 01/12] 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 02/12] 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 03/12] 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 04/12] 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): From 4a80456383bdc6bda7c8f9ac728dbd4b7658b04b Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 23:13:18 -0500 Subject: [PATCH 05/12] Drop cleanup_splits and reuse_splits; bump to v4.5.1 reuse_splits was broken (raised unconditionally regardless of whether the splits dir existed). cleanup_splits existed only to support it. Splits always go in the TemporaryDirectory that auto-cleans on exit. Both keys added to DEPRECATED_KEYS so old configs warn instead of silently no-oping. README examples, test fixtures, and test_options updated. Co-Authored-By: Claude Sonnet 4.6 --- ChangeLog.md | 15 ++++++ README.md | 4 -- neat/read_simulator/utils/options.py | 30 ++--------- pyproject.toml | 2 +- tests/test_compare_vcfs/test_integration.py | 4 +- .../test_gc_bias_pipeline.py | 1 - .../test_models_pipeline.py | 1 - tests/test_read_simulator/test_options.py | 53 ++++++++----------- tests/test_read_simulator/test_runner.py | 1 - 9 files changed, 43 insertions(+), 68 deletions(-) diff --git a/ChangeLog.md b/ChangeLog.md index 5cc4f683..bf730060 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,18 @@ +# NEAT v4.5.1 + +Removed the `cleanup_splits` and `reuse_splits` config options. + +`reuse_splits` was fully broken (it raised `FileNotFoundError` unconditionally +regardless of whether the splits directory existed). `cleanup_splits` existed +solely to support `reuse_splits`; without it the option has no value — the +splits directory always lives in a `TemporaryDirectory` that is cleaned up +automatically when the run exits. + +Both keys are now in `DEPRECATED_KEYS`: configs that still include them +receive a one-line deprecation warning and continue parsing cleanly, so no +user config breaks. The README examples and all internal test fixtures have +been updated to drop these keys. + # NEAT v4.5.0 New `neat compare-vcfs` subcommand: compares a downstream variant caller's VCF diff --git a/README.md b/README.md index f67d670d..73cb783f 100755 --- a/README.md +++ b/README.md @@ -205,8 +205,6 @@ More parameters are below: | `min_mutations` | Set the minimum number of mutations that NEAT should add, per contig. Default is 0. We recommend setting this to at least one for small chromosomes, so NEAT will produce at least one mutation per contig. | | `threads` | Number of threads to use. More than 1 will use multi-threading to speed up processing. With `threads > 1`, NEAT splits each contig into chunks; with `threads == 1`, one chunk per contig is used. | | `parallel_block_size` | Per-chunk size in bases when `threads > 1`. Default `0` (auto-tune from total genome length and thread count, targeting ~8 chunks per thread). Set to a positive integer to override. Ignored when `threads == 1`. | -| `cleanup_splits` | If running more than one simulation on the same input fasta, you can reuse splits files. By default, this will be set to `False`, and splits files will be deleted at the end of the run. | -| `reuse_splits` | If an existing splits file exists in the output folder, it will use those splits, if this value is set to `True`. | The command line options for NEAT are as follows: @@ -280,7 +278,6 @@ The inputs in single-threaded, contig-based mode most closely replicate the beha The configuration used: - `threads: 1` (NEAT processes one contig per chunk in single-thread mode) -- `cleanup_splits: True` | Organism | File size (bytes) | Avg. runtime (ms) | Avg. runtime (min) | |-----------------|-------------------|-------------------|--------------------| @@ -302,7 +299,6 @@ Here we enabled NEAT’s multi-threaded mode, which splits contigs into size-bas - `parallel_block_size: 500000` - `produce_bam: false` - `threads: 7` -- `cleanup_splits: True` | Organism | File size (bytes) | Avg. runtime (ms) | Avg. runtime (min) | |-----------------|-------------------|-------------------|--------------------| diff --git a/neat/read_simulator/utils/options.py b/neat/read_simulator/utils/options.py index 991857e7..193c3c4b 100644 --- a/neat/read_simulator/utils/options.py +++ b/neat/read_simulator/utils/options.py @@ -82,9 +82,7 @@ def __init__(self, produce_fastq: bool = True, min_mutations: int = 0, parallel_block_size: int = 0, - cleanup_splits: bool = True, splits_dir: Path | None = None, - reuse_splits: bool = False, gc_model: Path | None = None, **kwargs: Any ): @@ -136,9 +134,7 @@ def __init__(self, chunks of this size). Default 0 auto-tunes from total genome length and thread count, targeting ~8 chunks per thread. Specify a positive integer to override. Ignored when threads == 1, where one chunk per contig is used. - :param cleanup_splits: Set to False in order to preserve splits after run - :param reuse_splits: Attempts to reuse existing splits file - """ +""" super().__init__(**kwargs) self.reference: Path = reference self.output_dir: Path = output_dir @@ -173,9 +169,7 @@ def __init__(self, # threads > 1, contigs are split into chunks of `parallel_block_size`; with # threads == 1, each contig is processed as a single chunk. self.parallel_block_size: int = parallel_block_size - self.cleanup_splits: bool = cleanup_splits self.splits_dir: Path | None = splits_dir - self.reuse_splits: bool = reuse_splits self.gc_model: Path | None = Path(gc_model) if gc_model else None # Genome-wide mean GC bias weight, computed once at the runner level when # gc_model is loaded. cover_dataset divides per-chunk reads by this rather @@ -251,8 +245,6 @@ def from_cli(output_dir: Path, 'overwrite_output': (bool, False, None, None), 'parallel_block_size': (int, 0, None, None), 'threads': (int, 1, 1, 1000), - 'cleanup_splits': (bool, True, None, None), - 'reuse_splits': (bool, False, None, None), 'gc_model': (Path, None, 'exists', None) } @@ -312,6 +304,8 @@ def check_and_log_error(keyname, value_to_check, crit1, crit2): # Map: deprecated key -> short reason shown to the user. DEPRECATED_KEYS = { "parallel_mode": "splitting strategy is now derived from `threads`", + "cleanup_splits": "splits are always written to a temporary directory and cleaned up automatically", + "reuse_splits": "removed; splits are regenerated on each run", } def read_yaml(self, config_yaml: Path, args: dict): @@ -467,23 +461,7 @@ def log_configuration(self): else: _LOG.info('Single threading - 1 thread.') _LOG.info('Splitting input by contig.') - if self.reuse_splits: - splits_dir = Path(f'{self.output_dir}/splits/') - _LOG.info(f'Reusing existing splits {splits_dir}.') - if not splits_dir.is_dir(): - raise FileNotFoundError(f"reuse_splits=True but splits dir not found: {splits_dir}") - else: - if self.reuse_splits: - raise FileNotFoundError(f'reuse_splits=True') - else: - _LOG.warning(f'Reused splits set to True, but splits dir not found: {splits_dir}. Creating new splits') - _LOG.info(f'Preserving splits for next run in directory {self.splits_dir}.') - elif not self.cleanup_splits: - splits_dir = Path(f'{self.output_dir}/splits/') - _LOG.info(f'Preserving splits for next run in directory {self.splits_dir}.') - else: - splits_dir = self.temp_dir_path / "splits" - + splits_dir = self.temp_dir_path / "splits" validate_output_path(splits_dir, False) self.splits_dir = splits_dir diff --git a/pyproject.toml b/pyproject.toml index 6b5c73d7..f4581275 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "neat-genreads" -version = "4.5.0" +version = "4.5.1" description = "NGS Simulation toolkit" readme = "README.md" authors = ["Joshua Allen ", "Keshav Gandhi "] diff --git a/tests/test_compare_vcfs/test_integration.py b/tests/test_compare_vcfs/test_integration.py index cc903a16..8bb0b7b8 100644 --- a/tests/test_compare_vcfs/test_integration.py +++ b/tests/test_compare_vcfs/test_integration.py @@ -71,7 +71,6 @@ def _write_config(path: Path, ref_path: Path) -> Path: "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 @@ -250,8 +249,7 @@ def test_compare_vcfs_real_happy_with_chrom_mismatched_bed(tmp_path, happy_bin, "rng_seed: 42\n" "mutation_rate: 0.01\n" f"mutation_bed: {bed}\n" - "overwrite_output: true\n" - "cleanup_splits: true\n", + "overwrite_output: true\n", encoding="utf-8", ) read_simulator_runner(str(cfg_path), str(sim_out), "run") diff --git a/tests/test_read_simulator/test_gc_bias_pipeline.py b/tests/test_read_simulator/test_gc_bias_pipeline.py index d0c14ffb..7cd244d3 100644 --- a/tests/test_read_simulator/test_gc_bias_pipeline.py +++ b/tests/test_read_simulator/test_gc_bias_pipeline.py @@ -27,7 +27,6 @@ def _write_config(path: Path, ref_path: Path, **overrides) -> Path: "coverage": 2, "rng_seed": 42, "overwrite_output": "true", - "cleanup_splits": "true", } defaults.update(overrides) lines = "\n".join(f"{k}: {v}" for k, v in defaults.items()) diff --git a/tests/test_read_simulator/test_models_pipeline.py b/tests/test_read_simulator/test_models_pipeline.py index 6706b9f9..82b95f97 100644 --- a/tests/test_read_simulator/test_models_pipeline.py +++ b/tests/test_read_simulator/test_models_pipeline.py @@ -26,7 +26,6 @@ def _write_config(path: Path, ref_path: Path, **overrides) -> Path: "coverage": 2, "rng_seed": 42, "overwrite_output": "true", - "cleanup_splits": "true", } defaults.update(overrides) lines = "\n".join(f"{k}: {v}" for k, v in defaults.items()) diff --git a/tests/test_read_simulator/test_options.py b/tests/test_read_simulator/test_options.py index 6a82d5d6..86ebf78d 100644 --- a/tests/test_read_simulator/test_options.py +++ b/tests/test_read_simulator/test_options.py @@ -75,8 +75,6 @@ def test_from_cli_single_end_with_threads_and_splits(tmp_path: _PathAlias): parallel_block_size: 500000 threads: 2 - cleanup_splits: false - reuse_splits: false """ ).strip() + "\n" @@ -102,8 +100,8 @@ def test_from_cli_single_end_with_threads_and_splits(tmp_path: _PathAlias): assert opts.vcf is None assert opts.threads == 2 - assert opts.splits_dir == outdir / "splits" assert opts.splits_dir.is_dir() + assert opts.splits_dir.name == "splits" def test_from_cli_paired_end_fragments(tmp_path: _PathAlias): @@ -125,8 +123,6 @@ def test_from_cli_paired_end_fragments(tmp_path: _PathAlias): overwrite_output: true threads: 1 - cleanup_splits: true - reuse_splits: false """ ).strip() + "\n" @@ -160,8 +156,6 @@ def test_default_values(): # auto-tune logic. The splitting strategy itself is no longer a user-facing option — # it's derived from `threads` at runtime. assert opts.parallel_block_size == 0 - assert opts.cleanup_splits is True - assert opts.reuse_splits is False assert opts.overwrite_output is False assert opts.rescale_qualities is False assert opts.min_mutations == 0 @@ -286,6 +280,27 @@ def test_read_yaml_deprecated_parallel_mode_warns(tmp_path: _PathAlias, caplog): ) +def test_deprecated_cleanup_splits_warns(tmp_path: _PathAlias, caplog): + """cleanup_splits and reuse_splits must fire a deprecation warning and not crash.""" + ref = _project_root() / "data" / "H1N1.fa" + cfg = _textwrap.dedent( + f""" + reference: {ref} + cleanup_splits: true + reuse_splits: false + threads: 1 + """ + ).strip() + "\n" + yml_path = tmp_path / "old_config.yml" + yml_path.write_text(cfg, encoding="utf-8") + import logging as _logging + with caplog.at_level(_logging.WARNING): + Options.from_cli(tmp_path, "out", yml_path) + keys_warned = {rec.message.split("`")[1] for rec in caplog.records if "deprecated" in rec.message} + assert "cleanup_splits" in keys_warned + assert "reuse_splits" in keys_warned + + def test_log_configuration_fragment_mean_less_than_read_len_exits(tmp_path: _PathAlias): ref = _project_root() / "data" / "H1N1.fa" opts = Options(reference=ref, output_dir=tmp_path, output_prefix="out", @@ -313,27 +328,3 @@ def test_log_configuration_paired_without_model_or_mean_exits(tmp_path: _PathAli opts.log_configuration() -def test_from_cli_reuse_splits_missing_dir_raises(tmp_path: _PathAlias): - cfg = _textwrap.dedent( - f""" - reference: {(_project_root() / 'data' / 'H1N1.fa').as_posix()} - paired_ended: false - produce_fastq: true - produce_bam: false - produce_vcf: false - threads: 4 - parallel_block_size: 500000 - cleanup_splits: true - reuse_splits: true - overwrite_output: true - """ - ).strip() + "\n" - - yml_path = tmp_path / "neat_from_cli_reuse.yml" - yml_path.write_text(cfg, encoding="utf-8") - - outdir = tmp_path / "out" - outdir.mkdir(parents=True, exist_ok=True) - - with _pytest.raises(FileNotFoundError, match=r"reuse_splits=True"): - Options.from_cli(outdir, "reuse", yml_path) diff --git a/tests/test_read_simulator/test_runner.py b/tests/test_read_simulator/test_runner.py index d01c71d4..13d51086 100644 --- a/tests/test_read_simulator/test_runner.py +++ b/tests/test_read_simulator/test_runner.py @@ -263,7 +263,6 @@ def _write_config(path: Path, ref_path: Path, **overrides) -> Path: "coverage": 2, "rng_seed": 42, "overwrite_output": "true", - "cleanup_splits": "true", } defaults.update(overrides) lines = "\n".join(f"{k}: {v}" for k, v in defaults.items()) From 0a7247bce32e6e1753730a1a23338011b19e4f61 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 23:13:18 -0500 Subject: [PATCH 06/12] Drop cleanup_splits and reuse_splits; bump to v4.5.1 reuse_splits was broken (raised unconditionally regardless of whether the splits dir existed). cleanup_splits existed only to support it. Splits always go in the TemporaryDirectory that auto-cleans on exit. Both keys added to DEPRECATED_KEYS so old configs warn instead of silently no-oping. README examples, test fixtures, and test_options updated. Co-Authored-By: Claude Sonnet 4.6 --- config_template/simple_template.yml | 2 -- config_template/template_neat_config.yml | 13 +------------ 2 files changed, 1 insertion(+), 14 deletions(-) diff --git a/config_template/simple_template.yml b/config_template/simple_template.yml index 9eab62b4..e7a6ac27 100644 --- a/config_template/simple_template.yml +++ b/config_template/simple_template.yml @@ -29,5 +29,3 @@ overwrite_output: . parallel_block_size: . threads: . -cleanup_splits: . -reuse_splits: . diff --git a/config_template/template_neat_config.yml b/config_template/template_neat_config.yml index 342ccbb2..ac622704 100644 --- a/config_template/template_neat_config.yml +++ b/config_template/template_neat_config.yml @@ -1,4 +1,4 @@ -## Template for NEAT's read-simulator (as of version 4.3.5, parallelization-friendly) +## Template for NEAT's read-simulator (as of version 4.5.1, parallelization-friendly) ## Any parameter that is not required but has a default value will use the ## default value even if the variable is not included in the config. For ## required items, they must be included in the config and they must be given a value. @@ -152,14 +152,3 @@ parallel_block_size: . # Maximum number of concurrent NEAT jobs (threads or hyperthreads) to run # type = int | required: no | default = all available threads: . - -# Delete the 'splits' directory after stitching completes -# Note: If threads == 1, this option has no effect. -# type = bool | required: no | default = true -cleanup_splits: . - -# Reuse existing files in '/splits' and skip the split step. -# The directory must contain NEAT-generated files and must be in the output directory within "splits" -# Note: If threads == 1, this option has no effect. -# type = bool | required: no | default = False -reuse_splits: . From b82ba8f1411e92b16f38b65ab9a1aad3e4277916 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 23:26:10 -0500 Subject: [PATCH 07/12] Add 'Validating NEAT outputs' section to README Replaces the NEAT 2.1 validateFQ/validateBam scripts (Python 2, broken read-name assumption, no gzip support, BAM tool was a raw hex dump not a validator) with recommended standard tools: FastQC/fastp for FASTQ, samtools quickcheck/flagstat and Picard ValidateSamFile for BAM, and bcftools stats for VCF. Co-Authored-By: Claude Sonnet 4.6 --- README.md | 56 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/README.md b/README.md index 73cb783f..eced950f 100755 --- a/README.md +++ b/README.md @@ -664,6 +664,62 @@ 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. +## Validating NEAT Outputs + +NEAT does not ship its own validation utilities. Use the standard bioinformatics +tools below, which handle gzipped inputs, are actively maintained, and are +typically already present in any conda/bioconda environment. + +### FASTQ validation + +**FastQC** — format validation plus quality metrics (GC content, adapter +contamination, per-base quality scores): + +```bash +fastq read1.fastq.gz # single-end +fastqc read1.fastq.gz read2.fastq.gz # paired-end +``` + +**fastp** — lightweight format check without trimming or filtering: + +```bash +fastp --in1 read1.fastq.gz --in2 read2.fastq.gz \ + --disable_adapter_trimming --disable_quality_filtering \ + --disable_length_filtering --thread 4 +``` + +### BAM validation + +**samtools quickcheck** — fast EOF and header sanity check; exits non-zero on +failure, making it CI-friendly: + +```bash +samtools quickcheck output.bam +``` + +**samtools flagstat** — alignment statistics; a non-zero exit or obviously wrong +counts (e.g., 0 mapped reads when coverage > 0) indicate a problem: + +```bash +samtools flagstat output.bam +``` + +**Picard ValidateSamFile** — comprehensive structural validation (CIGAR +consistency, mate-pair pairing, flag conflicts): + +```bash +picard ValidateSamFile -I output.bam -MODE SUMMARY +``` + +### VCF validation + +**bcftools stats** — reports site counts, ts/tv ratio, and indel length +distribution; useful for a sanity check against expected mutation rates: + +```bash +bcftools stats output_golden.vcf.gz | grep "^SN" +``` + ## Tests We provide unit tests (e.g., mutation and sequencing error models) and basic integration tests for the CLI. From 858bd503060c4e406810f191f3f5ff992e4c1c98 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 19 May 2026 23:51:50 -0500 Subject: [PATCH 08/12] Fix seven doc gaps found in audit; update header to v4.5 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Version headers: README title/intro/ToC updated from v4.4 → v4.5 - read_len default: corrected from 101 → 151 (matches from_cli defs) - coverage type: corrected from "float or integer" → "integer only" - threads: removed duplicate row from config table - off_target_scalar: removed from template_neat_config.yml (key has never existed in options.py defs; silently ignored since at least v4.3) - Added three undocumented config keys to README table: quality_offset, overwrite_output, no_coverage_bias - Added neat bacterial-wrapper section to Utilities + ToC entry Co-Authored-By: Claude Sonnet 4.6 --- README.md | 41 ++++++++++++++++++------ config_template/template_neat_config.yml | 6 ---- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index eced950f..28a8c680 100755 --- a/README.md +++ b/README.md @@ -1,12 +1,12 @@ -# The NEAT Project v4.4 +# The NEAT Project v4.5 -Welcome to the NEAT project, the NExt-generation sequencing Analysis Toolkit, version 4.4. NEAT 4.4 is the official release of NEAT 4.0. It represents a lot of hard work from several contributors at NCSA and beyond. With the addition of parallel processing, we feel that the code is ready for production, and future releases will focus on compatibility, bug fixes, and testing. Future releases for the time being will be enumerations of 4.4.X. +Welcome to the NEAT project, the NExt-generation sequencing Analysis Toolkit, version 4.5. See the [ChangeLog](ChangeLog.md) for the full version history. -## NEAT v4.4 +## NEAT v4.5 -NEAT 4.4 is the current official release of NEAT 4.0, including parallel processing support and significant bug fixes to the sequencing error model. See the [ChangeLog](ChangeLog.md) for details. +NEAT 4.5 is the current release. It adds the `neat compare-vcfs` subcommand, removes the deprecated `cleanup_splits`/`reuse_splits` config keys, and ships a number of performance and correctness fixes accumulated across the 4.4.x line. See the [ChangeLog](ChangeLog.md) for details. -We have completed major revisions on NEAT since 3.4 and consider NEAT 4.4 to be a stable release, in that we will continue to update and provide bug fixes and support. We will consider new features and pull requests. Please include justification for major changes. See [contribute](CONTRIBUTING.md) for more information. If you'd like to use some of our code in your own, no problem! Just review the [license](LICENSE.md), first. +We consider NEAT 4.x to be a stable release line and will continue to provide bug fixes and support. We will consider new features and pull requests. Please include justification for major changes. See [contribute](CONTRIBUTING.md) for more information. If you'd like to use some of our code in your own, please review the [license](LICENSE.md) first. We've deprecated NEAT's command-line interface options for the most part, opting to simplify things with configuration files. If you require the CLI for legacy purposes, NEAT 3.4 was our last release to be fully supported via command-line interface. Please convert your CLI commands to the corresponding configuration file for future runs. @@ -26,8 +26,8 @@ To cite this work, please use both of the following: ## Table of Contents -* [The NEAT Project v4.4](#the-neat-project-v44) -* [NEAT v4.4](#neat-v44) +* [The NEAT Project v4.5](#the-neat-project-v45) +* [NEAT v4.5](#neat-v45) * [Table of Contents](#table-of-contents) * [Prerequisites](#prerequisites) * [Installation](#installation) @@ -48,6 +48,8 @@ To cite this work, please use both of the following: * [`neat model-qual-score`](#neat-model-qual-score) * [`neat model-gc-bias`](#neat-model-gc-bias) * [`neat compare-vcfs`](#neat-compare-vcfs) + * [`neat bacterial-wrapper`](#neat-bacterial-wrapper) + * [Validating NEAT Outputs](#validating-neat-outputs) * [Tests](#tests) * [Guide to run locally](#guide-to-run-locally) * [Note on Sensitive Patient Data](#note-on-sensitive-patient-data) @@ -165,9 +167,9 @@ To run the simulator in multithreaded mode, set the `threads` value in the confi `reference`: Full path to a FASTA file to generate reads from. -`read_len`: The length of the reads for the FASTQ (if using). _Integer value, default 101._ +`read_len`: The length of the reads for the FASTQ (if using). _Integer value, default 151._ -`coverage`: Desired coverage value. _Float or integer, default = 10._ +`coverage`: Desired coverage depth. _Integer value, default = 10._ `ploidy`: Desired value for ploidy (# of copies of each chromosome in the organism, where if ploidy > 2, "heterozygous" mutates floor(ploidy / 2) chromosomes). _Default is 2._ @@ -193,9 +195,9 @@ More parameters are below: | `mutation_model` | Full path to a mutation model generated by NEAT. Leave empty to use a default model (default model based on human data sequenced by Illumina). | | `fragment_model` | Full path to fragment length model generated by NEAT. Leave empty to use default model (default model based on human data sequenced by Illumina). | | `gc_model` | Full path to GC-bias model generated by NEAT. Leave empty for no GC bias. | -| `threads` | The number of threads for NEAT to use. Increasing the number will speed up read generation. | | `avg_seq_error` | Average sequencing error rate for the sequencing machine. Use to increase or decrease the rate of errors in the reads. Float between 0 and 0.3. Default is set by the error model. | | `rescale_qualities` | Rescale the quality scores to reflect the `avg_seq_error` rate above. Set `True` to activate if you notice issues with the sequencing error rates in your dataset. | +| `quality_offset` | ASCII offset for quality score encoding. Default `33` (Sanger/Illumina 1.8+ Phred+33). Only change this if your data uses a non-standard encoding (valid range 33–64). | | `include_vcf` | Full path to list of variants in VCF format to include in the simulation. These will be inserted as they appear in the input VCF into the final VCF, and the corresponding FASTQ and BAM files, if requested. | | `target_bed` | Full path to list of regions in BED format to target. All areas outside these regions will have coverage of 0. | | `discard_bed` | Full path to a list of regions to discard, in BED format. | @@ -203,6 +205,8 @@ More parameters are below: | `mutation_bed` | Full path to a list of regions with a column describing the mutation rate of that region, as a float with values between 0 and 0.3. The mutation rate must be in the third column as, e.g., `mut_rate`=0.00. | | `rng_seed` | Manually enter a seed for the random number generator. Used for repeating runs. Must be an integer. | | `min_mutations` | Set the minimum number of mutations that NEAT should add, per contig. Default is 0. We recommend setting this to at least one for small chromosomes, so NEAT will produce at least one mutation per contig. | +| `overwrite_output` | If `true`, existing output files with the same name will be overwritten. Default `false` (NEAT exits with an error if output files already exist). | +| `no_coverage_bias` | If `true`, disables statistical coverage models and forces uniform coverage across the genome. Intended for debugging. Default `false`. | | `threads` | Number of threads to use. More than 1 will use multi-threading to speed up processing. With `threads > 1`, NEAT splits each contig into chunks; with `threads == 1`, one chunk per contig is used. | | `parallel_block_size` | Per-chunk size in bases when `threads > 1`. Default `0` (auto-tune from total genome length and thread count, targeting ~8 chunks per thread). Set to a positive integer to override. Ignored when `threads == 1`. | @@ -664,6 +668,23 @@ 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. +### `neat bacterial-wrapper` + +Runs NEAT's read simulator twice on the same reference — once with a standard +mutation model ("Regular") and once with a higher mutation rate model +("Wrapped") representing a bacterium under selection pressure — then stitches +the two output sets together for downstream comparison. + +```bash +neat bacterial-wrapper reference.fa bacteria_name \ + -c config.yml \ + -o /path/to/output/dir +``` + +Outputs are written into `/Regular/` and `/Wrapped/` +subdirectories, then stitched into the parent output directory. The config file +follows the same format as `neat read-simulator`. + ## Validating NEAT Outputs NEAT does not ship its own validation utilities. Use the standard bioinformatics diff --git a/config_template/template_neat_config.yml b/config_template/template_neat_config.yml index ac622704..7692569b 100644 --- a/config_template/template_neat_config.yml +++ b/config_template/template_neat_config.yml @@ -57,12 +57,6 @@ include_vcf: . # type = string | required = no target_bed: . -# Scalar value for coverage in regions outside the targeted BED. Example: 0.5 -# would get you roughly half the coverage as the on target areas. Default is -# 0 coverage in off-target regions. Number should be a float in decimal -# type: float | required = no | default = 0.00 -off_target_scalar: . - # Absolute path to BED file containing reference regions that the simulation should discard # type = string | required = no discard_bed: . From 1be27edbc6a91de3744448dc4ce11cb31431061c Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Wed, 20 May 2026 00:57:33 -0500 Subject: [PATCH 09/12] Implement binned quality scoring for Illumina instruments (#293) TraditionalQualityModel now accepts an optional quality_bins parameter that snaps sampled scores to the nearest discrete bin using vectorized np.searchsorted, matching how NovaSeq/NextSeq 2000 instruments emit only a small set of Phred values (e.g. Q2/Q12/Q23/Q37). The runner passes allowed_bins through so both Markov and traditional models produce realistically binned output when bins are specified. - Add quality_bins to TraditionalQualityModel.__init__ and get_quality_scores() - Pass allowed_bins to TraditionalQualityModel in model_qual_score runner - Expose _down_bin_quality as public down_bin_quality in markov_utils - Add QUALITY_PRESETS dict in neat/quality_score_modeling/presets.py (novaseq, nextseq2000, nextseq500) - Add --quality-preset CLI flag to neat model-qual-score; implies --markov - Fix stale -Q help text in model_sequencing_error CLI - Add TraditionalQualityModel binning tests + preset data/CLI tests Co-Authored-By: Claude Sonnet 4.6 --- neat/cli/commands/model_qual_score.py | 28 +++++-- neat/cli/commands/model_sequencing_error.py | 5 +- neat/model_quality_score/runner.py | 1 + neat/models/error_models.py | 15 +++- neat/quality_score_modeling/markov_utils.py | 5 +- neat/quality_score_modeling/presets.py | 20 +++++ tests/test_models/test_markov_utils.py | 20 ++--- tests/test_models/test_qual_score_models.py | 86 ++++++++++++++++++++- tests/test_models/test_quality_presets.py | 71 +++++++++++++++++ 9 files changed, 229 insertions(+), 22 deletions(-) create mode 100644 neat/quality_score_modeling/presets.py create mode 100644 tests/test_models/test_quality_presets.py diff --git a/neat/cli/commands/model_qual_score.py b/neat/cli/commands/model_qual_score.py index aea15edf..422aa15d 100644 --- a/neat/cli/commands/model_qual_score.py +++ b/neat/cli/commands/model_qual_score.py @@ -1,6 +1,7 @@ import argparse from ...model_quality_score import model_qual_score_runner +from ...quality_score_modeling.presets import QUALITY_PRESETS from .base import BaseCommand from .options import output_group @@ -38,7 +39,20 @@ def add_arguments(self, parser: argparse.ArgumentParser): type=int, nargs="+", default=[42], - help="Max quality or explicit list of quality scores [42].", + help="Maximum quality score (single int) or explicit list of bin values " + "(e.g. -Q 2 12 23 37). A list enables Markov binning: all observed " + "scores are down-binned to the nearest value and simulation output is " + "constrained to those values. Overridden by --quality-preset. [42]", + ) + + parser.add_argument( + "--quality-preset", + dest="quality_preset", + choices=list(QUALITY_PRESETS), + default=None, + metavar="PRESET", + help="Named bin preset for common Illumina instruments. Implies --markov. " + f"Choices: {', '.join(QUALITY_PRESETS)}. Overrides -Q when set.", ) parser.add_argument( @@ -69,11 +83,15 @@ def add_arguments(self, parser: argparse.ArgumentParser): def execute(self, arguments: argparse.Namespace): - if len(arguments.quality_scores) == 1: - qual_scores: int | list[int] = arguments.quality_scores[0] - + if arguments.quality_preset: + qual_scores: int | list[int] = QUALITY_PRESETS[arguments.quality_preset] + use_markov = True + elif len(arguments.quality_scores) == 1: + qual_scores = arguments.quality_scores[0] + use_markov = arguments.use_markov else: qual_scores = arguments.quality_scores + use_markov = arguments.use_markov model_qual_score_runner( files=arguments.input_files, @@ -83,5 +101,5 @@ def execute(self, arguments: argparse.Namespace): overwrite=arguments.overwrite, output_dir=arguments.output_dir, output_prefix=arguments.prefix, - use_markov=arguments.use_markov, + use_markov=use_markov, ) diff --git a/neat/cli/commands/model_sequencing_error.py b/neat/cli/commands/model_sequencing_error.py index 4c3420b0..c9527219 100644 --- a/neat/cli/commands/model_sequencing_error.py +++ b/neat/cli/commands/model_sequencing_error.py @@ -49,8 +49,9 @@ def add_arguments(self, parser: argparse.ArgumentParser): nargs='+', required=False, default=42, - help="Quality score max. The default 42. The lowest possible score is 1. To used binned" - "scoring, enter a space separated list of scores, e.g., -Q 2 15 23 37") + help="Maximum quality score [42], or a space-separated list of bin values " + "(e.g., -Q 2 12 23 37) to constrain the model to those discrete levels. " + "For named instrument presets use `neat model-qual-score --quality-preset`.") parser.add_argument('-m', type=int, diff --git a/neat/model_quality_score/runner.py b/neat/model_quality_score/runner.py index 2c355c7b..0e9ca77a 100644 --- a/neat/model_quality_score/runner.py +++ b/neat/model_quality_score/runner.py @@ -133,6 +133,7 @@ def model_qual_score_runner( average_error=average_error, quality_scores=np.array(final_quality_scores), qual_score_probs=read_parameters[idx], + quality_bins=allowed_bins, ) markov_model: Optional[MarkovQualityModel] = None diff --git a/neat/models/error_models.py b/neat/models/error_models.py index c9ed0793..c3e44411 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -49,7 +49,8 @@ def __init__( transition_matrix: np.ndarray = default_error_transition_matrix, quality_scores: np.ndarray = default_quality_scores, qual_score_probs: np.ndarray = default_qual_score_probs, - is_uniform: bool = False + is_uniform: bool = False, + quality_bins: list[int] | None = None, ): self.transition_matrix = transition_matrix @@ -57,6 +58,10 @@ def __init__( self.quality_score_probabilities = qual_score_probs self.is_uniform = is_uniform self.average_error = average_error + # When set, generated scores are snapped to the nearest bin value ≤ the + # drawn score (down-binning), matching the discrete quality levels produced + # by instruments like NovaSeq (Q2/Q12/Q23/Q37). + self.quality_bins: list[int] | None = sorted(quality_bins) if quality_bins else None # pre-compute the error rate for each quality score. This is the inverse of the phred score equation self.quality_score_error_rate: dict[int, float] = {x: 10. ** (-x / 10) for x in self.quality_scores} @@ -101,7 +106,13 @@ def get_quality_scores( means = self.quality_score_probabilities[quality_index_map, 0] scales = self.quality_score_probabilities[quality_index_map, 1] scores = rng.normal(means, scales) - return np.clip(np.rint(scores).astype(int), 1, 42) + scores = np.clip(np.rint(scores).astype(int), 1, 42) + if self.quality_bins: + bins = np.array(self.quality_bins) + idx = np.searchsorted(bins, scores, side="right") - 1 + idx = np.clip(idx, 0, len(bins) - 1) + scores = bins[idx] + return scores class SequencingErrorModel(SnvModel, DeletionModel, InsertionModel): diff --git a/neat/quality_score_modeling/markov_utils.py b/neat/quality_score_modeling/markov_utils.py index d75c7a47..75f9088e 100644 --- a/neat/quality_score_modeling/markov_utils.py +++ b/neat/quality_score_modeling/markov_utils.py @@ -17,6 +17,7 @@ _LOG = logging.getLogger(__name__) __all__ = [ + "down_bin_quality", "read_quality_lists", "compute_initial_distribution", "compute_position_distributions", @@ -25,7 +26,7 @@ ] -def _down_bin_quality(q: int, allowed: List[int]) -> int: +def down_bin_quality(q: int, allowed: List[int]) -> int: """ Map q to the greatest allowed value <= q (down-binning). If q is below the smallest allowed, map to allowed[0]. @@ -110,7 +111,7 @@ def read_quality_lists( continue if allowed_sorted is not None: - qlist = [_down_bin_quality(q, allowed_sorted) for q in qlist] + qlist = [down_bin_quality(q, allowed_sorted) for q in qlist] qualities.append(qlist) reads_read += 1 diff --git a/neat/quality_score_modeling/presets.py b/neat/quality_score_modeling/presets.py new file mode 100644 index 00000000..b6b1f40f --- /dev/null +++ b/neat/quality_score_modeling/presets.py @@ -0,0 +1,20 @@ +""" +Named quality-bin presets for common Illumina instruments. + +Each preset maps to the discrete Phred score levels that the instrument emits. +Pass the preset name to ``neat model-qual-score --quality-preset`` instead of +specifying raw bin values with ``-Q``. + +References: + NovaSeq 6000 / NovaSeq X: Illumina RTA3 4-level binning (Q2, Q12, Q23, Q37) + NextSeq 2000: 4-level binning matching NovaSeq X (Q2, Q12, Q26, Q37) + NextSeq 500 / MiniSeq: 5-level binning (Q2, Q12, Q23, Q27, Q37) +""" + +__all__ = ["QUALITY_PRESETS"] + +QUALITY_PRESETS: dict[str, list[int]] = { + "novaseq": [2, 12, 23, 37], + "nextseq2000": [2, 12, 26, 37], + "nextseq500": [2, 12, 23, 27, 37], +} \ No newline at end of file diff --git a/tests/test_models/test_markov_utils.py b/tests/test_models/test_markov_utils.py index 333ccf36..c45389d9 100644 --- a/tests/test_models/test_markov_utils.py +++ b/tests/test_models/test_markov_utils.py @@ -5,7 +5,7 @@ import pytest from neat.quality_score_modeling.markov_utils import ( - _down_bin_quality, + down_bin_quality, read_quality_lists, compute_initial_distribution, compute_position_distributions, @@ -27,33 +27,33 @@ def _write_fastq(path, reads): # --------------------------------------------------------------------------- -# _down_bin_quality +# down_bin_quality # --------------------------------------------------------------------------- def test_down_bin_exact_match(): - assert _down_bin_quality(30, [10, 20, 30, 40]) == 30 + assert down_bin_quality(30, [10, 20, 30, 40]) == 30 def test_down_bin_between_bins_maps_down(): - assert _down_bin_quality(25, [10, 20, 30, 40]) == 20 + assert down_bin_quality(25, [10, 20, 30, 40]) == 20 def test_down_bin_below_min_maps_to_first_bin(): - assert _down_bin_quality(5, [10, 20, 30]) == 10 + assert down_bin_quality(5, [10, 20, 30]) == 10 def test_down_bin_above_max_maps_to_last_bin(): - assert _down_bin_quality(99, [10, 20, 30]) == 30 + assert down_bin_quality(99, [10, 20, 30]) == 30 def test_down_bin_empty_allowed_returns_q_unchanged(): - assert _down_bin_quality(25, []) == 25 + assert down_bin_quality(25, []) == 25 def test_down_bin_single_bin(): - assert _down_bin_quality(0, [20]) == 20 - assert _down_bin_quality(20, [20]) == 20 - assert _down_bin_quality(40, [20]) == 20 + assert down_bin_quality(0, [20]) == 20 + assert down_bin_quality(20, [20]) == 20 + assert down_bin_quality(40, [20]) == 20 # --------------------------------------------------------------------------- diff --git a/tests/test_models/test_qual_score_models.py b/tests/test_models/test_qual_score_models.py index 4e4eef9e..192e91ea 100644 --- a/tests/test_models/test_qual_score_models.py +++ b/tests/test_models/test_qual_score_models.py @@ -1,11 +1,14 @@ """ -Unit tests for MarkovQualityModel +Unit tests for MarkovQualityModel and TraditionalQualityModel (binned scoring). """ +import gzip +import pickle import pytest import numpy as np from numpy.random import default_rng +from neat.models.error_models import TraditionalQualityModel from neat.models.markov_quality_model import MarkovQualityModel # --------------------------------------------------------------------------- @@ -164,3 +167,84 @@ def test_position_index_for_length_single_position_read(): """length=1 always returns index 0 regardless of pos.""" qm = _simple_model(read_length=151) assert qm._position_index_for_length(0, 1) == 0 + + +# =========================================================================== +# TraditionalQualityModel — binned scoring +# =========================================================================== + +def _trad_model(quality_bins=None): + """Minimal TraditionalQualityModel with a flat μ=30, σ=5 distribution.""" + read_length = 151 + qual_score_probs = np.full((read_length, 2), [30.0, 5.0]) + quality_scores = np.arange(0, 43) + return TraditionalQualityModel( + quality_scores=quality_scores, + qual_score_probs=qual_score_probs, + quality_bins=quality_bins, + ) + + +def test_traditional_model_unbinned_produces_varied_scores(): + """Without bins, scores span a range (not locked to a small set).""" + rng = default_rng(42) + model = _trad_model() + scores = model.get_quality_scores(151, 151, rng) + assert len(scores) == 151 + assert len(set(scores.tolist())) > 4 + + +def test_traditional_model_binned_output_constrained(): + """All output scores must be members of the supplied bin set.""" + bins = [2, 12, 23, 37] + rng = default_rng(42) + model = _trad_model(quality_bins=bins) + scores = model.get_quality_scores(151, 151, rng) + assert set(scores.tolist()).issubset(set(bins)) + + +def test_traditional_model_binned_novaseq_bins(): + """NovaSeq preset bins produce only {2, 12, 23, 37} across many reads.""" + bins = [2, 12, 23, 37] + rng = default_rng(7) + model = _trad_model(quality_bins=bins) + all_scores = set() + for _ in range(20): + all_scores.update(model.get_quality_scores(151, 151, rng).tolist()) + assert all_scores.issubset(set(bins)) + + +def test_traditional_model_binned_below_min_maps_to_first_bin(): + """A score below the lowest bin should map to the first (lowest) bin.""" + # Force very low scores: μ=1, σ=0.1 → always clips to 1 → below Q2 → Q2 + read_length = 10 + qual_score_probs = np.full((read_length, 2), [1.0, 0.1]) + model = TraditionalQualityModel( + quality_scores=np.arange(0, 43), + qual_score_probs=qual_score_probs, + quality_bins=[2, 12, 23, 37], + ) + rng = default_rng(0) + scores = model.get_quality_scores(read_length, read_length, rng) + assert set(scores.tolist()) == {2} + + +def test_traditional_model_binned_persists_through_pickle(tmp_path): + """quality_bins must survive a pickle/unpickle round-trip.""" + bins = [2, 12, 23, 37] + model = _trad_model(quality_bins=bins) + path = tmp_path / "model.pkl" + with open(path, "wb") as f: + pickle.dump(model, f) + with open(path, "rb") as f: + loaded = pickle.load(f) + assert loaded.quality_bins == bins + rng = default_rng(1) + scores = loaded.get_quality_scores(151, 151, rng) + assert set(scores.tolist()).issubset(set(bins)) + + +def test_traditional_model_none_bins_is_unbinned(): + """Passing quality_bins=None (explicit) is identical to the default.""" + model = _trad_model(quality_bins=None) + assert model.quality_bins is None diff --git a/tests/test_models/test_quality_presets.py b/tests/test_models/test_quality_presets.py new file mode 100644 index 00000000..d234a878 --- /dev/null +++ b/tests/test_models/test_quality_presets.py @@ -0,0 +1,71 @@ +""" +Tests for neat/quality_score_modeling/presets.py and the --quality-preset +CLI wiring in neat/cli/commands/model_qual_score.py. +""" +import argparse + +import pytest + +from neat.quality_score_modeling.presets import QUALITY_PRESETS + + +# =========================================================================== +# Preset data integrity +# =========================================================================== + +def test_all_presets_are_sorted(): + for name, bins in QUALITY_PRESETS.items(): + assert bins == sorted(bins), f"preset {name!r} bins are not sorted" + + +def test_all_preset_bins_are_positive(): + for name, bins in QUALITY_PRESETS.items(): + assert all(b > 0 for b in bins), f"preset {name!r} has non-positive bin" + + +def test_novaseq_has_four_bins(): + assert len(QUALITY_PRESETS["novaseq"]) == 4 + + +def test_nextseq2000_has_four_bins(): + assert len(QUALITY_PRESETS["nextseq2000"]) == 4 + + +def test_nextseq500_has_five_bins(): + assert len(QUALITY_PRESETS["nextseq500"]) == 5 + + +def test_novaseq_bins_match_illumina_spec(): + assert QUALITY_PRESETS["novaseq"] == [2, 12, 23, 37] + + +def test_nextseq2000_bins_match_illumina_spec(): + assert QUALITY_PRESETS["nextseq2000"] == [2, 12, 26, 37] + + +# =========================================================================== +# CLI --quality-preset wiring +# =========================================================================== + +def _parse_qual_score_args(argv): + """Run the model-qual-score argument parser against argv.""" + from neat.cli.commands.model_qual_score import Command + parser = argparse.ArgumentParser() + Command(parser) + return parser.parse_args(argv) + + +def test_quality_preset_accepted_by_cli(): + args = _parse_qual_score_args(["-i", "f.fq", "-o", "/tmp", "--quality-preset", "novaseq"]) + assert args.quality_preset == "novaseq" + + +def test_unknown_preset_rejected_by_cli(): + with pytest.raises(SystemExit): + _parse_qual_score_args(["-i", "f.fq", "-o", "/tmp", "--quality-preset", "bad_preset"]) + + +def test_all_preset_names_accepted_by_cli(): + for name in QUALITY_PRESETS: + args = _parse_qual_score_args(["-i", "f.fq", "-o", "/tmp", "--quality-preset", name]) + assert args.quality_preset == name \ No newline at end of file From a46c227da2c7338b4e4d61c0393fca916f9bf43f Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Wed, 20 May 2026 00:58:28 -0500 Subject: [PATCH 10/12] docs: document --quality-preset and binned quality scoring in README Add a "Binned quality scoring" subsection to neat model-qual-score showing how to use --quality-preset (novaseq, nextseq2000, nextseq500) and explicit -Q bin lists for instruments that emit discrete Phred levels. Co-Authored-By: Claude Sonnet 4.6 --- README.md | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/README.md b/README.md index 28a8c680..2b676ee1 100755 --- a/README.md +++ b/README.md @@ -586,6 +586,42 @@ Similarly, use `-i2` to produce a model for paired-ended data. `-q` denotes the Finally, `-o` is the output directory for the model file and `-p` is the prefix for the output model, such that the file will be written as `.p.gz` inside the output folder. +#### Binned quality scoring for modern Illumina instruments + +Modern Illumina instruments (NovaSeq 6000, NovaSeq X, NextSeq 2000) compress +Phred quality scores into a small discrete set of bins rather than emitting a +continuous range. To train a model that faithfully reproduces this behaviour, +use either the named `--quality-preset` flag or an explicit bin list via `-Q`. + +**Named presets** (recommended): + +```bash +# NovaSeq 6000 / NovaSeq X — Q2, Q12, Q23, Q37 +neat model-qual-score -i reads.fastq.gz --quality-preset novaseq \ + -o /path/to/models -p novaseq_model + +# NextSeq 2000 — Q2, Q12, Q26, Q37 +neat model-qual-score -i reads.fastq.gz --quality-preset nextseq2000 \ + -o /path/to/models -p nextseq2000_model + +# NextSeq 500 / MiniSeq — Q2, Q12, Q23, Q27, Q37 +neat model-qual-score -i reads.fastq.gz --quality-preset nextseq500 \ + -o /path/to/models -p nextseq500_model +``` + +`--quality-preset` implies `--markov`; you do not need to pass both. + +**Explicit bin list** — if your instrument uses non-standard bins, pass them +directly with `-Q`: + +```bash +neat model-qual-score -i reads.fastq.gz -Q 2 12 23 37 --markov \ + -o /path/to/models -p custom_binned_model +``` + +When bins are specified, both the Markov and traditional quality models will +constrain simulation output to those exact Phred values. + ### `neat model-gc-bias` Computes GC-bias model from a BAM file and reference genome. It calculates the relative weight of fragments based on their GC content. From 3545f50aa290cd07e60391bb2ce4cdc14796abb5 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Wed, 20 May 2026 01:23:57 -0500 Subject: [PATCH 11/12] Updated notes in readme --- README.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/README.md b/README.md index 28a8c680..5b25e0a8 100755 --- a/README.md +++ b/README.md @@ -565,6 +565,8 @@ neat model-seq-err \ Please note that `-i2` can be used in place of `-i` to produce paired data. +> **Note on binned quality scores in `model-seq-err`:** the `-Q` flag accepts a space-separated list (e.g. `-Q 2 15 23 37`) and the CLI advertises this as a binning option, but binning is not yet wired up in the sequencing-error runner — only the maximum value in the list is currently used. If you need a binned quality distribution, train it with `neat model-qual-score --markov` (see below) and supply that model to `read-simulator`. + ### `neat model-qual-score` Typical usage: @@ -586,6 +588,22 @@ Similarly, use `-i2` to produce a model for paired-ended data. `-q` denotes the Finally, `-o` is the output directory for the model file and `-p` is the prefix for the output model, such that the file will be written as `.p.gz` inside the output folder. +#### Binned quality scores + +`-Q` accepts either a single integer (max Q, full 0..Q range learned) or a space-separated list of integers (bins). When `--markov` is set together with a list, the Markov quality model is trained with `allowed_quality_scores` restricted to those bins — observed qualities are snapped to the nearest allowed bin before fitting transition matrices, and the resulting model will only emit qualities drawn from that set. + +```bash +neat model-qual-score \ + -i input_reads.fastq(.gz) \ + -q 33 \ + -Q 2 11 24 37 \ + --markov \ + -o /path/to/models \ + -p my_binned_qual_model +``` + +Pass the resulting `*.p.gz` model to `read-simulator` (via the relevant model field in your config) and the simulated FASTQ qualities will be drawn from the binned distribution. There is no separate "quantize output qualities" knob in the simulator config — output binning is determined entirely by the quality model you train. Without `--markov`, the list form still trains over `0..max(bins)`; binning only takes effect on the Markov path. + ### `neat model-gc-bias` Computes GC-bias model from a BAM file and reference genome. It calculates the relative weight of fragments based on their GC content. From 90f41e45b2f6fe5ed8792309d383ceff65c366e2 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Mon, 8 Jun 2026 11:37:13 -0500 Subject: [PATCH 12/12] Updating gitignore --- .gitignore | 6 ------ AGENTS.md | 2 +- CLAUDE.md | 2 +- neat/read_simulator/utils/vcf_func.py | 5 ++++- tests/test_read_simulator/test_vcf_func.py | 4 ++++ 5 files changed, 10 insertions(+), 9 deletions(-) diff --git a/.gitignore b/.gitignore index 2452c556..cb0b0fd1 100755 --- a/.gitignore +++ b/.gitignore @@ -29,9 +29,3 @@ dist/ build/ *.egg-info/ *.egg - -# gitnexus stuff -/.gitnexus/ -/.claude/ -AGENTS.md -CLAUDE.md diff --git a/AGENTS.md b/AGENTS.md index 83bf1a02..7dc8ae5b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -1,7 +1,7 @@ # GitNexus — Code Intelligence -This project is indexed by GitNexus as **NEAT** (2742 symbols, 5450 relationships, 62 execution flows). Use the GitNexus MCP tools to understand code, assess impact, and navigate safely. +This project is indexed by GitNexus as **NEAT** (3389 symbols, 6682 relationships, 71 execution flows). Use the GitNexus MCP tools to understand code, assess impact, and navigate safely. > If any GitNexus tool warns the index is stale, run `npx gitnexus analyze` in terminal first. diff --git a/CLAUDE.md b/CLAUDE.md index 83bf1a02..7dc8ae5b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -1,7 +1,7 @@ # GitNexus — Code Intelligence -This project is indexed by GitNexus as **NEAT** (2742 symbols, 5450 relationships, 62 execution flows). Use the GitNexus MCP tools to understand code, assess impact, and navigate safely. +This project is indexed by GitNexus as **NEAT** (3389 symbols, 6682 relationships, 71 execution flows). Use the GitNexus MCP tools to understand code, assess impact, and navigate safely. > If any GitNexus tool warns the index is stale, run `npx gitnexus analyze` in terminal first. diff --git a/neat/read_simulator/utils/vcf_func.py b/neat/read_simulator/utils/vcf_func.py index 584550f0..f29fcc75 100755 --- a/neat/read_simulator/utils/vcf_func.py +++ b/neat/read_simulator/utils/vcf_func.py @@ -244,8 +244,11 @@ def parse_input_vcf( ) elif len(ref) > len(alt) and ref.startswith(alt): # type = deletion + # length is the full reference span (number of bases the deletion + # covers, including the shared anchor base), matching how deletions + # are reconstructed on output and generated by the mutation model. temp_variant = Deletion( - location, len(alt), temp_genotype, record[5], is_input=True, kwargs=data + location, len(ref), temp_genotype, record[5], is_input=True, kwargs=data ) elif len(alt) > len(ref) and alt.startswith(ref): # type = insertion diff --git a/tests/test_read_simulator/test_vcf_func.py b/tests/test_read_simulator/test_vcf_func.py index 8aa9a333..93bbb5ae 100644 --- a/tests/test_read_simulator/test_vcf_func.py +++ b/tests/test_read_simulator/test_vcf_func.py @@ -167,6 +167,10 @@ def test_parse_deletion(tmp_path, ref_fasta, empty_input_dict, opts): variants = empty_input_dict["chr1"].contig_variants.get(4, []) assert len(variants) == 1 assert isinstance(variants[0], Deletion) + # length is the full reference span (len(ref)), not len(alt). For ACGT -> A + # that is 4, so the deletion is reconstructed as ACGT -> A on output rather + # than being truncated to a single base. Regression test for issue #305. + assert variants[0].length == 4 def test_parse_insertion(tmp_path, ref_fasta, empty_input_dict, opts):