From cd589eab946d3abf6feb80b328289819e14e633a Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 18:20:42 +0200 Subject: [PATCH 1/9] Add an indexed regional-access arm for the retrieval comparison The whole-file retrieval comparison pits a SPARQL engine against cyvcf2 scanning the file, because none of Q1-Q13 is coordinate-restricted. That is internally consistent, but it measures the graph only against the one VCF access mode nobody uses for a selective question. Real VCF work seeks through a bgzip + tabix/CSI index. src/validation/regional_runner.py asks five region-restricted versions of the bioinformatic questions (R1 record count, R2 allele shapes, R3 Ti/Tv, R4 FILTER distribution, R5 per-sample genotype classes, expanded only) of every access path: cyvcf2 scanning, cyvcf2 and bcftools through the index, and the same SPARQL on QLever, Comunica, HDT and COTTAS. It reuses validation_runner's engine abstraction and its classification helpers, so a regional answer is classified by the same code as a whole-file one. Three design decisions: - A window means POS, not overlap. An index seek returns every record whose span overlaps the region, while a SPARQL ?pos filter does not; every VCF arm keeps the seek and drops records whose POS falls outside the window, so the arms agree on indels that straddle a window edge. - Windows are anchored on real records from a fixed seed, so a sparse benchmark VCF gives a selectivity ladder rather than a sparsity ladder. - The scan arm is timed on fewer windows, since it costs the same whichever window is asked; correctness is still checked on every window. Every arm's answer is compared against a cyvcf2 reference before any timing is reported; disagreements go to mismatches.json and the runner exits non-zero. The image gains the tabix package (bgzip and tabix): the Debian bcftools package links libhts but does not ship the htslib command-line tools. Queries under src/validation/queries/regional/, docs in docs/validation.md ("Indexed regional access"), 23 unit tests in test_regional_runner_unit.py. Co-Authored-By: Claude Opus 5.5 --- Dockerfile | 4 + docs/validation.md | 59 + .../common/r01_region_record_count.rq | 18 + .../common/r02_region_variant_shape_counts.rq | 41 + .../regional/common/r03_region_titv.rq | 18 + .../common/r04_region_filter_distribution.rq | 16 + .../r05_region_sample_genotype_counts.rq | 30 + src/validation/regional_runner.py | 1023 +++++++++++++++++ test/test_regional_runner_unit.py | 315 +++++ 9 files changed, 1524 insertions(+) create mode 100644 src/validation/queries/regional/common/r01_region_record_count.rq create mode 100644 src/validation/queries/regional/common/r02_region_variant_shape_counts.rq create mode 100644 src/validation/queries/regional/common/r03_region_titv.rq create mode 100644 src/validation/queries/regional/common/r04_region_filter_distribution.rq create mode 100644 src/validation/queries/regional/expanded/r05_region_sample_genotype_counts.rq create mode 100644 src/validation/regional_runner.py create mode 100644 test/test_regional_runner_unit.py diff --git a/Dockerfile b/Dockerfile index 3e811ed..ef995bb 100644 --- a/Dockerfile +++ b/Dockerfile @@ -164,6 +164,10 @@ RUN apt-get update \ python3 \ python3-venv \ raptor2-utils \ + # bgzip and tabix, for the indexed regional-access arm. bcftools alone + # cannot do it: the Debian package links libhts but does not install + # the htslib command-line tools. + tabix \ time \ && rm -rf /var/lib/apt/lists/* diff --git a/docs/validation.md b/docs/validation.md index d18bfd3..91d6afe 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -413,6 +413,65 @@ whether they resolve, and the validator reports that note if the engine cannot start. Comunica remains the default, so an image whose QLever binaries do not link still validates normally. +### Indexed regional access + +The timings above compare SPARQL against cyvcf2 **scanning the file**, because +none of Q1-Q13 is coordinate-restricted. That is internally consistent, and it +is the one VCF access mode nobody uses for a selective question: real VCF work +seeks, through a `bgzip` + `tabix` index. + +`validation/regional_runner.py` adds that arm. It asks five region-restricted +questions -- record count, allele shape, Ti/Tv, FILTER distribution, per-sample +genotype classes -- of every access path: + +| Arm | Access path | +| --- | --- | +| `cyvcf2-scan` | whole-file iteration with a POS filter (the status quo baseline) | +| `cyvcf2-indexed` | `VCF(path)(region)` over the bgzipped, indexed copy | +| `bcftools-indexed` | `bcftools query -r`, the only route to the exact FILTER text | +| `comunica`, `hdt`, `cottas`, `qlever` | the same question as SPARQL | + +```bash +/opt/pycottas-venv/bin/python /opt/vcf-rdfizer/validation/regional_runner.py \ + --vcf /data/vcf/sample.vcf.gz --rdf /data/rdf/sample.nt.gz --rdf-format nt.gz \ + --arms cyvcf2-indexed,bcftools-indexed,qlever \ + --results-dir /data/regional --dataset-id slice +``` + +It writes `regional.csv` (one row per arm, question, window and replicate), +`regional.json` (medians, plus each side's one-time setup), `windows.json` (the +regions, so a rerun measures the same ones) and `mismatches.json`. + +Three things about it are easy to get wrong, so they are worth stating. + +**A window means POS, not overlap.** A tabix or CSI seek returns every record +whose *span* overlaps the region, so a deletion beginning before the window and +reaching into it comes back from htslib while a SPARQL `?pos` filter excludes +it. Every VCF arm therefore keeps the seek and then drops records whose POS +falls outside the window. The seek still does the work; the post-filter only +reconciles the convention, and without it the arms would disagree on any window +whose left edge cuts an indel -- which would look like a graph defect rather +than a coordinate convention. + +**Windows are anchored on real records.** A GIAB benchmark VCF covers a few +percent of the genome, so uniformly drawn 1 kb windows would be almost all +empty and the experiment would measure the cost of returning nothing. Each +window starts at a randomly chosen record position, from a fixed seed, so the +four sizes are a selectivity ladder rather than a sparsity ladder and every arm +sees identical regions. + +**The scan arm is timed on fewer windows than the others.** Its cost does not +depend on the window -- it reads the file whichever region is asked for -- so +timing it on every window measures one number repeatedly at roughly 11 s a go. +Correctness is still checked on every window: the equality reference is a single +whole-file pass that fills all of them at once, and it is deliberately not +timed, because a pass that answers eighty windows is not what a one-question +user pays for. + +Results are compared for exact equality before any timing is reported. The +runner exits non-zero when arms disagree, because a speed number from arms that +disagree is a bug report rather than a result. + ### Engine equivalence, and one place they differed Every engine is held to producing identical results. That is verified, not diff --git a/src/validation/queries/regional/common/r01_region_record_count.rq b/src/validation/queries/regional/common/r01_region_record_count.rq new file mode 100644 index 0000000..cb6149f --- /dev/null +++ b/src/validation/queries/regional/common/r01_region_record_count.rq @@ -0,0 +1,18 @@ +PREFIX vcfc: +PREFIX xsd: + +# R1 - how many records fall inside one coordinate window. +# +# The window is POS-based and 1-based inclusive, matching the VCF arms after +# their post-seek POS filter. A tabix seek returns every record whose span +# OVERLAPS the region, so an indel starting before the window would come back +# from htslib but not from this pattern; the VCF arms drop those explicitly so +# both sides mean the same thing. See regional_runner.REGION_SEMANTICS. +# +# {{CHROM}}, {{START}} and {{END}} are substituted by the runner. +SELECT (COUNT(DISTINCT ?record) AS ?recordCount) +WHERE { + ?record a vcfc:VCFRecord ; vcfc:chrom ?chrom ; vcfc:pos ?pos . + FILTER(STR(?chrom) = "{{CHROM}}") + FILTER(xsd:integer(?pos) >= {{START}} && xsd:integer(?pos) <= {{END}}) +} diff --git a/src/validation/queries/regional/common/r02_region_variant_shape_counts.rq b/src/validation/queries/regional/common/r02_region_variant_shape_counts.rq new file mode 100644 index 0000000..06781c4 --- /dev/null +++ b/src/validation/queries/regional/common/r02_region_variant_shape_counts.rq @@ -0,0 +1,41 @@ +PREFIX vcfc: +PREFIX xsd: + +# R2 - allele-shape distribution inside one window. The classification is +# character-for-character the whole-file q02 so the two are comparable; only the +# coordinate restriction is new. +SELECT ?variantClass (COUNT(DISTINCT ?record) AS ?recordCount) +WHERE { + ?record a vcfc:VCFRecord ; vcfc:chrom ?chrom ; vcfc:pos ?pos ; + vcfc:ref ?refLiteral ; vcfc:alt ?altLiteral . + FILTER(STR(?chrom) = "{{CHROM}}") + FILTER(xsd:integer(?pos) >= {{START}} && xsd:integer(?pos) <= {{END}}) + BIND(UCASE(STR(?refLiteral)) AS ?ref) + BIND(UCASE(STR(?altLiteral)) AS ?alt) + BIND( + IF(?alt = ".", "NO_ALT", + IF(CONTAINS(?alt, ","), "MULTIALLELIC", + IF( + ?alt = "*" || CONTAINS(?alt, "[") || CONTAINS(?alt, "]") || + (STRSTARTS(?alt, "<") && STRENDS(?alt, ">")), + "SYMBOLIC_OR_BREAKEND", + IF( + !REGEX(?ref, "^[ACGTN]+$") || !REGEX(?alt, "^[ACGTN]+$"), + "OTHER", + IF( + STRLEN(?ref) = 1 && STRLEN(?alt) = 1, + "SNV", + IF( + STRLEN(?ref) = STRLEN(?alt), + "MNV_OR_EQUAL_LENGTH_SUBSTITUTION", + IF(STRLEN(?ref) < STRLEN(?alt), "INSERTION_SHAPE", "DELETION_SHAPE") + ) + ) + ) + ) + ) + ) AS ?variantClass + ) +} +GROUP BY ?variantClass +ORDER BY ?variantClass diff --git a/src/validation/queries/regional/common/r03_region_titv.rq b/src/validation/queries/regional/common/r03_region_titv.rq new file mode 100644 index 0000000..c0921ce --- /dev/null +++ b/src/validation/queries/regional/common/r03_region_titv.rq @@ -0,0 +1,18 @@ +PREFIX vcfc: +PREFIX xsd: + +# R3 - the standard regional QC statistic: biallelic SNV transitions and +# transversions inside one window. +SELECT + (COUNT(DISTINCT ?record) AS ?biallelicSnvCount) + (SUM(IF((?ref = "A" && ?alt = "G") || (?ref = "G" && ?alt = "A") || (?ref = "C" && ?alt = "T") || (?ref = "T" && ?alt = "C"), 1, 0)) AS ?transitionCount) + (SUM(IF((?ref = "A" && ?alt = "G") || (?ref = "G" && ?alt = "A") || (?ref = "C" && ?alt = "T") || (?ref = "T" && ?alt = "C"), 0, 1)) AS ?transversionCount) +WHERE { + ?record a vcfc:VCFRecord ; vcfc:chrom ?chrom ; vcfc:pos ?pos ; + vcfc:ref ?refLiteral ; vcfc:alt ?altLiteral . + FILTER(STR(?chrom) = "{{CHROM}}") + FILTER(xsd:integer(?pos) >= {{START}} && xsd:integer(?pos) <= {{END}}) + BIND(UCASE(STR(?refLiteral)) AS ?ref) + BIND(UCASE(STR(?altLiteral)) AS ?alt) + FILTER(REGEX(?ref, "^[ACGT]$") && REGEX(?alt, "^[ACGT]$") && ?ref != ?alt) +} diff --git a/src/validation/queries/regional/common/r04_region_filter_distribution.rq b/src/validation/queries/regional/common/r04_region_filter_distribution.rq new file mode 100644 index 0000000..f26582a --- /dev/null +++ b/src/validation/queries/regional/common/r04_region_filter_distribution.rq @@ -0,0 +1,16 @@ +PREFIX vcfc: +PREFIX xsd: + +# R4 - FILTER status and exact lexical value inside one window. The VCF side +# uses `bcftools query -r`, which is the only route to the exact FILTER text. +SELECT ?filterStatus ?filterLexical (COUNT(DISTINCT ?record) AS ?recordCount) +WHERE { + ?record a vcfc:VCFRecord ; vcfc:chrom ?chrom ; vcfc:pos ?pos ; vcfc:hasCall ?call . + FILTER(STR(?chrom) = "{{CHROM}}") + FILTER(xsd:integer(?pos) >= {{START}} && xsd:integer(?pos) <= {{END}}) + ?call vcfc:filter ?filterLiteral . + BIND(STR(?filterLiteral) AS ?filterLexical) + BIND(IF(?filterLexical = "PASS", "PASS", IF(?filterLexical = ".", "NOT_APPLIED", "FAILED")) AS ?filterStatus) +} +GROUP BY ?filterStatus ?filterLexical +ORDER BY ?filterStatus ?filterLexical diff --git a/src/validation/queries/regional/expanded/r05_region_sample_genotype_counts.rq b/src/validation/queries/regional/expanded/r05_region_sample_genotype_counts.rq new file mode 100644 index 0000000..ed34824 --- /dev/null +++ b/src/validation/queries/regional/expanded/r05_region_sample_genotype_counts.rq @@ -0,0 +1,30 @@ +PREFIX vcfc: +PREFIX xsd: + +# R5 - per-sample genotype classes inside one window. +# +# Unlike the whole-file q05, this starts from the record rather than the +# SampleCall, because the coordinate restriction lives on the record. The +# genotype classification below is identical to q05's. +SELECT ?sampleId ?genotypeClass (COUNT(DISTINCT ?sampleCall) AS ?callCount) +WHERE { + ?record a vcfc:VCFRecord ; vcfc:chrom ?chrom ; vcfc:pos ?pos ; vcfc:hasCall ?call . + FILTER(STR(?chrom) = "{{CHROM}}") + FILTER(xsd:integer(?pos) >= {{START}} && xsd:integer(?pos) <= {{END}}) + ?call vcfc:hasSampleCall ?sampleCall . + ?sampleCall vcfc:sampleId ?sampleId . + OPTIONAL { + ?sampleCall vcfc:hasFormatValue ?gtValueNode . + FILTER(STRENDS(STR(?gtValueNode), "/fmt/GT")) + ?gtValueNode vcfc:fieldValue ?gtLiteral . + } + BIND(IF(BOUND(?gtLiteral), REPLACE(STR(?gtLiteral), "[|]", "/"), "") AS ?gt) + BIND(IF(!BOUND(?gtLiteral), "NO_GT_FIELD", + IF(CONTAINS(?gt, "."), "MISSING", + IF(REGEX(?gt, "^[0-9]+$"), IF(?gt = "0", "HAPLOID_REF", "HAPLOID_ALT"), + IF(REGEX(?gt, "^[0-9]+/[0-9]+$"), + IF(STRBEFORE(?gt, "/") = STRAFTER(?gt, "/"), IF(STRBEFORE(?gt, "/") = "0", "HOM_REF", "HOM_ALT"), "HET"), + "OTHER_PLOIDY")))) AS ?genotypeClass) +} +GROUP BY ?sampleId ?genotypeClass +ORDER BY ?sampleId ?genotypeClass diff --git a/src/validation/regional_runner.py b/src/validation/regional_runner.py new file mode 100644 index 0000000..f6c3a4e --- /dev/null +++ b/src/validation/regional_runner.py @@ -0,0 +1,1023 @@ +#!/usr/bin/env python3 +"""Indexed regional access: SPARQL against a coordinate-indexed VCF. + +The whole-file retrieval comparison in ``validation_runner`` pits a SPARQL +engine against cyvcf2 *scanning the file*, because none of Q1-Q13 is +coordinate-restricted. That is internally consistent, but it measures the graph +against the one VCF access mode nobody uses for a selective question. Real VCF +work seeks: bgzip + tabix, then ``bcftools -r`` or cyvcf2's region iterator, +which reads a handful of BGZF blocks instead of the file. + +This runner adds that arm. It asks the same five bioinformatic questions, +restricted to a coordinate window, of every access path: + + cyvcf2-scan whole-file iteration with a position filter (the status + quo, kept so the indexed arms have a baseline) + cyvcf2-indexed VCF(path)(region) over the bgzipped, indexed file + bcftools-indexed bcftools query -r + comunica|hdt|cottas|qlever the same question as SPARQL + +Results are compared for exact equality before any timing is reported, exactly +as the whole-file suite does. A speed number from arms that disagree is not a +result, it is a bug. + +Run it through ``benchmarks/14_regional_access.sh``, which mounts the inputs +and pins the image; the raw entry point is documented in ``--help``. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import random +import shutil +import statistics +import subprocess +import sys +import tempfile +import time +from collections import Counter +from pathlib import Path +from typing import Any + +SCRIPT_DIR = Path(__file__).resolve().parent +if str(SCRIPT_DIR) not in sys.path: + sys.path.insert(0, str(SCRIPT_DIR)) + +import validation_runner as V # noqa: E402 + +# --------------------------------------------------------------------------- +# What a window means, on both sides +# --------------------------------------------------------------------------- +REGION_SEMANTICS = """\ +A window is POS-based, 1-based and inclusive: a record belongs to it when +START <= POS <= END on the named contig. + +This has to be said explicitly because htslib does not mean that. A tabix or +CSI seek returns every record whose *span* overlaps the region, so a deletion +beginning before START and reaching into the window comes back from +`bcftools -r` and from cyvcf2's region iterator, while a POS filter in SPARQL +excludes it. Left alone, that difference would make the arms disagree on any +window whose left edge cuts an indel -- and it would look like a graph defect +rather than a coordinate convention. + +So both indexed arms keep the index seek (the thing being measured) and then +drop records whose POS falls outside the window. The seek still does the work; +the post-filter only reconciles the convention, costs one integer comparison +per returned record, and is applied by every VCF arm alike. +""" + +#: The regional questions, in report order. +REGIONAL_QUERIES = ( + "r01_region_record_count", + "r02_region_variant_shape_counts", + "r03_region_titv", + "r04_region_filter_distribution", + "r05_region_sample_genotype_counts", +) + +#: Same shape as ``validation_runner.QUERY_SCHEMAS``: fields, integer fields, +#: sort fields. Reused by the SPARQL normalizer so a regional answer is +#: canonicalised the same way a whole-file answer is. +REGIONAL_QUERY_SCHEMAS = { + "r01_region_record_count": (("recordCount",), {"recordCount"}, ()), + "r02_region_variant_shape_counts": (("variantClass", "recordCount"), {"recordCount"}, ("variantClass",)), + "r03_region_titv": ( + ("biallelicSnvCount", "transitionCount", "transversionCount"), + {"biallelicSnvCount", "transitionCount", "transversionCount"}, + (), + ), + "r04_region_filter_distribution": (("filterStatus", "filterLexical", "recordCount"), {"recordCount"}, ("filterStatus", "filterLexical")), + "r05_region_sample_genotype_counts": (("sampleId", "genotypeClass", "callCount"), {"callCount"}, ("sampleId", "genotypeClass")), +} + +#: Aggregates over an empty window still return one row, so these are shaped as +#: a single dict rather than a list. +REGIONAL_SINGLE_ROW = frozenset({"r01_region_record_count", "r03_region_titv"}) + +#: Queries that need the sample layer. On a single-sample file this is noise; +#: on a cohort it is most of the work, and it is the reason an arm's ranking can +#: differ between R1-R4 and R5. +REGIONAL_SAMPLE_LEVEL_QUERIES = frozenset({"r05_region_sample_genotype_counts"}) + +VCF_ARMS = ("cyvcf2-scan", "cyvcf2-indexed", "bcftools-indexed") +DEFAULT_WINDOW_SIZES = (1_000, 100_000, 1_000_000, 10_000_000) +DEFAULT_WINDOWS_PER_SIZE = 20 +DEFAULT_SEED = 20260923 + + +# --------------------------------------------------------------------------- +# Window selection +# --------------------------------------------------------------------------- +def contig_record_positions(vcf_path: Path) -> dict[str, list[int]]: + """One pass to learn where the records actually are. + + Windows are anchored on real record positions rather than drawn uniformly + across the contig. A GIAB benchmark VCF covers a few percent of the genome, + so uniform 1 kb windows would be almost all empty and the experiment would + measure the cost of returning nothing. Anchoring guarantees every window + holds at least one record, which is what makes the four window sizes a + selectivity ladder instead of a sparsity ladder. + + This is setup, not measurement: it runs once, before any timing. + """ + if V.VCF is None: + raise RuntimeError("cyvcf2 is unavailable; run this inside the VCF-RDFizer image") + positions: dict[str, list[int]] = {} + reader = V.VCF(str(vcf_path)) + try: + for variant in reader: + positions.setdefault(str(variant.CHROM), []).append(int(variant.POS)) + finally: + reader.close() + return positions + + +def draw_windows( + positions: dict[str, list[int]], + sizes: tuple[int, ...], + per_size: int, + seed: int, +) -> list[dict[str, Any]]: + """Draw the same windows for every arm, reproducibly. + + Each window is anchored at a randomly chosen record and extended rightwards, + so a window of any size contains that record. The draw is seeded, and the + resulting list is written to the results directory, so a rerun measures the + identical regions and a reader can check which ones they were. + """ + if not positions: + raise RuntimeError("the VCF has no records, so there is nothing to window") + rng = random.Random(seed) + contigs = sorted(positions) + windows: list[dict[str, Any]] = [] + for size in sizes: + for index in range(per_size): + contig = rng.choice(contigs) + anchor = rng.choice(positions[contig]) + # Anchor at the record, extend right. Starting at the anchor rather + # than centring on it keeps START >= 1 without a clamp that would + # quietly shrink the smallest windows. + start = anchor + end = start + size - 1 + windows.append({ + "window_id": f"w{size}_{index:02d}", + "chrom": contig, + "start": start, + "end": end, + "size": size, + "anchor_pos": anchor, + }) + return windows + + +def windows_expected_counts( + positions: dict[str, list[int]], windows: list[dict[str, Any]] +) -> None: + """Record how many records each window holds, for the report's x-axis. + + Selectivity is what the experiment varies, and window *size* is only a + proxy for it -- a 1 Mb window in a dense region can hold more records than + a 10 Mb window in a sparse one. Reporting the realised record count lets the + result be plotted against actual selectivity. + """ + for window in windows: + contig_positions = positions.get(window["chrom"], []) + start, end = window["start"], window["end"] + window["records_in_window"] = sum(1 for pos in contig_positions if start <= pos <= end) + + +# --------------------------------------------------------------------------- +# Index preparation (one-time cost, timed and reported separately) +# --------------------------------------------------------------------------- +def is_bgzf(path: Path) -> bool: + """True when the file is BGZF, which is what an index needs. + + A plain gzip .vcf.gz is not seekable, and tabix refuses it. Detecting this + up front turns a confusing downstream failure into one clear message. + """ + try: + with path.open("rb") as handle: + header = handle.read(18) + except OSError: + return False + if len(header) < 18 or header[:3] != b"\x1f\x8b\x08": + return False + # BGZF marks itself with an BC extra subfield in the gzip header. + return header[12:14] == b"BC" + + +def prepare_indexed_vcf( + vcf_path: Path, workdir: Path, *, index_kind: str = "auto" +) -> dict[str, Any]: + """Produce a bgzipped, coordinate-indexed copy, timing each step. + + Returns the paths and the one-time seconds. This is the VCF side's + equivalent of building an HDT or a QLever index, and it is reported the same + way: separately from query time, so neither side's setup is smuggled into a + per-question number. + + `index_kind`: "tbi" is the familiar tabix index but cannot address a contig + longer than 2^29-1 bp; "csi" has no such limit. "auto" picks csi when any + contig needs it, which is the only safe default for a whole-genome file. + """ + workdir.mkdir(parents=True, exist_ok=True) + report: dict[str, Any] = {"source": str(vcf_path), "indexKind": index_kind} + + target = workdir / (vcf_path.name if vcf_path.name.endswith(".gz") else vcf_path.name + ".gz") + compress_seconds = 0.0 + if is_bgzf(vcf_path): + # Already seekable: copy rather than recompress, and say so, because a + # recompression would inflate the one-time cost with work a real user + # would not repeat. + started = time.monotonic() + shutil.copy2(vcf_path, target) + compress_seconds = time.monotonic() - started + report["compression"] = "already-bgzf (copied)" + else: + started = time.monotonic() + _bgzip(vcf_path, target) + compress_seconds = time.monotonic() - started + report["compression"] = "bgzip" + report["bgzipSeconds"] = compress_seconds + report["bgzipBytes"] = target.stat().st_size + + if index_kind == "auto": + index_kind = "csi" if _needs_csi(target) else "tbi" + report["indexKind"] = index_kind + report["indexKindReason"] = ( + "a contig exceeds the 2^29-1 bp that a .tbi can address" + if index_kind == "csi" else "every contig fits a .tbi" + ) + + started = time.monotonic() + index_path = _index(target, index_kind) + report["indexSeconds"] = time.monotonic() - started + report["indexedVcf"] = str(target) + report["indexPath"] = str(index_path) + report["indexBytes"] = index_path.stat().st_size + report["totalSetupSeconds"] = report["bgzipSeconds"] + report["indexSeconds"] + return report + + +def _bgzip(source: Path, target: Path) -> None: + """Compress with bgzip, falling back to bcftools when bgzip is absent. + + The image installs `tabix`, which provides bgzip. The bcftools fallback + keeps the runner usable on an image built before that was added, and + produces the same BGZF container. + """ + if shutil.which("bgzip"): + with target.open("wb") as handle: + subprocess.run(["bgzip", "-c", str(source)], stdout=handle, check=True) + return + if shutil.which("bcftools"): + subprocess.run( + ["bcftools", "view", "--no-version", "-Oz", "-o", str(target), str(source)], + check=True, capture_output=True, + ) + return + raise RuntimeError("neither bgzip nor bcftools is available to produce a BGZF file") + + +def _needs_csi(bgzf_path: Path) -> bool: + """True when any contig is longer than a .tbi can address.""" + limit = 2 ** 29 - 1 + try: + header = V.read_vcf_header_text(bgzf_path) + except (OSError, ValueError): + return False + for line in header.splitlines(): + if not line.startswith("##contig="): + continue + for field in line[len("##contig=<"):].rstrip(">").split(","): + key, _, value = field.partition("=") + if key.strip() == "length": + try: + if int(value) > limit: + return True + except ValueError: + continue + return False + + +def _index(bgzf_path: Path, index_kind: str) -> Path: + suffix = ".csi" if index_kind == "csi" else ".tbi" + expected = Path(str(bgzf_path) + suffix) + if shutil.which("tabix"): + args = ["tabix", "-f", "-p", "vcf"] + if index_kind == "csi": + args.append("-C") + subprocess.run([*args, str(bgzf_path)], check=True, capture_output=True) + elif shutil.which("bcftools"): + flag = "-c" if index_kind == "csi" else "-t" + subprocess.run(["bcftools", "index", flag, "-f", str(bgzf_path)], + check=True, capture_output=True) + else: + raise RuntimeError("neither tabix nor bcftools is available to build an index") + if not expected.is_file(): + raise RuntimeError(f"index was not produced at {expected}") + return expected + + +# --------------------------------------------------------------------------- +# The VCF arms +# --------------------------------------------------------------------------- +# Every arm returns the same canonical structure for a given question, so the +# comparator is a plain equality test rather than a per-arm special case. The +# classification helpers are imported from validation_runner, not reimplemented, +# so a regional answer is classified character-for-character as the whole-file +# oracle classifies it. + +def _empty_answer(query_id: str) -> Any: + if query_id == "r01_region_record_count": + return {"recordCount": 0} + if query_id == "r03_region_titv": + return {"biallelicSnvCount": 0, "transitionCount": 0, "transversionCount": 0} + return [] + + +class _Accumulator: + """Per-question counters for one window.""" + + def __init__(self, query_id: str, samples: list[str]): + self.query_id = query_id + self.samples = samples + self.records = 0 + self.shapes: Counter[str] = Counter() + self.filters: Counter[tuple[str, str]] = Counter() + self.genotypes: Counter[tuple[str, str]] = Counter() + self.biallelic_snvs = 0 + self.transitions = 0 + self.transversions = 0 + + def add_variant(self, variant: Any) -> None: + query_id = self.query_id + if query_id == "r01_region_record_count": + self.records += 1 + return + + alt = V.alt_lexical(variant) + ref = str(variant.REF) + + if query_id == "r02_region_variant_shape_counts": + self.shapes[V.classify_variant_shape(ref, alt)] += 1 + return + + if query_id == "r03_region_titv": + ref_upper, alt_upper = ref.upper(), alt.upper() + if (len(ref_upper) == 1 and len(alt_upper) == 1 + and ref_upper in "ACGT" and alt_upper in "ACGT" + and ref_upper != alt_upper): + self.biallelic_snvs += 1 + if {ref_upper, alt_upper} in ({"A", "G"}, {"C", "T"}): + self.transitions += 1 + else: + self.transversions += 1 + return + + if query_id == "r04_region_filter_distribution": + raw = _filter_lexical(variant) + self.filters[(V.filter_status(raw), raw)] += 1 + return + + if query_id == "r05_region_sample_genotype_counts": + self.add_genotypes(variant) + + def add_genotypes(self, variant: Any) -> None: + has_gt = "GT" in (variant.FORMAT or []) + alleles: list[Any] = [None] * len(self.samples) + if self.samples and has_gt: + raw = list(variant.genotypes) + if len(raw) != len(self.samples): + raise ValueError( + f"sample/genotype length mismatch at {variant.CHROM}:{variant.POS}" + ) + alleles = [V.genotype_alleles(entry) for entry in raw] + for sample, call in zip(self.samples, alleles): + self.genotypes[(sample, V.classify_genotype(call, has_gt=has_gt))] += 1 + + def answer(self) -> Any: + query_id = self.query_id + if query_id == "r01_region_record_count": + return {"recordCount": self.records} + if query_id == "r02_region_variant_shape_counts": + return [{"variantClass": name, "recordCount": int(count)} + for name, count in sorted(self.shapes.items())] + if query_id == "r03_region_titv": + return {"biallelicSnvCount": self.biallelic_snvs, + "transitionCount": self.transitions, + "transversionCount": self.transversions} + if query_id == "r04_region_filter_distribution": + return [{"filterStatus": status, "filterLexical": lexical, + "recordCount": int(count)} + for (status, lexical), count in sorted(self.filters.items())] + return [{"sampleId": sample, "genotypeClass": genotype_class, + "callCount": int(count)} + for (sample, genotype_class), count in sorted(self.genotypes.items())] + + +def _filter_lexical(variant: Any) -> str: + """The FILTER column exactly as the file spells it. + + cyvcf2 reports PASS as None, so the mapping back to the lexical form has to + be explicit; `.` means no filter was applied and is a different state again. + """ + raw = variant.FILTER + if raw is None: + return "PASS" + return str(raw) + + +def vcf_sample_names(vcf_path: Path) -> list[str]: + header = V.read_vcf_header_text(vcf_path) + for line in header.splitlines(): + if line.startswith("#CHROM"): + fields = line.rstrip("\r\n").split("\t") + return fields[9:] if len(fields) > 9 else [] + return [] + + +def answer_cyvcf2( + vcf_path: Path, + query_id: str, + windows: list[dict[str, Any]], + samples: list[str], + *, + indexed: bool, +) -> dict[str, Any]: + """Answer one question for one or more windows with cyvcf2. + + ``indexed=True`` uses the region iterator, which seeks through the index; + ``indexed=False`` iterates the whole file and filters on POS, which is the + status quo. Both apply the POS filter, for the reason in REGION_SEMANTICS. + """ + if V.VCF is None: + raise RuntimeError("cyvcf2 is unavailable; run this inside the VCF-RDFizer image") + accumulators = {w["window_id"]: _Accumulator(query_id, samples) for w in windows} + + if indexed: + reader = V.VCF(str(vcf_path)) + try: + for window in windows: + region = f"{window['chrom']}:{window['start']}-{window['end']}" + accumulator = accumulators[window["window_id"]] + for variant in reader(region): + # htslib returns overlap; the window means POS. + if window["start"] <= int(variant.POS) <= window["end"]: + accumulator.add_variant(variant) + finally: + reader.close() + else: + by_contig: dict[str, list[dict[str, Any]]] = {} + for window in windows: + by_contig.setdefault(window["chrom"], []).append(window) + reader = V.VCF(str(vcf_path)) + try: + for variant in reader: + contig_windows = by_contig.get(str(variant.CHROM)) + if not contig_windows: + continue + pos = int(variant.POS) + for window in contig_windows: + if window["start"] <= pos <= window["end"]: + accumulators[window["window_id"]].add_variant(variant) + finally: + reader.close() + + return {window_id: acc.answer() for window_id, acc in accumulators.items()} + + +#: One bcftools format string per question: ask for the columns that question +#: needs and nothing else, so the arm is not penalised for fetching fields it +#: will discard. +#: +#: CHROM is in every one of them even though `-r` already restricts the contig. +#: Without it the folder would be trusting the region argument to have done its +#: job, and a record at the same coordinate on another contig would be counted +#: silently. One extra column makes the arm check itself. +BCFTOOLS_FORMATS = { + "r01_region_record_count": "%CHROM\t%POS\n", + "r02_region_variant_shape_counts": "%CHROM\t%POS\t%REF\t%ALT\n", + "r03_region_titv": "%CHROM\t%POS\t%REF\t%ALT\n", + "r04_region_filter_distribution": "%CHROM\t%POS\t%FILTER\n", + "r05_region_sample_genotype_counts": "%CHROM\t%POS[\t%GT]\n", +} + + +def answer_bcftools( + vcf_path: Path, query_id: str, window: dict[str, Any], samples: list[str] +) -> Any: + """Answer one question for one window with `bcftools query -r`. + + This is the arm a bioinformatician would actually reach for, and the only + one that reads the FILTER column's exact text without a Python parser in + the way. + """ + region = f"{window['chrom']}:{window['start']}-{window['end']}" + completed = subprocess.run( + ["bcftools", "query", "-r", region, "-f", BCFTOOLS_FORMATS[query_id], str(vcf_path)], + check=True, capture_output=True, text=True, + ) + return _fold_bcftools_output(completed.stdout, query_id, window, samples) + + +def _fold_bcftools_output( + output: str, query_id: str, window: dict[str, Any], samples: list[str] +) -> Any: + records = 0 + shapes: Counter[str] = Counter() + filters: Counter[tuple[str, str]] = Counter() + genotypes: Counter[tuple[str, str]] = Counter() + biallelic = transitions = transversions = 0 + + for line in output.splitlines(): + if not line: + continue + fields = line.split("\t") + chrom, pos = fields[0], int(fields[1]) + # The contig is checked rather than assumed, and POS is reconciled the + # same way the cyvcf2 arms reconcile it (see REGION_SEMANTICS). + if chrom != str(window["chrom"]): + continue + if not (window["start"] <= pos <= window["end"]): + continue + + if query_id == "r01_region_record_count": + records += 1 + elif query_id == "r02_region_variant_shape_counts": + shapes[V.classify_variant_shape(fields[2], fields[3])] += 1 + elif query_id == "r03_region_titv": + ref, alt = fields[2].upper(), fields[3].upper() + if (len(ref) == 1 and len(alt) == 1 and ref in "ACGT" + and alt in "ACGT" and ref != alt): + biallelic += 1 + if {ref, alt} in ({"A", "G"}, {"C", "T"}): + transitions += 1 + else: + transversions += 1 + elif query_id == "r04_region_filter_distribution": + raw = fields[2] + filters[(V.filter_status(raw), raw)] += 1 + elif query_id == "r05_region_sample_genotype_counts": + calls = fields[2:] + for sample, raw_gt in zip(samples, calls): + genotypes[(sample, _classify_bcftools_gt(raw_gt))] += 1 + + if query_id == "r01_region_record_count": + return {"recordCount": records} + if query_id == "r02_region_variant_shape_counts": + return [{"variantClass": name, "recordCount": int(count)} + for name, count in sorted(shapes.items())] + if query_id == "r03_region_titv": + return {"biallelicSnvCount": biallelic, "transitionCount": transitions, + "transversionCount": transversions} + if query_id == "r04_region_filter_distribution": + return [{"filterStatus": status, "filterLexical": lexical, "recordCount": int(count)} + for (status, lexical), count in sorted(filters.items())] + return [{"sampleId": sample, "genotypeClass": genotype_class, "callCount": int(count)} + for (sample, genotype_class), count in sorted(genotypes.items())] + + +def _classify_bcftools_gt(raw_gt: str) -> str: + """Map a `%GT` string onto the same classes cyvcf2's oracle produces. + + `bcftools query` hands back the genotype as text, so this reproduces + `classify_genotype` from the lexical form rather than from allele indices. + A record with no GT field prints `.` here, which is indistinguishable from + a missing call -- so a file whose records lack GT entirely would make this + arm disagree with the cyvcf2 arms, and the comparator would catch it. + """ + normalized = raw_gt.replace("|", "/") + if normalized in ("", "."): + return "MISSING" + if "." in normalized: + return "MISSING" + parts = normalized.split("/") + if len(parts) == 1: + return "HAPLOID_REF" if parts[0] == "0" else "HAPLOID_ALT" + if len(parts) == 2: + if parts[0] == parts[1]: + return "HOM_REF" if parts[0] == "0" else "HOM_ALT" + return "HET" + return "OTHER_PLOIDY" + + +# --------------------------------------------------------------------------- +# The SPARQL arm +# --------------------------------------------------------------------------- +REGIONAL_QUERY_ROOT = SCRIPT_DIR / "queries" / "regional" + + +def regional_query_path(representation: str, query_id: str) -> Path: + """Regional templates live beside the whole-file ones, same lookup order.""" + candidate = REGIONAL_QUERY_ROOT / representation / f"{query_id}.rq" + if candidate.is_file(): + return candidate + return REGIONAL_QUERY_ROOT / "common" / f"{query_id}.rq" + + +def render_query(template_path: Path, window: dict[str, Any]) -> str: + """Substitute one window into a template. + + The contig name is substituted into a SPARQL string literal, so a name + containing a quote or a backslash would change the query's meaning. VCF + contig names cannot contain whitespace and in practice are alphanumeric, + but this refuses rather than silently emitting a broken -- or differently + scoped -- query. + """ + chrom = str(window["chrom"]) + if any(character in chrom for character in '"\\\n\r\t'): + raise ValueError(f"contig name is not safe to substitute into SPARQL: {chrom!r}") + return (template_path.read_text(encoding="utf-8") + .replace("{{CHROM}}", chrom) + .replace("{{START}}", str(int(window["start"]))) + .replace("{{END}}", str(int(window["end"])))) + + +def normalize_regional(query_id: str, raw_path: Path) -> Any: + """Canonicalise SPARQL Results JSON into the arms' shared shape. + + One deliberate canonicalisation: for the single-row aggregates, a field the + engine leaves unbound is read as 0. An empty window makes `SUM` unbound on + some engines and 0 on others, which is an engine difference rather than a + disagreement about the data -- and without this, every empty window would + be reported as a mismatch between engines that in fact agree. + """ + fields, integer_fields, sort_fields = REGIONAL_QUERY_SCHEMAS[query_id] + rows: list[dict[str, Any]] = [] + for number, binding in enumerate(V.bindings(raw_path), start=1): + row: dict[str, Any] = {} + for field in fields: + entry = binding.get(field) + if entry is None: + if field in integer_fields and query_id in REGIONAL_SINGLE_ROW: + row[field] = 0 + continue + raise ValueError(f"row {number} of {query_id} has no {field!r}") + lexical = entry["value"] + row[field] = int(lexical) if field in integer_fields else str(lexical) + rows.append(row) + + if query_id in REGIONAL_SINGLE_ROW: + if not rows: + return _empty_answer(query_id) + if len(rows) != 1: + raise ValueError(f"{query_id} must return one row, got {len(rows)}") + return rows[0] + + if sort_fields: + rows.sort(key=lambda row: tuple(row[field] for field in sort_fields)) + return rows + + +# --------------------------------------------------------------------------- +# Comparison +# --------------------------------------------------------------------------- +def canonical(answer: Any) -> str: + """A stable string for exact comparison between arms.""" + return json.dumps(answer, sort_keys=True, separators=(",", ":")) + + +# --------------------------------------------------------------------------- +# Driver +# --------------------------------------------------------------------------- +def timed_vcf_arm( + arm: str, + vcf_path: Path, + indexed_vcf: Path, + query_id: str, + window: dict[str, Any], + samples: list[str], +) -> tuple[Any, float]: + """Run one question, one window, one arm, and time exactly that.""" + started = time.monotonic() + if arm == "bcftools-indexed": + answer = answer_bcftools(indexed_vcf, query_id, window, samples) + elif arm == "cyvcf2-indexed": + answer = answer_cyvcf2(indexed_vcf, query_id, [window], samples, indexed=True)[window["window_id"]] + elif arm == "cyvcf2-scan": + answer = answer_cyvcf2(vcf_path, query_id, [window], samples, indexed=False)[window["window_id"]] + else: + raise ValueError(f"not a VCF arm: {arm}") + return answer, time.monotonic() - started + + +def run(args: argparse.Namespace) -> int: + results_dir = args.results_dir + results_dir.mkdir(parents=True, exist_ok=True) + raw_dir = results_dir / "raw" + raw_dir.mkdir(parents=True, exist_ok=True) + + arms = [arm.strip() for arm in args.arms.split(",") if arm.strip()] + unknown = [a for a in arms if a not in VCF_ARMS and a not in V.SPARQL_ENGINES] + if unknown: + raise SystemExit(f"unknown arm(s): {', '.join(unknown)}") + queries = [q.strip() for q in args.queries.split(",") if q.strip()] + unknown_queries = [q for q in queries if q not in REGIONAL_QUERIES] + if unknown_queries: + raise SystemExit(f"unknown question(s): {', '.join(unknown_queries)}") + + samples = vcf_sample_names(args.vcf) + sizes = tuple(int(s) for s in args.window_sizes.split(",")) + + print(f"[setup] scanning {args.vcf.name} to place windows on real records") + setup_started = time.monotonic() + positions = contig_record_positions(args.vcf) + position_scan_seconds = time.monotonic() - setup_started + total_records = sum(len(v) for v in positions.values()) + print(f"[setup] {total_records:,} records across {len(positions)} contigs " + f"in {position_scan_seconds:.1f}s") + + windows = draw_windows(positions, sizes, args.windows_per_size, args.seed) + windows_expected_counts(positions, windows) + V.write_json(results_dir / "windows.json", { + "seed": args.seed, + "sizes": list(sizes), + "windowsPerSize": args.windows_per_size, + "anchoring": "each window starts at a randomly chosen record position", + "regionSemantics": REGION_SEMANTICS, + "positionScanSeconds": position_scan_seconds, + "totalRecords": total_records, + "windows": windows, + }) + + setup: dict[str, Any] = {"positionScanSeconds": position_scan_seconds} + + needs_index = any(a in ("cyvcf2-indexed", "bcftools-indexed") for a in arms) + indexed_vcf = args.vcf + if needs_index: + print("[setup] building the bgzip + index (one-time VCF-side cost)") + index_report = prepare_indexed_vcf( + args.vcf, args.scratch_dir / "indexed", index_kind=args.index_kind + ) + indexed_vcf = Path(index_report["indexedVcf"]) + setup["vcfIndex"] = index_report + print(f"[setup] bgzip {index_report['bgzipSeconds']:.1f}s + " + f"{index_report['indexKind']} index {index_report['indexSeconds']:.1f}s") + + # Reference answers: one whole-file pass per question covering EVERY window + # at once. This is the equality reference, and it is deliberately not timed + # -- a pass that fills 80 windows is not what a one-question user pays. + print("[reference] computing expected answers for every window") + reference: dict[str, dict[str, Any]] = {} + for query_id in queries: + reference[query_id] = answer_cyvcf2(args.vcf, query_id, windows, samples, indexed=False) + + rows: list[dict[str, Any]] = [] + mismatches: list[dict[str, Any]] = [] + + def record(arm: str, query_id: str, window: dict[str, Any], replicate: int, + answer: Any, seconds: float, status: str = "OK", + error: str | None = None) -> None: + expected = reference[query_id][window["window_id"]] + agrees = canonical(answer) == canonical(expected) if status == "OK" else None + if agrees is False: + mismatches.append({ + "arm": arm, "query_id": query_id, "window_id": window["window_id"], + "expected": expected, "observed": answer, + }) + rows.append({ + "arm": arm, + "query_id": query_id, + "window_id": window["window_id"], + "chrom": window["chrom"], + "start": window["start"], + "end": window["end"], + "window_size": window["size"], + "records_in_window": window["records_in_window"], + "replicate": replicate, + "wall_seconds": seconds, + "status": status, + "agrees_with_reference": agrees, + "error": error, + }) + + # The scan arm's cost does not depend on the window -- it reads the whole + # file whichever region is asked for. Timing it on all 80 windows would + # measure one number 80 times, at roughly 11 s a go. It is sampled instead, + # and the sample size is reported so the thinness is visible. + scan_windows = _sample_scan_windows(windows, args.scan_windows_per_size) + + for arm in [a for a in arms if a in VCF_ARMS]: + targets = scan_windows if arm == "cyvcf2-scan" else windows + print(f"[arm] {arm}: {len(targets)} windows x {len(queries)} questions " + f"x {args.replicates} replicates") + for query_id in queries: + for window in targets: + for replicate in range(1, args.replicates + 1): + try: + answer, seconds = timed_vcf_arm( + arm, args.vcf, indexed_vcf, query_id, window, samples) + record(arm, query_id, window, replicate, answer, seconds) + except (subprocess.CalledProcessError, OSError, ValueError) as error: + record(arm, query_id, window, replicate, None, float("nan"), + status="FAILED", error=str(error)) + + engine_arms = [a for a in arms if a in V.SPARQL_ENGINES] + if engine_arms and args.rdf is None: + raise SystemExit("--rdf is required when a SPARQL engine is among --arms") + + for arm in engine_arms: + print(f"[arm] {arm}: preparing engine") + engine_raw = raw_dir / arm + engine_raw.mkdir(parents=True, exist_ok=True) + options = { + "query_timeout": args.query_timeout, + "rdf_format": args.rdf_format, + "qlever_memory_gb": args.qlever_memory_gb, + "qlever_port": args.qlever_port, + "qlever_startup_timeout": args.qlever_startup_timeout, + } + try: + engine = V.build_engine(arm, args.rdf, raw_dir=engine_raw, + scratch=args.scratch_dir, options=options) + except (ValueError, RuntimeError) as error: + setup.setdefault("engineErrors", {})[arm] = str(error) + print(f"[arm] {arm}: unavailable ({error})") + continue + + try: + with engine: + setup.setdefault("engineSetupSeconds", {})[arm] = engine.setup_seconds + print(f"[arm] {arm}: setup {engine.setup_seconds:.2f}s; " + f"{len(windows)} windows x {len(queries)} questions " + f"x {args.replicates} replicates") + with tempfile.TemporaryDirectory(dir=str(args.scratch_dir)) as rendered_dir: + rendered_root = Path(rendered_dir) + for query_id in queries: + template = regional_query_path(args.representation, query_id) + for window in windows: + rendered = rendered_root / f"{query_id}__{window['window_id']}.rq" + rendered.write_text(render_query(template, window), encoding="utf-8") + for replicate in range(1, args.replicates + 1): + envelope = engine.execute( + f"{query_id}__{window['window_id']}__r{replicate}", rendered) + if envelope["status"] != "PASS": + record(arm, query_id, window, replicate, None, + envelope["wallSeconds"], status="FAILED", + error=envelope.get("error") or "engine execution failed") + continue + try: + answer = normalize_regional( + query_id, Path(envelope["rawResult"])) + except (ValueError, KeyError, OSError) as error: + record(arm, query_id, window, replicate, None, + envelope["wallSeconds"], status="UNREADABLE", + error=str(error)) + continue + record(arm, query_id, window, replicate, answer, + envelope["wallSeconds"]) + except RuntimeError as error: + setup.setdefault("engineErrors", {})[arm] = str(error) + print(f"[arm] {arm}: failed ({error})") + + _write_outputs(results_dir, args, rows, mismatches, setup, windows, queries, arms) + + failures = sum(1 for row in rows if row["status"] != "OK") + disagreements = sum(1 for row in rows if row["agrees_with_reference"] is False) + print(f"\n{len(rows)} timed executions, {failures} failed, " + f"{disagreements} disagreed with the reference") + if disagreements: + print("Arms disagree, so no speed comparison from this run is reportable.") + print(f"See {results_dir / 'mismatches.json'}") + return 1 + return 0 + + +def _sample_scan_windows( + windows: list[dict[str, Any]], per_size: int +) -> list[dict[str, Any]]: + """The first N windows of each size, in the order they were drawn.""" + seen: Counter[int] = Counter() + sampled = [] + for window in windows: + if seen[window["size"]] < per_size: + sampled.append(window) + seen[window["size"]] += 1 + return sampled + + +# --------------------------------------------------------------------------- +# Output +# --------------------------------------------------------------------------- +REGIONAL_CSV_HEADER = [ + "arm", "query_id", "window_id", "chrom", "start", "end", + "window_size", "records_in_window", "replicate", + "wall_seconds", "status", "agrees_with_reference", "error", +] + + +def _write_outputs( + results_dir: Path, + args: argparse.Namespace, + rows: list[dict[str, Any]], + mismatches: list[dict[str, Any]], + setup: dict[str, Any], + windows: list[dict[str, Any]], + queries: list[str], + arms: list[str], +) -> None: + csv_path = results_dir / "regional.csv" + with csv_path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.DictWriter(handle, fieldnames=REGIONAL_CSV_HEADER) + writer.writeheader() + for row in rows: + writer.writerow({key: row.get(key) for key in REGIONAL_CSV_HEADER}) + + # Medians per (arm, question, window size): replicate noise collapsed, but + # the raw rows stay in the CSV so the collapse can be checked. + summary: dict[str, Any] = {} + grouped: dict[tuple[str, str, int], list[float]] = {} + for row in rows: + if row["status"] != "OK": + continue + grouped.setdefault( + (row["arm"], row["query_id"], row["window_size"]), [] + ).append(row["wall_seconds"]) + for (arm, query_id, size), seconds in sorted(grouped.items()): + summary.setdefault(arm, {}).setdefault(query_id, {})[str(size)] = { + "medianSeconds": statistics.median(seconds), + "minSeconds": min(seconds), + "maxSeconds": max(seconds), + "executions": len(seconds), + } + + V.write_json(results_dir / "regional.json", { + "datasetId": args.dataset_id, + "vcf": str(args.vcf), + "vcfSha256": V.sha256_file(args.vcf), + "rdf": str(args.rdf) if args.rdf else None, + "representation": args.representation, + "arms": arms, + "queries": queries, + "replicates": args.replicates, + "scanWindowsPerSize": args.scan_windows_per_size, + "windowCount": len(windows), + "setup": setup, + "regionSemantics": REGION_SEMANTICS, + "summary": summary, + "executions": len(rows), + "failures": sum(1 for row in rows if row["status"] != "OK"), + "disagreements": sum(1 for row in rows if row["agrees_with_reference"] is False), + }) + + V.write_json(results_dir / "mismatches.json", { + "count": len(mismatches), + "note": "Any entry here invalidates the speed comparison for that arm.", + "mismatches": mismatches[:200], + }) + + print(f"\nwrote {csv_path}") + print(f"wrote {results_dir / 'regional.json'}") + print(f"wrote {results_dir / 'windows.json'}") + + +def build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + parser.add_argument("--vcf", type=Path, required=True, + help="source VCF; the scan arm and the reference read this") + parser.add_argument("--rdf", type=Path, default=None, + help="RDF artifact for the SPARQL arms (.nt/.nt.gz/.hdt/.cottas)") + parser.add_argument("--rdf-format", default="nt", + help="format of --rdf, as validation_runner names it") + parser.add_argument("--representation", choices=("expanded", "condensed"), + default="expanded") + parser.add_argument("--arms", default=",".join(VCF_ARMS) + ",qlever", + help="comma-separated: " + ", ".join(VCF_ARMS + V.SPARQL_ENGINES)) + parser.add_argument("--queries", default=",".join(REGIONAL_QUERIES)) + parser.add_argument("--window-sizes", + default=",".join(str(s) for s in DEFAULT_WINDOW_SIZES)) + parser.add_argument("--windows-per-size", type=int, default=DEFAULT_WINDOWS_PER_SIZE) + parser.add_argument("--scan-windows-per-size", type=int, default=3, + help="windows timed for cyvcf2-scan, whose cost does not " + "depend on the window (default: 3)") + parser.add_argument("--replicates", type=int, default=3) + parser.add_argument("--seed", type=int, default=DEFAULT_SEED) + parser.add_argument("--index-kind", choices=("auto", "tbi", "csi"), default="auto") + parser.add_argument("--results-dir", type=Path, required=True) + parser.add_argument("--dataset-id", required=True) + parser.add_argument("--scratch-dir", type=Path, default=Path("/work")) + parser.add_argument("--query-timeout", type=int, default=V.DEFAULT_QUERY_TIMEOUT) + parser.add_argument("--qlever-memory-gb", type=int, default=4) + parser.add_argument("--qlever-port", type=int, default=7019) + parser.add_argument("--qlever-startup-timeout", type=int, default=900) + return parser + + +def main(argv: list[str] | None = None) -> int: + args = build_parser().parse_args(argv) + args.scratch_dir.mkdir(parents=True, exist_ok=True) + try: + return run(args) + except (RuntimeError, ValueError, OSError) as error: + print(f"error: {error}", file=sys.stderr) + return 2 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/test/test_regional_runner_unit.py b/test/test_regional_runner_unit.py new file mode 100644 index 0000000..64872a8 --- /dev/null +++ b/test/test_regional_runner_unit.py @@ -0,0 +1,315 @@ +"""Unit tests for the indexed regional-access arm. + +The point of the regional experiment is a speed comparison between access +paths, and a speed comparison is only meaningful once the paths agree on the +answer. So most of what is worth testing here is agreement: that the bcftools +arm folds its text output into the same structure the cyvcf2 arm builds from +parsed records, that the SPARQL normalizer produces that structure too, and +that the window convention is applied identically everywhere. + +The bcftools and index-building paths need binaries that are present in the +VCF-RDFizer image but not necessarily on a developer machine, so those tests +skip rather than fail when the binary is absent. The folding logic itself is a +pure function and is always tested. +""" + +from __future__ import annotations + +import json +import shutil +import sys +import unittest +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "src" / "validation")) + +from test.helpers import VerboseTestCase # noqa: E402 + +try: + import regional_runner as R +except ImportError: # pragma: no cover - the module must import to test it + R = None + +FIXTURES = Path(__file__).resolve().parent / "test_vcf_files" +SMALL_VCF = FIXTURES / "test-1k.vcf" + +cyvcf2_available = R is not None and R.V.VCF is not None + + +def _bcftools_text(vcf_path: Path, query_id: str, samples: list[str]) -> str: + """Reproduce `bcftools query -f ` output from the VCF text. + + This is what makes the cross-check meaningful: the bcftools arm is fed the + exact columns bcftools would have printed, without needing bcftools on the + machine running the test. + """ + lines = [] + for line in vcf_path.read_text().splitlines(): + if line.startswith("#"): + continue + fields = line.split("\t") + chrom, pos, ref, alt, filt = fields[0], fields[1], fields[3], fields[4], fields[6] + if query_id == "r01_region_record_count": + lines.append(f"{chrom}\t{pos}") + elif query_id in ("r02_region_variant_shape_counts", "r03_region_titv"): + lines.append(f"{chrom}\t{pos}\t{ref}\t{alt}") + elif query_id == "r04_region_filter_distribution": + lines.append(f"{chrom}\t{pos}\t{filt}") + elif query_id == "r05_region_sample_genotype_counts": + gts = [] + format_keys = fields[8].split(":") if len(fields) > 8 else [] + gt_index = format_keys.index("GT") if "GT" in format_keys else None + for payload in fields[9:9 + len(samples)]: + parts = payload.split(":") + gts.append(parts[gt_index] if gt_index is not None + and gt_index < len(parts) else ".") + lines.append("\t".join([chrom, pos, *gts])) + return "\n".join(lines) + "\n" + + +@unittest.skipUnless(cyvcf2_available, "cyvcf2 is required") +class ArmAgreementTests(VerboseTestCase): + """The two VCF arms must produce byte-identical canonical answers.""" + + @classmethod + def setUpClass(cls): + cls.positions = R.contig_record_positions(SMALL_VCF) + cls.samples = R.vcf_sample_names(SMALL_VCF) + cls.windows = R.draw_windows(cls.positions, (1_000, 100_000), 3, R.DEFAULT_SEED) + R.windows_expected_counts(cls.positions, cls.windows) + + def test_bcftools_folding_agrees_with_the_cyvcf2_arm(self): + for query_id in R.REGIONAL_QUERIES: + text = _bcftools_text(SMALL_VCF, query_id, self.samples) + expected = R.answer_cyvcf2( + SMALL_VCF, query_id, self.windows, self.samples, indexed=False) + for window in self.windows: + with self.subTest(query=query_id, window=window["window_id"]): + folded = R._fold_bcftools_output( + text, query_id, window, self.samples) + self.assertEqual( + R.canonical(folded), + R.canonical(expected[window["window_id"]]), + ) + + def test_every_window_holds_at_least_one_record(self): + """Anchoring is what makes the sizes a selectivity ladder.""" + for window in self.windows: + with self.subTest(window=window["window_id"]): + self.assertGreater(window["records_in_window"], 0) + + def test_the_record_count_question_matches_the_windows_own_count(self): + answers = R.answer_cyvcf2( + SMALL_VCF, "r01_region_record_count", self.windows, self.samples, indexed=False) + for window in self.windows: + self.assertEqual(answers[window["window_id"]]["recordCount"], + window["records_in_window"]) + + def test_windows_are_reproducible_from_the_seed(self): + """Every arm has to see the same regions, across runs and machines.""" + again = R.draw_windows(self.positions, (1_000, 100_000), 3, R.DEFAULT_SEED) + self.assertEqual([w["window_id"] for w in again], + [w["window_id"] for w in self.windows]) + self.assertEqual([(w["chrom"], w["start"]) for w in again], + [(w["chrom"], w["start"]) for w in self.windows]) + + def test_a_different_seed_draws_different_windows(self): + other = R.draw_windows(self.positions, (1_000,), 3, R.DEFAULT_SEED + 1) + mine = [w for w in self.windows if w["size"] == 1_000] + self.assertNotEqual([(w["chrom"], w["start"]) for w in other], + [(w["chrom"], w["start"]) for w in mine]) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class RegionConventionTests(VerboseTestCase): + """A window means POS, on every arm.""" + + WINDOW = {"window_id": "w", "chrom": "1", "start": 100, "end": 200, "size": 101} + + def test_a_record_outside_the_window_is_dropped_by_the_folder(self): + """htslib returns overlap; the post-filter reconciles it to POS.""" + text = "\n".join(f"1\t{pos}" for pos in (50, 100, 150, 200, 201)) + "\n" + folded = R._fold_bcftools_output(text, "r01_region_record_count", self.WINDOW, []) + self.assertEqual(folded, {"recordCount": 3}) + + def test_the_window_is_inclusive_at_both_edges(self): + text = "1\t100\n1\t200\n" + folded = R._fold_bcftools_output(text, "r01_region_record_count", self.WINDOW, []) + self.assertEqual(folded, {"recordCount": 2}) + + def test_a_record_on_another_contig_is_not_counted(self): + """Regression: the folder used to trust `-r` and filter on POS alone. + + A record at the same coordinate on a different contig was then counted, + which `-r` happens to prevent -- so the arm was correct only by relying + on an argument it never checked. Every format string now carries CHROM. + """ + text = "\n".join(["1\t150", "2\t150", "X\t150"]) + "\n" + folded = R._fold_bcftools_output(text, "r01_region_record_count", self.WINDOW, []) + self.assertEqual(folded, {"recordCount": 1}) + + def test_every_bcftools_format_asks_for_the_contig(self): + for query_id, fmt in R.BCFTOOLS_FORMATS.items(): + with self.subTest(query=query_id): + self.assertTrue(fmt.startswith("%CHROM\t%POS"), fmt) + + def test_region_semantics_are_documented_next_to_the_code(self): + """The convention is the easiest thing here to get silently wrong.""" + self.assertIn("POS", R.REGION_SEMANTICS) + self.assertIn("overlap", R.REGION_SEMANTICS.lower()) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class GenotypeClassTests(VerboseTestCase): + """`%GT` text must land in the same classes the allele-index oracle uses.""" + + def test_the_lexical_classifier_matches_the_allele_classifier(self): + cases = { + "0/0": "HOM_REF", "0|0": "HOM_REF", + "1/1": "HOM_ALT", "2|2": "HOM_ALT", + "0/1": "HET", "1|2": "HET", + "0": "HAPLOID_REF", "1": "HAPLOID_ALT", + "./.": "MISSING", ".": "MISSING", "./1": "MISSING", "": "MISSING", + "0/1/1": "OTHER_PLOIDY", + } + for raw, expected in cases.items(): + with self.subTest(gt=raw): + self.assertEqual(R._classify_bcftools_gt(raw), expected) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class NormalizerTests(VerboseTestCase): + """SPARQL Results JSON folds into the same shape the VCF arms produce.""" + + def _raw(self, tmp: Path, bindings: list[dict]) -> Path: + path = tmp / "raw.json" + path.write_text(json.dumps({"results": {"bindings": bindings}})) + return path + + def test_a_single_row_aggregate_becomes_a_dict(self): + import tempfile + with tempfile.TemporaryDirectory() as td: + raw = self._raw(Path(td), [{"recordCount": {"value": "42"}}]) + self.assertEqual( + R.normalize_regional("r01_region_record_count", raw), + {"recordCount": 42}, + ) + + def test_an_unbound_sum_on_an_empty_window_reads_as_zero(self): + """Engines disagree on SUM over nothing; that is not a data mismatch.""" + import tempfile + with tempfile.TemporaryDirectory() as td: + raw = self._raw(Path(td), [{"biallelicSnvCount": {"value": "0"}}]) + self.assertEqual( + R.normalize_regional("r03_region_titv", raw), + {"biallelicSnvCount": 0, "transitionCount": 0, "transversionCount": 0}, + ) + + def test_no_rows_at_all_is_still_a_well_formed_empty_answer(self): + import tempfile + with tempfile.TemporaryDirectory() as td: + raw = self._raw(Path(td), []) + self.assertEqual( + R.normalize_regional("r01_region_record_count", raw), + {"recordCount": 0}, + ) + self.assertEqual( + R.normalize_regional("r02_region_variant_shape_counts", raw), []) + + def test_grouped_rows_are_sorted_into_the_canonical_order(self): + import tempfile + with tempfile.TemporaryDirectory() as td: + raw = self._raw(Path(td), [ + {"variantClass": {"value": "SNV"}, "recordCount": {"value": "2"}}, + {"variantClass": {"value": "DELETION_SHAPE"}, "recordCount": {"value": "1"}}, + ]) + self.assertEqual( + R.normalize_regional("r02_region_variant_shape_counts", raw), + [{"variantClass": "DELETION_SHAPE", "recordCount": 1}, + {"variantClass": "SNV", "recordCount": 2}], + ) + + def test_a_missing_field_on_a_grouped_row_is_an_error_not_a_zero(self): + """Only the single-row aggregates get the unbound-means-zero rule.""" + import tempfile + with tempfile.TemporaryDirectory() as td: + raw = self._raw(Path(td), [{"variantClass": {"value": "SNV"}}]) + with self.assertRaises(ValueError): + R.normalize_regional("r02_region_variant_shape_counts", raw) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class QueryTemplateTests(VerboseTestCase): + def test_every_question_has_a_template(self): + for query_id in R.REGIONAL_QUERIES: + with self.subTest(query=query_id): + self.assertTrue(R.regional_query_path("expanded", query_id).is_file()) + + def test_rendering_substitutes_all_three_placeholders(self): + window = {"chrom": "chr7", "start": 100, "end": 200} + for query_id in R.REGIONAL_QUERIES: + with self.subTest(query=query_id): + text = R.render_query(R.regional_query_path("expanded", query_id), window) + self.assertNotIn("{{", text) + self.assertIn('"chr7"', text) + self.assertIn("100", text) + self.assertIn("200", text) + + def test_a_contig_name_that_would_break_out_of_the_literal_is_refused(self): + """Substitution into a SPARQL string is the one injection risk here.""" + for bad in ('1" || (1=1) || "', "a\\b", "a\nb"): + with self.subTest(chrom=bad): + with self.assertRaises(ValueError): + R.render_query( + R.regional_query_path("expanded", "r01_region_record_count"), + {"chrom": bad, "start": 1, "end": 2}, + ) + + def test_every_template_restricts_on_both_contig_and_position(self): + """A template missing a bound would silently answer a different question.""" + for query_id in R.REGIONAL_QUERIES: + with self.subTest(query=query_id): + text = R.regional_query_path("expanded", query_id).read_text() + self.assertIn("{{CHROM}}", text) + self.assertIn("{{START}}", text) + self.assertIn("{{END}}", text) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class ScanSamplingTests(VerboseTestCase): + def test_the_scan_arm_is_sampled_evenly_across_window_sizes(self): + """Its cost does not vary with the window, so it is measured thinly.""" + windows = [{"window_id": f"w{size}_{i:02d}", "size": size} + for size in (1_000, 100_000) for i in range(20)] + sampled = R._sample_scan_windows(windows, 3) + self.assertEqual(len(sampled), 6) + self.assertEqual(sum(1 for w in sampled if w["size"] == 1_000), 3) + self.assertEqual(sum(1 for w in sampled if w["size"] == 100_000), 3) + + +@unittest.skipUnless( + R is not None and (shutil.which("bgzip") or shutil.which("bcftools")), + "bgzip or bcftools is required to build a BGZF file", +) +class IndexPreparationTests(VerboseTestCase): + def test_building_the_index_reports_its_one_time_cost(self): + import tempfile + with tempfile.TemporaryDirectory() as td: + report = R.prepare_indexed_vcf(SMALL_VCF, Path(td), index_kind="tbi") + self.assertGreaterEqual(report["bgzipSeconds"], 0.0) + self.assertGreaterEqual(report["indexSeconds"], 0.0) + self.assertTrue(Path(report["indexPath"]).is_file()) + self.assertTrue(R.is_bgzf(Path(report["indexedVcf"]))) + + def test_a_plain_gzip_file_is_not_mistaken_for_bgzf(self): + """tabix cannot seek a plain gzip stream; detecting it early is kinder.""" + import gzip + import tempfile + with tempfile.TemporaryDirectory() as td: + plain = Path(td) / "plain.vcf.gz" + plain.write_bytes(gzip.compress(SMALL_VCF.read_bytes())) + self.assertFalse(R.is_bgzf(plain)) + + +if __name__ == "__main__": + unittest.main() From eba43714df4b9a72848df7c044a6f5e7f867fcf7 Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 20:03:41 +0200 Subject: [PATCH 2/9] Test the regional runner's driver, index preparation and failure paths The existing regional tests cover agreement between the arms' folding and classification logic, which left the code around it untested: building the bgzip + index copy, dispatching a timed execution to the right arm, and the driver itself -- 49% of regional_runner.py. test_regional_runner_driver_unit.py adds 37 tests: - BGZF detection, including a real bgzip file and a plain gzip that tabix would refuse. - Index kind: a contig past 2^29-1 bp needs CSI; an unparseable length or an unreadable header falls back to TBI. - The tool fallbacks: clear errors when neither bgzip/bcftools or tabix/bcftools exists, an index that a tool claims to have written but did not, and bcftools building both .tbi and .csi when tabix is absent. - prepare_indexed_vcf built for real: plain input is bgzipped, BGZF input is copied rather than recompressed, auto picks CSI for a long contig. - The indexed arms against a real index: cyvcf2's seek and bcftools query -r agree with the scan on every question and window. - The driver end to end: a scan-only run writes all four outputs and agrees; a disagreeing arm exits 1 and is written to mismatches.json; a failing arm is recorded without aborting; unknown arms and questions, a SPARQL arm without a graph, and a record-less VCF are refused cleanly; all three VCF arms agree. - The engine plumbing through a stand-in engine: agreement on every window, per-execution failure, unreadable output, and an engine that cannot start. The binary-dependent tests skip where bgzip, tabix or bcftools are absent and run in full inside the VCF-RDFizer image. Co-Authored-By: Claude Opus 5.5 --- test/test_regional_runner_driver_unit.py | 527 +++++++++++++++++++++++ 1 file changed, 527 insertions(+) create mode 100644 test/test_regional_runner_driver_unit.py diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver_unit.py new file mode 100644 index 0000000..53833fa --- /dev/null +++ b/test/test_regional_runner_driver_unit.py @@ -0,0 +1,527 @@ +"""The regional runner's driver, index preparation and failure paths. + +test_regional_runner_unit.py covers agreement between the arms' folding and +classification logic. This module covers what surrounds it: building the +bgzip + index copy the indexed arms seek through, dispatching one timed +execution to the right arm, and the driver end to end -- the files it writes, +the exit status it returns when arms disagree, and how a failing arm or engine +is recorded rather than aborting the run. + +The index and bcftools tests use the real binaries. They are present in the +VCF-RDFizer image, which is where these tests are meant to run in full; on a +machine without them those tests skip rather than pass vacuously. The SPARQL +engine plumbing is driven through a stand-in engine, because what is under test +there is the runner's handling of an engine's envelope, not the engine. +""" + +from __future__ import annotations + +import csv +import gzip +import json +import shutil +import sys +import tempfile +import unittest +from contextlib import redirect_stdout +from io import StringIO +from pathlib import Path +from unittest import mock + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "src" / "validation")) + +from test.helpers import VerboseTestCase # noqa: E402 + +try: + import regional_runner as R +except ImportError: # pragma: no cover - the module must import to test it + R = None + +FIXTURES = Path(__file__).resolve().parent / "test_vcf_files" +SMALL_VCF = FIXTURES / "test-1k.vcf" + +cyvcf2_available = R is not None and R.V.VCF is not None +have_bgzip = shutil.which("bgzip") is not None +have_tabix = shutil.which("tabix") is not None +have_bcftools = shutil.which("bcftools") is not None + +HEADER = ( + "##fileformat=VCFv4.2\n" + "##contig=\n" + "##FORMAT=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1\n" +) + + +def write_vcf(path: Path, *, length: int = 1_000_000, records: int = 3) -> Path: + lines = [HEADER.format(length=length)] + for index in range(records): + lines.append(f"1\t{100 + index * 10}\t.\tA\tG\t50\tPASS\t.\tGT\t0/1\n") + path.write_text("".join(lines), encoding="utf-8") + return path + + +def run_main(argv: list[str]) -> tuple[int, str]: + out = StringIO() + with redirect_stdout(out): + status = R.main(argv) + return status, out.getvalue() + + +# --------------------------------------------------------------------------- +# Index preparation +# --------------------------------------------------------------------------- +@unittest.skipUnless(R is not None, "regional_runner must import") +class BgzfDetectionTests(VerboseTestCase): + def test_a_missing_file_is_not_bgzf(self): + self.assertFalse(R.is_bgzf(Path("/nonexistent/input.vcf.gz"))) + + def test_a_file_shorter_than_a_gzip_header_is_not_bgzf(self): + with tempfile.TemporaryDirectory() as td: + path = Path(td) / "short.gz" + path.write_bytes(b"\x1f\x8b\x08") + self.assertFalse(R.is_bgzf(path)) + + def test_plain_text_is_not_bgzf(self): + self.assertFalse(R.is_bgzf(SMALL_VCF)) + + def test_plain_gzip_is_not_bgzf(self): + """A plain .vcf.gz cannot be seeked, and tabix refuses it.""" + with tempfile.TemporaryDirectory() as td: + path = Path(td) / "plain.vcf.gz" + with gzip.open(path, "wb") as handle: + handle.write(SMALL_VCF.read_bytes()) + self.assertFalse(R.is_bgzf(path)) + + @unittest.skipUnless(have_bgzip, "bgzip is required to produce a real BGZF file") + def test_real_bgzip_output_is_bgzf(self): + with tempfile.TemporaryDirectory() as td: + target = Path(td) / "small.vcf.gz" + R._bgzip(SMALL_VCF, target) + self.assertTrue(R.is_bgzf(target)) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class IndexKindTests(VerboseTestCase): + """A .tbi cannot address a contig longer than 2^29-1 bp; a .csi can.""" + + def test_a_contig_within_the_tbi_limit_does_not_need_csi(self): + with tempfile.TemporaryDirectory() as td: + path = write_vcf(Path(td) / "short.vcf", length=2 ** 29 - 1) + self.assertFalse(R._needs_csi(path)) + + def test_a_contig_beyond_the_tbi_limit_needs_csi(self): + with tempfile.TemporaryDirectory() as td: + path = write_vcf(Path(td) / "long.vcf", length=2 ** 29) + self.assertTrue(R._needs_csi(path)) + + def test_an_unparseable_length_is_ignored_rather_than_fatal(self): + with tempfile.TemporaryDirectory() as td: + path = Path(td) / "odd.vcf" + path.write_text(HEADER.replace("length={length}", "length=unknown"), + encoding="utf-8") + self.assertFalse(R._needs_csi(path)) + + def test_an_unreadable_header_falls_back_to_tbi(self): + self.assertFalse(R._needs_csi(Path("/nonexistent/input.vcf.gz"))) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class ToolFallbackTests(VerboseTestCase): + """What the runner does when a binary it prefers is missing.""" + + def test_bgzip_without_bgzip_or_bcftools_is_a_clear_error(self): + with tempfile.TemporaryDirectory() as td, \ + mock.patch.object(R.shutil, "which", return_value=None): + with self.assertRaisesRegex(RuntimeError, "neither bgzip nor bcftools"): + R._bgzip(SMALL_VCF, Path(td) / "out.vcf.gz") + + def test_indexing_without_tabix_or_bcftools_is_a_clear_error(self): + with tempfile.TemporaryDirectory() as td, \ + mock.patch.object(R.shutil, "which", return_value=None): + with self.assertRaisesRegex(RuntimeError, "neither tabix nor bcftools"): + R._index(Path(td) / "in.vcf.gz", "tbi") + + def test_an_index_that_was_not_written_is_reported(self): + """A tool can exit 0 and still leave nothing behind; that must not pass.""" + with tempfile.TemporaryDirectory() as td, \ + mock.patch.object(R.shutil, "which", return_value="/usr/bin/tabix"), \ + mock.patch.object(R.subprocess, "run") as run: + with self.assertRaisesRegex(RuntimeError, "index was not produced"): + R._index(Path(td) / "in.vcf.gz", "csi") + self.assertIn("-C", run.call_args.args[0]) + + @unittest.skipUnless(have_bcftools and have_bgzip, "bcftools and bgzip are required") + def test_bcftools_builds_the_index_when_tabix_is_absent(self): + """The fallback keeps the runner usable on an image without tabix.""" + real_which = shutil.which + with tempfile.TemporaryDirectory() as td: + target = Path(td) / "small.vcf.gz" + R._bgzip(SMALL_VCF, target) + with mock.patch.object( + R.shutil, "which", + side_effect=lambda name: None if name == "tabix" else real_which(name)): + for kind, suffix in (("tbi", ".tbi"), ("csi", ".csi")): + with self.subTest(kind=kind): + index = R._index(target, kind) + self.assertEqual(index, Path(str(target) + suffix)) + self.assertGreater(index.stat().st_size, 0) + + +@unittest.skipUnless(R is not None and have_bgzip and have_tabix, + "bgzip and tabix are required") +class PrepareIndexedVcfTests(VerboseTestCase): + """The VCF side's one-time cost, built for real and reported separately.""" + + def test_a_plain_vcf_is_bgzipped_and_indexed(self): + with tempfile.TemporaryDirectory() as td: + report = R.prepare_indexed_vcf(SMALL_VCF, Path(td) / "indexed") + self.assertEqual(report["compression"], "bgzip") + self.assertEqual(report["indexKind"], "tbi") + self.assertEqual(report["indexKindReason"], "every contig fits a .tbi") + self.assertTrue(R.is_bgzf(Path(report["indexedVcf"]))) + self.assertTrue(Path(report["indexPath"]).is_file()) + self.assertGreater(report["indexBytes"], 0) + self.assertAlmostEqual(report["totalSetupSeconds"], + report["bgzipSeconds"] + report["indexSeconds"]) + + def test_an_already_bgzf_input_is_copied_not_recompressed(self): + """Recompressing would charge the VCF side for work no user repeats.""" + with tempfile.TemporaryDirectory() as td: + source = Path(td) / "source.vcf.gz" + R._bgzip(SMALL_VCF, source) + report = R.prepare_indexed_vcf(source, Path(td) / "indexed") + self.assertEqual(report["compression"], "already-bgzf (copied)") + self.assertEqual(Path(report["indexedVcf"]).read_bytes(), source.read_bytes()) + + def test_auto_picks_csi_for_a_contig_a_tbi_cannot_address(self): + with tempfile.TemporaryDirectory() as td: + source = write_vcf(Path(td) / "long.vcf", length=2 ** 29 + 10) + report = R.prepare_indexed_vcf(source, Path(td) / "indexed") + self.assertEqual(report["indexKind"], "csi") + self.assertTrue(report["indexPath"].endswith(".csi")) + self.assertIn("2^29-1", report["indexKindReason"]) + + def test_an_explicit_index_kind_is_honoured(self): + with tempfile.TemporaryDirectory() as td: + report = R.prepare_indexed_vcf(SMALL_VCF, Path(td) / "indexed", index_kind="csi") + self.assertEqual(report["indexKind"], "csi") + self.assertNotIn("indexKindReason", report) + + +# --------------------------------------------------------------------------- +# The indexed arms against a real index +# --------------------------------------------------------------------------- +@unittest.skipUnless(cyvcf2_available and have_bgzip and have_tabix, + "cyvcf2, bgzip and tabix are required") +class IndexedArmTests(VerboseTestCase): + """The seek path must return exactly what the scan path returns.""" + + @classmethod + def setUpClass(cls): + cls._tmp = tempfile.TemporaryDirectory() + report = R.prepare_indexed_vcf(SMALL_VCF, Path(cls._tmp.name) / "indexed") + cls.indexed = Path(report["indexedVcf"]) + cls.samples = R.vcf_sample_names(SMALL_VCF) + positions = R.contig_record_positions(SMALL_VCF) + cls.windows = R.draw_windows(positions, (1_000, 100_000, 1_000_000), 3, R.DEFAULT_SEED) + R.windows_expected_counts(positions, cls.windows) + + @classmethod + def tearDownClass(cls): + cls._tmp.cleanup() + + def test_the_cyvcf2_seek_agrees_with_the_scan_on_every_question(self): + for query_id in R.REGIONAL_QUERIES: + scan = R.answer_cyvcf2(SMALL_VCF, query_id, self.windows, self.samples, indexed=False) + seek = R.answer_cyvcf2(self.indexed, query_id, self.windows, self.samples, indexed=True) + for window in self.windows: + with self.subTest(query=query_id, window=window["window_id"]): + self.assertEqual(R.canonical(seek[window["window_id"]]), + R.canonical(scan[window["window_id"]])) + + @unittest.skipUnless(have_bcftools, "bcftools is required") + def test_bcftools_query_agrees_with_the_scan_on_every_question(self): + for query_id in R.REGIONAL_QUERIES: + scan = R.answer_cyvcf2(SMALL_VCF, query_id, self.windows, self.samples, indexed=False) + for window in self.windows: + with self.subTest(query=query_id, window=window["window_id"]): + answer = R.answer_bcftools(self.indexed, query_id, window, self.samples) + self.assertEqual(R.canonical(answer), + R.canonical(scan[window["window_id"]])) + + +# --------------------------------------------------------------------------- +# Dispatch and small helpers +# --------------------------------------------------------------------------- +@unittest.skipUnless(R is not None, "regional_runner must import") +class DispatchTests(VerboseTestCase): + WINDOW = {"window_id": "w1", "chrom": "1", "start": 1, "end": 100, "size": 100} + + def test_an_unknown_arm_is_refused(self): + with self.assertRaisesRegex(ValueError, "not a VCF arm"): + R.timed_vcf_arm("qlever", SMALL_VCF, SMALL_VCF, "r01_region_record_count", + self.WINDOW, []) + + def test_each_arm_reads_the_file_it_is_meant_to(self): + """The scan reads the source; both indexed arms read the indexed copy.""" + source, indexed = Path("source.vcf"), Path("indexed.vcf.gz") + with mock.patch.object(R, "answer_cyvcf2", return_value={"w1": "a"}) as cy, \ + mock.patch.object(R, "answer_bcftools", return_value="b") as bc: + self.assertEqual(R.timed_vcf_arm("cyvcf2-scan", source, indexed, + "r01_region_record_count", self.WINDOW, [])[0], "a") + self.assertEqual(cy.call_args.args[0], source) + self.assertFalse(cy.call_args.kwargs["indexed"]) + R.timed_vcf_arm("cyvcf2-indexed", source, indexed, + "r01_region_record_count", self.WINDOW, []) + self.assertEqual(cy.call_args.args[0], indexed) + self.assertTrue(cy.call_args.kwargs["indexed"]) + self.assertEqual(R.timed_vcf_arm("bcftools-indexed", source, indexed, + "r01_region_record_count", self.WINDOW, [])[0], "b") + self.assertEqual(bc.call_args.args[0], indexed) + + def test_the_timing_covers_only_the_execution(self): + with mock.patch.object(R, "answer_bcftools", return_value={}): + _, seconds = R.timed_vcf_arm("bcftools-indexed", SMALL_VCF, SMALL_VCF, + "r01_region_record_count", self.WINDOW, []) + self.assertGreaterEqual(seconds, 0.0) + + def test_an_empty_window_has_a_well_formed_answer_for_every_question(self): + self.assertEqual(R._empty_answer("r01_region_record_count"), {"recordCount": 0}) + self.assertEqual(R._empty_answer("r03_region_titv"), + {"biallelicSnvCount": 0, "transitionCount": 0, "transversionCount": 0}) + for query_id in ("r02_region_variant_shape_counts", "r04_region_filter_distribution", + "r05_region_sample_genotype_counts"): + self.assertEqual(R._empty_answer(query_id), []) + + def test_a_sites_only_header_has_no_samples(self): + with tempfile.TemporaryDirectory() as td: + path = Path(td) / "sites.vcf" + path.write_text("##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n", + encoding="utf-8") + self.assertEqual(R.vcf_sample_names(path), []) + + def test_a_genotype_row_of_the_wrong_width_is_an_error(self): + """A sample/genotype mismatch must not be silently truncated by zip().""" + accumulator = R._Accumulator("r05_region_sample_genotype_counts", ["A", "B"]) + variant = mock.Mock(FORMAT=["GT"], genotypes=[[0, 1, False]], CHROM="1", POS=5) + with self.assertRaisesRegex(ValueError, "length mismatch"): + accumulator.add_genotypes(variant) + + def test_a_missing_filter_reads_back_as_pass(self): + """cyvcf2 reports PASS as None; the lexical form has to be restored.""" + self.assertEqual(R._filter_lexical(mock.Mock(FILTER=None)), "PASS") + self.assertEqual(R._filter_lexical(mock.Mock(FILTER="q10;s50")), "q10;s50") + + +# --------------------------------------------------------------------------- +# The driver end to end +# --------------------------------------------------------------------------- +def base_argv(results: Path, scratch: Path, arms: str, **extra) -> list[str]: + argv = [ + "--vcf", str(SMALL_VCF), + "--arms", arms, + "--window-sizes", "1000,100000", + "--windows-per-size", "2", + "--scan-windows-per-size", "1", + "--replicates", "1", + "--results-dir", str(results), + "--scratch-dir", str(scratch), + "--dataset-id", "test-1k", + ] + for key, value in extra.items(): + argv += [f"--{key.replace('_', '-')}", str(value)] + return argv + + +@unittest.skipUnless(cyvcf2_available, "cyvcf2 is required") +class DriverTests(VerboseTestCase): + def setUp(self): + self._tmp = tempfile.TemporaryDirectory() + self.root = Path(self._tmp.name) + self.results = self.root / "results" + self.scratch = self.root / "scratch" + + def tearDown(self): + self._tmp.cleanup() + + def read(self, name: str): + return json.loads((self.results / name).read_text(encoding="utf-8")) + + def test_a_scan_only_run_writes_every_output_and_agrees(self): + status, _ = run_main(base_argv(self.results, self.scratch, "cyvcf2-scan")) + self.assertEqual(status, 0) + for name in ("regional.csv", "regional.json", "windows.json", "mismatches.json"): + self.assertTrue((self.results / name).is_file(), name) + + report = self.read("regional.json") + # One scan window per size, five questions, one replicate. + self.assertEqual(report["executions"], 2 * 5) + self.assertEqual(report["failures"], 0) + self.assertEqual(report["disagreements"], 0) + self.assertEqual(report["windowCount"], 4) + self.assertEqual(set(report["summary"]["cyvcf2-scan"]), set(R.REGIONAL_QUERIES)) + self.assertNotIn("vcfIndex", report["setup"]) # no indexed arm, no index built + + windows = self.read("windows.json") + self.assertEqual(windows["seed"], R.DEFAULT_SEED) + self.assertEqual(len(windows["windows"]), 4) + self.assertEqual(self.read("mismatches.json")["count"], 0) + + with (self.results / "regional.csv").open(encoding="utf-8") as handle: + rows = list(csv.DictReader(handle)) + self.assertEqual(len(rows), 10) + self.assertTrue(all(row["agrees_with_reference"] == "True" for row in rows)) + + def test_a_disagreeing_arm_fails_the_run_and_is_written_down(self): + """A speed from arms that disagree is a bug, not a result.""" + with mock.patch.object(R, "timed_vcf_arm", return_value=({"recordCount": -1}, 0.01)): + status, output = run_main(base_argv(self.results, self.scratch, "cyvcf2-scan", + queries="r01_region_record_count")) + self.assertEqual(status, 1) + self.assertIn("no speed comparison", output) + mismatches = self.read("mismatches.json") + self.assertEqual(mismatches["count"], 2) + self.assertEqual(mismatches["mismatches"][0]["observed"], {"recordCount": -1}) + + def test_a_failing_arm_is_recorded_without_aborting_the_run(self): + with mock.patch.object(R, "timed_vcf_arm", side_effect=OSError("disk vanished")): + status, _ = run_main(base_argv(self.results, self.scratch, "cyvcf2-scan", + queries="r01_region_record_count")) + # Failures are not disagreements: nothing was compared, so exit 0. + self.assertEqual(status, 0) + report = self.read("regional.json") + self.assertEqual(report["failures"], 2) + self.assertEqual(report["disagreements"], 0) + with (self.results / "regional.csv").open(encoding="utf-8") as handle: + rows = list(csv.DictReader(handle)) + self.assertEqual({row["status"] for row in rows}, {"FAILED"}) + self.assertEqual({row["error"] for row in rows}, {"disk vanished"}) + + def test_an_unknown_arm_or_question_is_refused_before_any_work(self): + with self.assertRaises(SystemExit): + run_main(base_argv(self.results, self.scratch, "cyvcf2-scan,duckdb")) + with self.assertRaises(SystemExit): + run_main(base_argv(self.results, self.scratch, "cyvcf2-scan", queries="r99_nothing")) + + def test_a_sparql_arm_without_a_graph_is_refused(self): + with self.assertRaisesRegex(SystemExit, "--rdf is required"): + run_main(base_argv(self.results, self.scratch, "qlever")) + + def test_a_vcf_without_records_is_an_error_exit_not_a_traceback(self): + path = self.root / "empty.vcf" + path.write_text("##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n", + encoding="utf-8") + argv = base_argv(self.results, self.scratch, "cyvcf2-scan") + argv[argv.index("--vcf") + 1] = str(path) + from contextlib import redirect_stderr + err = StringIO() + with redirect_stderr(err): + status, _ = run_main(argv) + self.assertEqual(status, 2) + self.assertIn("no records", err.getvalue()) + + @unittest.skipUnless(have_bgzip and have_tabix and have_bcftools, + "bgzip, tabix and bcftools are required") + def test_all_three_vcf_arms_agree_end_to_end(self): + status, _ = run_main(base_argv( + self.results, self.scratch, "cyvcf2-scan,cyvcf2-indexed,bcftools-indexed")) + self.assertEqual(status, 0) + report = self.read("regional.json") + self.assertEqual(report["disagreements"], 0) + self.assertEqual(report["failures"], 0) + self.assertIn("vcfIndex", report["setup"]) + # scan: 2 windows; the indexed arms: all 4 windows; five questions each. + self.assertEqual(report["executions"], (2 + 4 + 4) * 5) + self.assertEqual(set(report["summary"]), + {"cyvcf2-scan", "cyvcf2-indexed", "bcftools-indexed"}) + + +class _StandInEngine: + """Answers each rendered query from the VCF, so agreement is testable.""" + + def __init__(self, results: Path, mode: str): + self.results = results + self.mode = mode + self.setup_seconds = 0.25 + + def __enter__(self): + return self + + def __exit__(self, *exc): + return False + + def execute(self, name: str, query_path: Path) -> dict: + query_id, window_id, _ = name.split("__") + if self.mode == "fail": + return {"status": "FAILED", "wallSeconds": 0.01, "error": "engine said no"} + raw = query_path.with_suffix(".json") + if self.mode == "garbage": + raw.write_text("{}", encoding="utf-8") + else: + windows = json.loads((self.results / "windows.json").read_text())["windows"] + window = next(w for w in windows if w["window_id"] == window_id) + answer = R.answer_cyvcf2(SMALL_VCF, query_id, [window], [], indexed=False)[window_id] + raw.write_text(json.dumps({"head": {"vars": list(answer)}, "results": { + "bindings": [{k: {"type": "literal", "value": str(v)} for k, v in answer.items()}] + }}), encoding="utf-8") + return {"status": "PASS", "wallSeconds": 0.02, "rawResult": str(raw)} + + +@unittest.skipUnless(cyvcf2_available, "cyvcf2 is required") +class EngineArmTests(VerboseTestCase): + """How the driver handles an engine: agreement, failure, unreadable output.""" + + QUERIES = "r01_region_record_count,r03_region_titv" + + def setUp(self): + self._tmp = tempfile.TemporaryDirectory() + self.root = Path(self._tmp.name) + self.results = self.root / "results" + self.scratch = self.root / "scratch" + self.graph = self.root / "graph.nt" + self.graph.write_text("", encoding="utf-8") + + def tearDown(self): + self._tmp.cleanup() + + def run_with(self, mode: str, build=None): + build = build or (lambda *a, **k: _StandInEngine(self.results, mode)) + with mock.patch.object(R.V, "build_engine", side_effect=build): + status, _ = run_main(base_argv(self.results, self.scratch, "qlever", + queries=self.QUERIES, rdf=self.graph)) + report = json.loads((self.results / "regional.json").read_text()) + with (self.results / "regional.csv").open(encoding="utf-8") as handle: + rows = list(csv.DictReader(handle)) + return status, report, rows + + def test_an_agreeing_engine_is_timed_on_every_window(self): + status, report, rows = self.run_with("agree") + self.assertEqual(status, 0) + self.assertEqual(report["setup"]["engineSetupSeconds"], {"qlever": 0.25}) + self.assertEqual(len(rows), 4 * 2) # every window, two questions + self.assertTrue(all(row["agrees_with_reference"] == "True" for row in rows)) + + def test_an_engine_failure_is_recorded_per_execution(self): + status, report, rows = self.run_with("fail") + self.assertEqual(status, 0) + self.assertEqual(report["failures"], 8) + self.assertEqual({row["error"] for row in rows}, {"engine said no"}) + + def test_an_unreadable_result_is_not_mistaken_for_an_answer(self): + status, report, rows = self.run_with("garbage") + self.assertEqual({row["status"] for row in rows}, {"UNREADABLE"}) + self.assertEqual(report["disagreements"], 0) + + def test_an_engine_that_cannot_start_is_reported_and_skipped(self): + def refuse(*args, **kwargs): + raise RuntimeError("qlever-index not on PATH") + + status, report, rows = self.run_with("agree", build=refuse) + self.assertEqual(status, 0) + self.assertEqual(report["setup"]["engineErrors"], {"qlever": "qlever-index not on PATH"}) + self.assertEqual(rows, []) + + +if __name__ == "__main__": + unittest.main() From c3815a487d43fe40526633bcd5ebf06e6aedec0f Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 20:10:29 +0200 Subject: [PATCH 3/9] Sort an unsorted VCF before indexing it for the regional arms Found by running the regional tests inside the image on vcf-bench-1, where bgzip, tabix and bcftools exist: every attempt to index test-1k failed with "[E::hts_idx_push] Chromosome blocks not continuous", from tabix and from the bcftools fallback alike. test-1k interleaves its two contigs record by record, and an index needs coordinate order; a VCF is not obliged to be in it. The existing index test would have caught this, but it skips off a machine without those binaries and so had never run. prepare_indexed_vcf now recognises that failure, coordinate-sorts with bcftools sort, and indexes the sorted copy. The sort is what a user would have to do before seeking, so it is charged to the VCF side's one-time cost -- as its own sortSeconds, not folded into the index time -- and the report says why it happened. Only the sorted copy is kept. Any other index failure is raised as before rather than sorted around, and a missing bcftools is a clear error. The benchmark inputs (the HG005 slices) are already sorted and are unaffected. Tests: the unsorted-input decision is covered without binaries; the sort, the already-sorted path and the fallbacks run against the real tools, now on a sorted synthetic VCF where the test is about something other than sorting. Co-Authored-By: Claude Opus 5.5 --- src/validation/regional_runner.py | 54 +++++++++++++++++- test/test_regional_runner_driver_unit.py | 70 ++++++++++++++++++++++-- 2 files changed, 116 insertions(+), 8 deletions(-) diff --git a/src/validation/regional_runner.py b/src/validation/regional_runner.py index f6c3a4e..3515413 100644 --- a/src/validation/regional_runner.py +++ b/src/validation/regional_runner.py @@ -251,16 +251,66 @@ def prepare_indexed_vcf( if index_kind == "csi" else "every contig fits a .tbi" ) + report["sortSeconds"] = 0.0 started = time.monotonic() - index_path = _index(target, index_kind) + try: + index_path = _index(target, index_kind) + except subprocess.CalledProcessError as error: + if not _is_unsorted_error(error): + raise + # An index needs coordinate order, and a VCF is not required to be in + # it. Sorting is what a user would have to do before seeking, so it is + # done here and charged to the VCF side's one-time cost -- reported as + # its own line rather than folded into the index time. + sorted_target = target.with_name(target.name.replace(".vcf.gz", "") + ".sorted.vcf.gz") + sort_started = time.monotonic() + _sort(vcf_path, sorted_target, workdir) + report["sortSeconds"] = time.monotonic() - sort_started + report["sorted"] = "input was not coordinate-sorted; sorted with bcftools sort" + target.unlink(missing_ok=True) + target = sorted_target + report["bgzipBytes"] = target.stat().st_size + started = time.monotonic() + index_path = _index(target, index_kind) report["indexSeconds"] = time.monotonic() - started report["indexedVcf"] = str(target) report["indexPath"] = str(index_path) report["indexBytes"] = index_path.stat().st_size - report["totalSetupSeconds"] = report["bgzipSeconds"] + report["indexSeconds"] + report["totalSetupSeconds"] = ( + report["bgzipSeconds"] + report["sortSeconds"] + report["indexSeconds"] + ) return report +#: What tabix and bcftools print when records are not in coordinate order: +#: "Chromosome blocks not continuous" when a contig recurs after another one, +#: "unsorted positions" when POS goes backwards within a contig. +UNSORTED_MARKERS = ("not continuous", "unsorted") + + +def _is_unsorted_error(error: subprocess.CalledProcessError) -> bool: + stderr = error.stderr or b"" + if isinstance(stderr, bytes): + stderr = stderr.decode("utf-8", "replace") + return any(marker in stderr.lower() for marker in UNSORTED_MARKERS) + + +def _sort(source: Path, target: Path, workdir: Path) -> None: + """Coordinate-sort into a BGZF file, the only way an unsorted VCF can be indexed.""" + if not shutil.which("bcftools"): + raise RuntimeError( + "the VCF is not coordinate-sorted, so it cannot be indexed, and bcftools " + "is not available to sort it" + ) + sort_tmp = workdir / "bcftools-sort-tmp" + sort_tmp.mkdir(parents=True, exist_ok=True) + subprocess.run( + ["bcftools", "sort", "-Oz", "-T", str(sort_tmp), "-o", str(target), str(source)], + check=True, capture_output=True, + ) + shutil.rmtree(sort_tmp, ignore_errors=True) + + def _bgzip(source: Path, target: Path) -> None: """Compress with bgzip, falling back to bcftools when bgzip is absent. diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver_unit.py index 53833fa..4cee0ef 100644 --- a/test/test_regional_runner_driver_unit.py +++ b/test/test_regional_runner_driver_unit.py @@ -157,7 +157,7 @@ def test_bcftools_builds_the_index_when_tabix_is_absent(self): real_which = shutil.which with tempfile.TemporaryDirectory() as td: target = Path(td) / "small.vcf.gz" - R._bgzip(SMALL_VCF, target) + R._bgzip(write_vcf(Path(td) / "small.vcf"), target) with mock.patch.object( R.shutil, "which", side_effect=lambda name: None if name == "tabix" else real_which(name)): @@ -168,11 +168,66 @@ def test_bcftools_builds_the_index_when_tabix_is_absent(self): self.assertGreater(index.stat().st_size, 0) +@unittest.skipUnless(R is not None, "regional_runner must import") +class UnsortedInputTests(VerboseTestCase): + """An index needs coordinate order; a VCF is not obliged to have it.""" + + def error(self, stderr): + return R.subprocess.CalledProcessError(1, ["tabix"], stderr=stderr) + + def test_both_tools_unsorted_messages_are_recognised(self): + self.assertTrue(R._is_unsorted_error(self.error(b"[E::hts_idx_push] Chromosome blocks not continuous"))) + self.assertTrue(R._is_unsorted_error(self.error("[E::hts_idx_push] Unsorted positions on sequence #1"))) + + def test_any_other_index_failure_is_not_read_as_unsorted(self): + self.assertFalse(R._is_unsorted_error(self.error(b"[E::hts_open] fail to open file"))) + self.assertFalse(R._is_unsorted_error(self.error(None))) + + def test_an_unrelated_index_failure_is_raised_not_sorted_around(self): + with tempfile.TemporaryDirectory() as td, \ + mock.patch.object(R, "_bgzip", side_effect=lambda s, t: t.write_bytes(b"x")), \ + mock.patch.object(R, "_needs_csi", return_value=False), \ + mock.patch.object(R, "_index", side_effect=self.error(b"fail to open file")), \ + mock.patch.object(R, "_sort") as sort: + with self.assertRaises(R.subprocess.CalledProcessError): + R.prepare_indexed_vcf(SMALL_VCF, Path(td) / "indexed") + sort.assert_not_called() + + def test_sorting_without_bcftools_is_a_clear_error(self): + with tempfile.TemporaryDirectory() as td, \ + mock.patch.object(R.shutil, "which", return_value=None): + with self.assertRaisesRegex(RuntimeError, "not coordinate-sorted"): + R._sort(SMALL_VCF, Path(td) / "out.vcf.gz", Path(td)) + + @unittest.skipUnless(R is not None and have_bgzip and have_tabix, "bgzip and tabix are required") class PrepareIndexedVcfTests(VerboseTestCase): """The VCF side's one-time cost, built for real and reported separately.""" + def test_a_sorted_vcf_is_indexed_without_sorting(self): + with tempfile.TemporaryDirectory() as td: + report = R.prepare_indexed_vcf(write_vcf(Path(td) / "sorted.vcf"), Path(td) / "indexed") + self.assertNotIn("sorted", report) + self.assertEqual(report["sortSeconds"], 0.0) + + @unittest.skipUnless(have_bcftools, "bcftools is required to sort") + def test_an_unsorted_vcf_is_sorted_before_indexing(self): + """test-1k interleaves contigs, which tabix refuses to index as it stands.""" + with tempfile.TemporaryDirectory() as td: + report = R.prepare_indexed_vcf(SMALL_VCF, Path(td) / "indexed") + self.assertIn("not coordinate-sorted", report["sorted"]) + self.assertGreater(report["sortSeconds"], 0.0) + self.assertTrue(report["indexedVcf"].endswith(".sorted.vcf.gz")) + self.assertTrue(Path(report["indexPath"]).is_file()) + self.assertAlmostEqual( + report["totalSetupSeconds"], + report["bgzipSeconds"] + report["sortSeconds"] + report["indexSeconds"]) + # Only the sorted copy is kept, so a reader cannot pick up the wrong one. + self.assertEqual(sorted(p.name for p in (Path(td) / "indexed").glob("*.vcf.gz")), + ["test-1k.sorted.vcf.gz"]) + + @unittest.skipUnless(have_bcftools, "bcftools is required to sort") def test_a_plain_vcf_is_bgzipped_and_indexed(self): with tempfile.TemporaryDirectory() as td: report = R.prepare_indexed_vcf(SMALL_VCF, Path(td) / "indexed") @@ -183,15 +238,17 @@ def test_a_plain_vcf_is_bgzipped_and_indexed(self): self.assertTrue(Path(report["indexPath"]).is_file()) self.assertGreater(report["indexBytes"], 0) self.assertAlmostEqual(report["totalSetupSeconds"], - report["bgzipSeconds"] + report["indexSeconds"]) + report["bgzipSeconds"] + report["sortSeconds"] + + report["indexSeconds"]) def test_an_already_bgzf_input_is_copied_not_recompressed(self): """Recompressing would charge the VCF side for work no user repeats.""" with tempfile.TemporaryDirectory() as td: source = Path(td) / "source.vcf.gz" - R._bgzip(SMALL_VCF, source) + R._bgzip(write_vcf(Path(td) / "source.vcf"), source) report = R.prepare_indexed_vcf(source, Path(td) / "indexed") self.assertEqual(report["compression"], "already-bgzf (copied)") + self.assertNotIn("sorted", report) self.assertEqual(Path(report["indexedVcf"]).read_bytes(), source.read_bytes()) def test_auto_picks_csi_for_a_contig_a_tbi_cannot_address(self): @@ -204,7 +261,8 @@ def test_auto_picks_csi_for_a_contig_a_tbi_cannot_address(self): def test_an_explicit_index_kind_is_honoured(self): with tempfile.TemporaryDirectory() as td: - report = R.prepare_indexed_vcf(SMALL_VCF, Path(td) / "indexed", index_kind="csi") + report = R.prepare_indexed_vcf(write_vcf(Path(td) / "sorted.vcf"), + Path(td) / "indexed", index_kind="csi") self.assertEqual(report["indexKind"], "csi") self.assertNotIn("indexKindReason", report) @@ -212,8 +270,8 @@ def test_an_explicit_index_kind_is_honoured(self): # --------------------------------------------------------------------------- # The indexed arms against a real index # --------------------------------------------------------------------------- -@unittest.skipUnless(cyvcf2_available and have_bgzip and have_tabix, - "cyvcf2, bgzip and tabix are required") +@unittest.skipUnless(cyvcf2_available and have_bgzip and have_tabix and have_bcftools, + "cyvcf2, bgzip, tabix and bcftools are required") class IndexedArmTests(VerboseTestCase): """The seek path must return exactly what the scan path returns.""" From 0ca2ff6db9fc270b02486cee11f33a6bc942d54a Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 20:15:52 +0200 Subject: [PATCH 4/9] Hand the SPARQL arms N-Triples, the engine options they read, and isolation Found by running the runner end to end in the image on vcf-bench-1, with every arm against a real test-1k conversion. With the .nt.gz the harness reuses from 13_query_cost as --rdf, three of the four engines failed: qlever-index exited 2, Comunica could not read it (HTTP 400), and the COTTAS builder raised KeyError: '.gz' -- which the runner did not catch, so it ended the whole run. Only HDT, which builds its own artifact, coped. Three fixes: - The artifact is materialized as N-Triples once, before any engine, with validation_runner.materialize_ntriples -- the step the validation stage has always taken. It is timed and reported as setup (rdfMaterialization); a plain .nt is used in place. - engine_options() passes the names the engines read (memory_gb, port, startup_timeout, query_timeout) instead of the runner's own flag names, which the engines had silently ignored, and offers the supplied artifact as artifact_path/artifact_format so HDT and COTTAS can query it natively. - An engine that fails to start, whatever the exception type, is recorded in setup.engineErrors against that engine while the other arms' results stand. The engine loop moves into _run_engine_arm; its behaviour is otherwise unchanged. Four new tests cover the decode, the plain-file path, the option names and a non-RuntimeError start-up failure. Co-Authored-By: Claude Opus 5.5 --- src/validation/regional_runner.py | 156 +++++++++++++++-------- test/test_regional_runner_driver_unit.py | 66 +++++++++- 2 files changed, 167 insertions(+), 55 deletions(-) diff --git a/src/validation/regional_runner.py b/src/validation/regional_runner.py index 3515413..53589c7 100644 --- a/src/validation/regional_runner.py +++ b/src/validation/regional_runner.py @@ -869,59 +869,28 @@ def record(arm: str, query_id: str, window: dict[str, Any], replicate: int, if engine_arms and args.rdf is None: raise SystemExit("--rdf is required when a SPARQL engine is among --arms") - for arm in engine_arms: - print(f"[arm] {arm}: preparing engine") - engine_raw = raw_dir / arm - engine_raw.mkdir(parents=True, exist_ok=True) - options = { - "query_timeout": args.query_timeout, - "rdf_format": args.rdf_format, - "qlever_memory_gb": args.qlever_memory_gb, - "qlever_port": args.qlever_port, - "qlever_startup_timeout": args.qlever_startup_timeout, - } - try: - engine = V.build_engine(arm, args.rdf, raw_dir=engine_raw, - scratch=args.scratch_dir, options=options) - except (ValueError, RuntimeError) as error: - setup.setdefault("engineErrors", {})[arm] = str(error) - print(f"[arm] {arm}: unavailable ({error})") - continue - - try: - with engine: - setup.setdefault("engineSetupSeconds", {})[arm] = engine.setup_seconds - print(f"[arm] {arm}: setup {engine.setup_seconds:.2f}s; " - f"{len(windows)} windows x {len(queries)} questions " - f"x {args.replicates} replicates") - with tempfile.TemporaryDirectory(dir=str(args.scratch_dir)) as rendered_dir: - rendered_root = Path(rendered_dir) - for query_id in queries: - template = regional_query_path(args.representation, query_id) - for window in windows: - rendered = rendered_root / f"{query_id}__{window['window_id']}.rq" - rendered.write_text(render_query(template, window), encoding="utf-8") - for replicate in range(1, args.replicates + 1): - envelope = engine.execute( - f"{query_id}__{window['window_id']}__r{replicate}", rendered) - if envelope["status"] != "PASS": - record(arm, query_id, window, replicate, None, - envelope["wallSeconds"], status="FAILED", - error=envelope.get("error") or "engine execution failed") - continue - try: - answer = normalize_regional( - query_id, Path(envelope["rawResult"])) - except (ValueError, KeyError, OSError) as error: - record(arm, query_id, window, replicate, None, - envelope["wallSeconds"], status="UNREADABLE", - error=str(error)) - continue - record(arm, query_id, window, replicate, answer, - envelope["wallSeconds"]) - except RuntimeError as error: - setup.setdefault("engineErrors", {})[arm] = str(error) - print(f"[arm] {arm}: failed ({error})") + with tempfile.TemporaryDirectory(prefix="regional-rdf-", dir=str(args.scratch_dir)) as rdf_scratch: + ntriples = None + if engine_arms: + # The engines take N-Triples, exactly as in the validation stage: a + # packaged artifact (.nt.gz, .nt.br, .hdt, .cottas) is materialized + # once, here, and the engines are handed the plain file. Passing the + # package straight through made qlever-index, Comunica and the COTTAS + # builder all fail on an .nt.gz, which is the artifact the harness + # reuses from 13_query_cost. The decode is one-time setup and is + # reported as such. + print(f"[setup] materializing {args.rdf_format} as N-Triples for the SPARQL arms") + started = time.monotonic() + ntriples, materialization = V.materialize_ntriples( + args.rdf, args.rdf_format, Path(rdf_scratch), + log_dir=raw_dir / "materialization", + ) + materialization["wallSeconds"] = time.monotonic() - started + materialization["ntriplesPath"] = str(ntriples) + setup["rdfMaterialization"] = materialization + + for arm in engine_arms: + _run_engine_arm(arm, args, ntriples, raw_dir, setup, windows, queries, record) _write_outputs(results_dir, args, rows, mismatches, setup, windows, queries, arms) @@ -936,6 +905,87 @@ def record(arm: str, query_id: str, window: dict[str, Any], replicate: int, return 0 +def engine_options(args: argparse.Namespace) -> dict[str, Any]: + """The option names the validation engines read, not the runner's own flags. + + The engines are validation_runner's, and they look up ``memory_gb``, + ``port``, ``startup_timeout``, ``query_timeout`` and -- for HDT and COTTAS -- + ``artifact_path`` / ``artifact_format``, which lets them query the supplied + artifact natively instead of rebuilding it from N-Triples. Passing the + runner's flag names instead made every one of those settings a silent no-op. + """ + return { + "query_timeout": args.query_timeout, + "memory_gb": args.qlever_memory_gb, + "port": args.qlever_port, + "startup_timeout": args.qlever_startup_timeout, + "artifact_path": str(args.rdf), + "artifact_format": args.rdf_format, + } + + +def _run_engine_arm( + arm: str, + args: argparse.Namespace, + ntriples: Path, + raw_dir: Path, + setup: dict[str, Any], + windows: list[dict[str, Any]], + queries: list[str], + record, +) -> None: + """Time one SPARQL engine on every window; record, never raise, its failures.""" + print(f"[arm] {arm}: preparing engine") + engine_raw = raw_dir / arm + engine_raw.mkdir(parents=True, exist_ok=True) + try: + engine = V.build_engine(arm, ntriples, raw_dir=engine_raw, + scratch=args.scratch_dir, options=engine_options(args)) + except (ValueError, RuntimeError) as error: + setup.setdefault("engineErrors", {})[arm] = str(error) + print(f"[arm] {arm}: unavailable ({error})") + return + + try: + with engine: + setup.setdefault("engineSetupSeconds", {})[arm] = engine.setup_seconds + print(f"[arm] {arm}: setup {engine.setup_seconds:.2f}s; " + f"{len(windows)} windows x {len(queries)} questions " + f"x {args.replicates} replicates") + with tempfile.TemporaryDirectory(dir=str(args.scratch_dir)) as rendered_dir: + rendered_root = Path(rendered_dir) + for query_id in queries: + template = regional_query_path(args.representation, query_id) + for window in windows: + rendered = rendered_root / f"{query_id}__{window['window_id']}.rq" + rendered.write_text(render_query(template, window), encoding="utf-8") + for replicate in range(1, args.replicates + 1): + envelope = engine.execute( + f"{query_id}__{window['window_id']}__r{replicate}", rendered) + if envelope["status"] != "PASS": + record(arm, query_id, window, replicate, None, + envelope["wallSeconds"], status="FAILED", + error=envelope.get("error") or "engine execution failed") + continue + try: + answer = normalize_regional( + query_id, Path(envelope["rawResult"])) + except (ValueError, KeyError, OSError) as error: + record(arm, query_id, window, replicate, None, + envelope["wallSeconds"], status="UNREADABLE", + error=str(error)) + continue + record(arm, query_id, window, replicate, answer, + envelope["wallSeconds"]) + except Exception as error: # noqa: BLE001 - one engine must not end the run + # Start-up is where engines fail, and they fail in their own ways: the + # COTTAS builder raised KeyError on an .nt.gz, for instance. Whatever + # the type, it is this engine's failure -- recorded against it, while + # the other arms' results stand. + setup.setdefault("engineErrors", {})[arm] = f"{type(error).__name__}: {error}" + print(f"[arm] {arm}: failed ({type(error).__name__}: {error})") + + def _sample_scan_windows( windows: list[dict[str, Any]], per_size: int ) -> list[dict[str, Any]]: diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver_unit.py index 4cee0ef..40f1f84 100644 --- a/test/test_regional_runner_driver_unit.py +++ b/test/test_regional_runner_driver_unit.py @@ -543,11 +543,12 @@ def setUp(self): def tearDown(self): self._tmp.cleanup() - def run_with(self, mode: str, build=None): + def run_with(self, mode: str, build=None, *, rdf=None, rdf_format="nt"): build = build or (lambda *a, **k: _StandInEngine(self.results, mode)) with mock.patch.object(R.V, "build_engine", side_effect=build): status, _ = run_main(base_argv(self.results, self.scratch, "qlever", - queries=self.QUERIES, rdf=self.graph)) + queries=self.QUERIES, rdf=rdf or self.graph, + rdf_format=rdf_format)) report = json.loads((self.results / "regional.json").read_text()) with (self.results / "regional.csv").open(encoding="utf-8") as handle: rows = list(csv.DictReader(handle)) @@ -580,6 +581,67 @@ def refuse(*args, **kwargs): self.assertEqual(report["setup"]["engineErrors"], {"qlever": "qlever-index not on PATH"}) self.assertEqual(rows, []) + def test_any_start_up_failure_is_the_engines_not_the_runs(self): + """Found on bench-1: the COTTAS builder raised KeyError on an .nt.gz.""" + class Crashes(_StandInEngine): + def __enter__(self): + raise KeyError(".gz") + + status, report, rows = self.run_with( + "agree", build=lambda *a, **k: Crashes(self.results, "agree")) + self.assertEqual(status, 0) + self.assertEqual(report["setup"]["engineErrors"], {"qlever": "KeyError: '.gz'"}) + self.assertEqual(rows, []) + + def test_a_packaged_graph_is_decoded_once_before_any_engine_sees_it(self): + """qlever-index, Comunica and the COTTAS builder all refuse an .nt.gz.""" + packed = self.root / "graph.nt.gz" + with gzip.open(packed, "wt", encoding="utf-8") as handle: + handle.write(' "o" .\n') + seen = {} + + def build(name, source, **kwargs): + seen["source"] = Path(source) + seen["readable"] = Path(source).read_text(encoding="utf-8") + seen["options"] = kwargs["options"] + return _StandInEngine(self.results, "agree") + + status, report, _ = self.run_with("agree", build=build, rdf=packed, rdf_format="nt.gz") + self.assertEqual(status, 0) + self.assertEqual(seen["source"].suffix, ".nt") + self.assertIn("", seen["readable"]) + self.assertTrue(report["setup"]["rdfMaterialization"]["materialized"]) + self.assertGreaterEqual(report["setup"]["rdfMaterialization"]["wallSeconds"], 0.0) + # The native artifact is still offered, so HDT and COTTAS can use one. + self.assertEqual(seen["options"]["artifact_path"], str(packed)) + self.assertEqual(seen["options"]["artifact_format"], "nt.gz") + + def test_a_plain_graph_is_used_in_place(self): + seen = {} + + def build(name, source, **kwargs): + seen["source"] = Path(source) + return _StandInEngine(self.results, "agree") + + self.run_with("agree", build=build) + self.assertEqual(seen["source"], self.graph) + + +@unittest.skipUnless(R is not None, "regional_runner must import") +class EngineOptionTests(VerboseTestCase): + def test_the_engines_are_given_the_option_names_they_read(self): + """The runner's own flag names were silently ignored by the engines.""" + args = R.build_parser().parse_args([ + "--vcf", "x.vcf", "--results-dir", "r", "--dataset-id", "d", + "--rdf", "g.hdt", "--rdf-format", "hdt", + "--qlever-memory-gb", "12", "--qlever-port", "7100", + "--qlever-startup-timeout", "60", "--query-timeout", "30", + ]) + self.assertEqual(R.engine_options(args), { + "query_timeout": 30, "memory_gb": 12, "port": 7100, "startup_timeout": 60, + "artifact_path": "g.hdt", "artifact_format": "hdt", + }) + if __name__ == "__main__": unittest.main() From 6ecd528a84bc51dda74b228e84385105160b47f3 Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 21:10:09 +0200 Subject: [PATCH 5/9] Decompress plain-gzip VCFs before bgzip; report tool failures cleanly The benchmark's derived slices are plain gzip. bgzip -c on one compressed it a second time and tabix could not parse the result, so the first real-slice run on bench-1 ended in a traceback. A gzip input that is not BGZF is now streamed through gzip into bgzip, and a failing external tool is reported as an error exit naming the command. Co-Authored-By: Claude Opus 5.5 --- src/validation/regional_runner.py | 37 ++++++++++++++++++++++- test/test_regional_runner_driver_unit.py | 38 ++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 1 deletion(-) diff --git a/src/validation/regional_runner.py b/src/validation/regional_runner.py index 53589c7..1765caa 100644 --- a/src/validation/regional_runner.py +++ b/src/validation/regional_runner.py @@ -29,6 +29,7 @@ import argparse import csv +import gzip import json import random import shutil @@ -320,7 +321,24 @@ def _bgzip(source: Path, target: Path) -> None: """ if shutil.which("bgzip"): with target.open("wb") as handle: - subprocess.run(["bgzip", "-c", str(source)], stdout=handle, check=True) + if _is_gzip(source): + # A plain-gzip VCF has to be decompressed first. `bgzip -c` on + # it compresses the gzip bytes a second time, and tabix then + # fails to parse the result ("was wrong -p [type] used?"). The + # benchmark's derived slices are plain gzip, so this is the + # common case, not an edge one. + with gzip.open(source, "rb") as src: + proc = subprocess.Popen(["bgzip", "-c"], stdin=subprocess.PIPE, + stdout=handle) + assert proc.stdin is not None + try: + shutil.copyfileobj(src, proc.stdin, length=1024 * 1024) + finally: + proc.stdin.close() + if proc.wait() != 0: + raise subprocess.CalledProcessError(proc.returncode, ["bgzip", "-c"]) + else: + subprocess.run(["bgzip", "-c", str(source)], stdout=handle, check=True) return if shutil.which("bcftools"): subprocess.run( @@ -331,6 +349,15 @@ def _bgzip(source: Path, target: Path) -> None: raise RuntimeError("neither bgzip nor bcftools is available to produce a BGZF file") +def _is_gzip(path: Path) -> bool: + """True for any gzip stream, BGZF or not (BGZF is checked by is_bgzf).""" + try: + with path.open("rb") as handle: + return handle.read(2) == b"\x1f\x8b" + except OSError: + return False + + def _needs_csi(bgzf_path: Path) -> bool: """True when any contig is longer than a .tbi can address.""" limit = 2 ** 29 - 1 @@ -1114,6 +1141,14 @@ def main(argv: list[str] | None = None) -> int: args.scratch_dir.mkdir(parents=True, exist_ok=True) try: return run(args) + except subprocess.CalledProcessError as error: + # A tool the runner shells out to failed: say which one and what it + # printed, instead of ending in a traceback. + detail = error.stderr.decode("utf-8", "replace").strip() if isinstance( + error.stderr, bytes) else (error.stderr or "").strip() + print(f"error: {' '.join(map(str, error.cmd))} exited {error.returncode}" + + (f": {detail.splitlines()[-1]}" if detail else ""), file=sys.stderr) + return 2 except (RuntimeError, ValueError, OSError) as error: print(f"error: {error}", file=sys.stderr) return 2 diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver_unit.py index 40f1f84..3e4f481 100644 --- a/test/test_regional_runner_driver_unit.py +++ b/test/test_regional_runner_driver_unit.py @@ -93,6 +93,33 @@ def test_plain_gzip_is_not_bgzf(self): handle.write(SMALL_VCF.read_bytes()) self.assertFalse(R.is_bgzf(path)) + @unittest.skipUnless(have_bgzip and have_tabix, "bgzip and tabix are required") + def test_a_plain_gzip_input_is_decompressed_before_bgzip(self): + """Found on bench-1: the derived slices are plain gzip, and bgzip -c on + one compressed it twice, which tabix could not parse.""" + with tempfile.TemporaryDirectory() as td: + plain = Path(td) / "plain.vcf.gz" + with gzip.open(plain, "wb") as handle: + handle.write(write_vcf(Path(td) / "src.vcf").read_bytes()) + target = Path(td) / "out" / "plain.vcf.gz" + target.parent.mkdir() + R._bgzip(plain, target) + self.assertTrue(R.is_bgzf(target)) + with gzip.open(target, "rt", encoding="utf-8") as handle: + self.assertTrue(handle.readline().startswith("##fileformat=VCFv4.2")) + report = R.prepare_indexed_vcf(plain, Path(td) / "indexed") + self.assertEqual(report["compression"], "bgzip") + self.assertTrue(Path(report["indexPath"]).is_file()) + + def test_gzip_is_recognised_whether_or_not_it_is_bgzf(self): + with tempfile.TemporaryDirectory() as td: + plain = Path(td) / "plain.vcf.gz" + with gzip.open(plain, "wb") as handle: + handle.write(b"x") + self.assertTrue(R._is_gzip(plain)) + self.assertFalse(R._is_gzip(SMALL_VCF)) + self.assertFalse(R._is_gzip(Path(td) / "missing.gz")) + @unittest.skipUnless(have_bgzip, "bgzip is required to produce a real BGZF file") def test_real_bgzip_output_is_bgzf(self): with tempfile.TemporaryDirectory() as td: @@ -466,6 +493,17 @@ def test_a_sparql_arm_without_a_graph_is_refused(self): with self.assertRaisesRegex(SystemExit, "--rdf is required"): run_main(base_argv(self.results, self.scratch, "qlever")) + def test_a_failing_tool_is_an_error_exit_not_a_traceback(self): + error = R.subprocess.CalledProcessError( + 1, ["tabix", "-p", "vcf", "x.vcf.gz"], stderr=b"[E::get_intv] Failed to parse TBX_VCF") + from contextlib import redirect_stderr + err = StringIO() + with mock.patch.object(R, "run", side_effect=error), redirect_stderr(err): + status, _ = run_main(base_argv(self.results, self.scratch, "cyvcf2-scan")) + self.assertEqual(status, 2) + self.assertIn("tabix -p vcf x.vcf.gz exited 1: [E::get_intv] Failed to parse TBX_VCF", + err.getvalue()) + def test_a_vcf_without_records_is_an_error_exit_not_a_traceback(self): path = self.root / "empty.vcf" path.write_text("##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n", From 325e3483459efd08fd33f0d9fe1ff2f8f1d4df51 Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 21:27:48 +0200 Subject: [PATCH 6/9] Time slow SPARQL engines on sampled windows, like the scan arm Comunica, HDT and COTTAS cost tens of seconds a question on the 10k-record input; timed on every window, three replicates, they would take about a day and a half per scale. --thin-arms names arms timed on only --scan-windows-per-size windows of each size, and regional.json records which were. The default keeps the old behaviour: only the scan arm is thin. Also corrects the documented scan cost, measured on vcf-bench-1 at 0.64 s a question on the HG005 slice rather than the estimated 11 s. Co-Authored-By: Claude Opus 5.5 --- docs/validation.md | 15 ++++++++----- src/validation/regional_runner.py | 27 ++++++++++++++++++------ test/test_regional_runner_driver_unit.py | 26 +++++++++++++++++++++++ 3 files changed, 56 insertions(+), 12 deletions(-) diff --git a/docs/validation.md b/docs/validation.md index 91d6afe..627080e 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -462,11 +462,16 @@ sees identical regions. **The scan arm is timed on fewer windows than the others.** Its cost does not depend on the window -- it reads the file whichever region is asked for -- so -timing it on every window measures one number repeatedly at roughly 11 s a go. -Correctness is still checked on every window: the equality reference is a single -whole-file pass that fills all of them at once, and it is deliberately not -timed, because a pass that answers eighty windows is not what a one-question -user pays for. +timing it on every window measures one number repeatedly (0.64 s a question on +the 100,000-record HG005 slice, on vcf-bench-1). `--scan-windows-per-size` sets +how many windows of each size it gets. `--thin-arms` names other arms to time +the same way; the benchmark uses it for Comunica, HDT and COTTAS, which cost tens +of seconds a question, so timing them on all eighty windows would take about +a day and a half. The equality reference still covers every window: it +is a single whole-file pass that fills all of them at once, and it is +deliberately not timed, because a pass that answers eighty windows is not what +a one-question user pays for. A thin arm's own agreement is checked on the +windows it is timed on. Results are compared for exact equality before any timing is reported. The runner exits non-zero when arms disagree, because a speed number from arms that diff --git a/src/validation/regional_runner.py b/src/validation/regional_runner.py index 1765caa..f4bdc1c 100644 --- a/src/validation/regional_runner.py +++ b/src/validation/regional_runner.py @@ -791,6 +791,10 @@ def run(args: argparse.Namespace) -> int: unknown = [a for a in arms if a not in VCF_ARMS and a not in V.SPARQL_ENGINES] if unknown: raise SystemExit(f"unknown arm(s): {', '.join(unknown)}") + thin_arms = {arm.strip() for arm in args.thin_arms.split(",") if arm.strip()} + unknown_thin = sorted(a for a in thin_arms if a not in VCF_ARMS and a not in V.SPARQL_ENGINES) + if unknown_thin: + raise SystemExit(f"unknown arm(s) in --thin-arms: {', '.join(unknown_thin)}") queries = [q.strip() for q in args.queries.split(",") if q.strip()] unknown_queries = [q for q in queries if q not in REGIONAL_QUERIES] if unknown_queries: @@ -873,12 +877,14 @@ def record(arm: str, query_id: str, window: dict[str, Any], replicate: int, # The scan arm's cost does not depend on the window -- it reads the whole # file whichever region is asked for. Timing it on all 80 windows would - # measure one number 80 times, at roughly 11 s a go. It is sampled instead, - # and the sample size is reported so the thinness is visible. - scan_windows = _sample_scan_windows(windows, args.scan_windows_per_size) + # measure one number 80 times. It is sampled instead, and the sample size is + # reported so the thinness is visible. The same applies to any arm named in + # --thin-arms: Comunica, HDT and COTTAS cost tens of seconds a question, so + # timing them on every window would take a day and a half per scale. + thin_windows = _sample_scan_windows(windows, args.scan_windows_per_size) for arm in [a for a in arms if a in VCF_ARMS]: - targets = scan_windows if arm == "cyvcf2-scan" else windows + targets = thin_windows if arm in thin_arms else windows print(f"[arm] {arm}: {len(targets)} windows x {len(queries)} questions " f"x {args.replicates} replicates") for query_id in queries: @@ -917,7 +923,8 @@ def record(arm: str, query_id: str, window: dict[str, Any], replicate: int, setup["rdfMaterialization"] = materialization for arm in engine_arms: - _run_engine_arm(arm, args, ntriples, raw_dir, setup, windows, queries, record) + targets = thin_windows if arm in thin_arms else windows + _run_engine_arm(arm, args, ntriples, raw_dir, setup, targets, queries, record) _write_outputs(results_dir, args, rows, mismatches, setup, windows, queries, arms) @@ -1081,6 +1088,8 @@ def _write_outputs( "queries": queries, "replicates": args.replicates, "scanWindowsPerSize": args.scan_windows_per_size, + "thinArms": [a for a in arms if a in { + t.strip() for t in args.thin_arms.split(",") if t.strip()}], "windowCount": len(windows), "setup": setup, "regionSemantics": REGION_SEMANTICS, @@ -1121,8 +1130,12 @@ def build_parser() -> argparse.ArgumentParser: default=",".join(str(s) for s in DEFAULT_WINDOW_SIZES)) parser.add_argument("--windows-per-size", type=int, default=DEFAULT_WINDOWS_PER_SIZE) parser.add_argument("--scan-windows-per-size", type=int, default=3, - help="windows timed for cyvcf2-scan, whose cost does not " - "depend on the window (default: 3)") + help="windows timed for each arm in --thin-arms " + "(default: 3)") + parser.add_argument("--thin-arms", default="cyvcf2-scan", + help="arms timed on only --scan-windows-per-size windows " + "of each size; cyvcf2-scan by default, whose cost " + "does not depend on the window") parser.add_argument("--replicates", type=int, default=3) parser.add_argument("--seed", type=int, default=DEFAULT_SEED) parser.add_argument("--index-kind", choices=("auto", "tbi", "csi"), default="auto") diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver_unit.py index 3e4f481..5e89b0a 100644 --- a/test/test_regional_runner_driver_unit.py +++ b/test/test_regional_runner_driver_unit.py @@ -599,6 +599,32 @@ def test_an_agreeing_engine_is_timed_on_every_window(self): self.assertEqual(len(rows), 4 * 2) # every window, two questions self.assertTrue(all(row["agrees_with_reference"] == "True" for row in rows)) + def test_a_thin_engine_is_timed_on_the_sampled_windows_only(self): + """Comunica, HDT and COTTAS cost tens of seconds a question, so the + harness times them as thinly as the scan arm.""" + with mock.patch.object(R.V, "build_engine", + side_effect=lambda *a, **k: _StandInEngine(self.results, "agree")): + status, _ = run_main(base_argv(self.results, self.scratch, "cyvcf2-scan,qlever", + queries=self.QUERIES, rdf=self.graph, + thin_arms="qlever")) + self.assertEqual(status, 0) + report = json.loads((self.results / "regional.json").read_text()) + self.assertEqual(report["thinArms"], ["qlever"]) + with (self.results / "regional.csv").open(encoding="utf-8") as handle: + rows = list(csv.DictReader(handle)) + per_arm = {arm: {row["window_id"] for row in rows if row["arm"] == arm} + for arm in ("qlever", "cyvcf2-scan")} + # One window per size for the thin arm; the scan arm, not named as + # thin here, is timed on both. + self.assertEqual(len(per_arm["qlever"]), 2) + self.assertEqual(len(per_arm["cyvcf2-scan"]), 4) + self.assertLess(per_arm["qlever"], per_arm["cyvcf2-scan"]) + self.assertTrue(all(row["agrees_with_reference"] == "True" for row in rows)) + + def test_an_unknown_thin_arm_is_refused(self): + with self.assertRaisesRegex(SystemExit, "--thin-arms"): + run_main(base_argv(self.results, self.scratch, "cyvcf2-scan", thin_arms="duckdb")) + def test_an_engine_failure_is_recorded_per_execution(self): status, report, rows = self.run_with("fail") self.assertEqual(status, 0) From 8ba2a4fc1c20e2d5f27d88b9dc5b25e140ba171f Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 21:41:51 +0200 Subject: [PATCH 7/9] Correct the stated cost of the SPARQL engines on regional questions The previous commit justified --thin-arms with the 23-45 s these engines take on 13's whole-file questions. Measured on vcf-bench-1, a region-restricted question costs Comunica, HDT and COTTAS about one second on the 10k-record fixture, so the benchmark times every engine on every window. The option stays for a genuinely slow arm; its default is unchanged. Co-Authored-By: Claude Opus 5.5 --- docs/validation.md | 8 +++++--- src/validation/regional_runner.py | 5 ++--- test/test_regional_runner_driver_unit.py | 3 +-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/validation.md b/docs/validation.md index 627080e..a4fb001 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -465,9 +465,11 @@ depend on the window -- it reads the file whichever region is asked for -- so timing it on every window measures one number repeatedly (0.64 s a question on the 100,000-record HG005 slice, on vcf-bench-1). `--scan-windows-per-size` sets how many windows of each size it gets. `--thin-arms` names other arms to time -the same way; the benchmark uses it for Comunica, HDT and COTTAS, which cost tens -of seconds a question, so timing them on all eighty windows would take about -a day and a half. The equality reference still covers every window: it +the same way, for an arm too slow to time on every window. The benchmark does +not need it: a region-restricted question costs Comunica, HDT and COTTAS about +one second on the 10,000-record fixture (vcf-bench-1), against 23-45 s for the +whole-file questions above, so every engine is timed on every window. The +equality reference still covers every window: it is a single whole-file pass that fills all of them at once, and it is deliberately not timed, because a pass that answers eighty windows is not what a one-question user pays for. A thin arm's own agreement is checked on the diff --git a/src/validation/regional_runner.py b/src/validation/regional_runner.py index f4bdc1c..e692b6d 100644 --- a/src/validation/regional_runner.py +++ b/src/validation/regional_runner.py @@ -878,9 +878,8 @@ def record(arm: str, query_id: str, window: dict[str, Any], replicate: int, # The scan arm's cost does not depend on the window -- it reads the whole # file whichever region is asked for. Timing it on all 80 windows would # measure one number 80 times. It is sampled instead, and the sample size is - # reported so the thinness is visible. The same applies to any arm named in - # --thin-arms: Comunica, HDT and COTTAS cost tens of seconds a question, so - # timing them on every window would take a day and a half per scale. + # reported so the thinness is visible. --thin-arms extends the same + # treatment to any arm too slow to time on every window. thin_windows = _sample_scan_windows(windows, args.scan_windows_per_size) for arm in [a for a in arms if a in VCF_ARMS]: diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver_unit.py index 5e89b0a..235283e 100644 --- a/test/test_regional_runner_driver_unit.py +++ b/test/test_regional_runner_driver_unit.py @@ -600,8 +600,7 @@ def test_an_agreeing_engine_is_timed_on_every_window(self): self.assertTrue(all(row["agrees_with_reference"] == "True" for row in rows)) def test_a_thin_engine_is_timed_on_the_sampled_windows_only(self): - """Comunica, HDT and COTTAS cost tens of seconds a question, so the - harness times them as thinly as the scan arm.""" + """--thin-arms gives any slow arm the scan arm's sampled windows.""" with mock.patch.object(R.V, "build_engine", side_effect=lambda *a, **k: _StandInEngine(self.results, "agree")): status, _ = run_main(base_argv(self.results, self.scratch, "cyvcf2-scan,qlever", From c0b07d8717c4041b05ce8552e4c38d1400d8a756 Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 22:21:35 +0200 Subject: [PATCH 8/9] Keep the regional runner's tests out of the CI suite The regional runner is a performance investigation, not part of conversion or validation, so its tests should not run in the normal test workflow. They are renamed to miss CI's test_*_unit.py pattern, as cross_engine_agreement.py and test_cottas_tool.py already do, and test/README.md says how to run them. Co-Authored-By: Claude Opus 5.5 --- test/README.md | 11 +++++++++++ ...egional_runner_unit.py => test_regional_runner.py} | 0 ..._driver_unit.py => test_regional_runner_driver.py} | 2 +- 3 files changed, 12 insertions(+), 1 deletion(-) rename test/{test_regional_runner_unit.py => test_regional_runner.py} (100%) rename test/{test_regional_runner_driver_unit.py => test_regional_runner_driver.py} (99%) diff --git a/test/README.md b/test/README.md index 4b6fa46..5c8e117 100644 --- a/test/README.md +++ b/test/README.md @@ -62,6 +62,17 @@ This repository uses `unittest` (Python standard library) to isolate orchestrati must agree with the Python oracle and not merely with the other engines. - Pass a comma-separated subset as the first argument to narrow it. +- `test/test_regional_runner.py` and `test/test_regional_runner_driver.py` + - Cover `src/validation/regional_runner.py`, the indexed regional-access + performance investigation (see "Indexed regional access" in + `docs/validation.md`): agreement between the access paths, index + preparation, and the driver end to end. + - Not run in CI or by the normal test workflow. The names deliberately miss + the `test_*_unit.py` pattern, because the runner is a performance tool, not + part of conversion or validation. Run them explicitly when changing it, + inside the image so the index and bcftools tests do not skip: + `python -m unittest test.test_regional_runner test.test_regional_runner_driver`. + - `test/test_validation_logic_unit.py` - Mutation tests over the validator's pure comparison layer, run on the host without cyvcf2 or Docker. diff --git a/test/test_regional_runner_unit.py b/test/test_regional_runner.py similarity index 100% rename from test/test_regional_runner_unit.py rename to test/test_regional_runner.py diff --git a/test/test_regional_runner_driver_unit.py b/test/test_regional_runner_driver.py similarity index 99% rename from test/test_regional_runner_driver_unit.py rename to test/test_regional_runner_driver.py index 235283e..eca2af7 100644 --- a/test/test_regional_runner_driver_unit.py +++ b/test/test_regional_runner_driver.py @@ -1,6 +1,6 @@ """The regional runner's driver, index preparation and failure paths. -test_regional_runner_unit.py covers agreement between the arms' folding and +test_regional_runner.py covers agreement between the arms' folding and classification logic. This module covers what surrounds it: building the bgzip + index copy the indexed arms seek through, dispatching one timed execution to the right arm, and the driver end to end -- the files it writes, From e13f974b73cc13d2c569bc8b889cc93772efc67e Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 24 Sep 2026 22:21:45 +0200 Subject: [PATCH 9/9] Say in the docs that the regional runner is a standalone investigation Co-Authored-By: Claude Opus 5.5 --- docs/validation.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/validation.md b/docs/validation.md index a4fb001..cb31e7d 100644 --- a/docs/validation.md +++ b/docs/validation.md @@ -420,8 +420,10 @@ none of Q1-Q13 is coordinate-restricted. That is internally consistent, and it is the one VCF access mode nobody uses for a selective question: real VCF work seeks, through a `bgzip` + `tabix` index. -`validation/regional_runner.py` adds that arm. It asks five region-restricted -questions -- record count, allele shape, Ti/Tv, FILTER distribution, per-sample +`validation/regional_runner.py` adds that arm as a standalone performance +investigation: nothing in conversion or `--mode validation` calls it, and its +tests (`test/test_regional_runner*.py`) are outside the CI suite. It asks five +region-restricted questions -- record count, allele shape, Ti/Tv, FILTER distribution, per-sample genotype classes -- of every access path: | Arm | Access path |