Add an indexed regional-access arm for the retrieval comparison

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/*
 

docs/validation.md

@@ -413,6 +413,74 @@ 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 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 |
+| --- | --- |
+| `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 (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, 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
+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
+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

src/validation/queries/regional/common/r01_region_record_count.rq

@@ -0,0 +1,18 @@
+PREFIX vcfc: <https://w3id.org/vcf-core/vocab#>
+PREFIX xsd:  <http://www.w3.org/2001/XMLSchema#>
+
+# 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}})
+}

src/validation/queries/regional/common/r02_region_variant_shape_counts.rq

@@ -0,0 +1,41 @@
+PREFIX vcfc: <https://w3id.org/vcf-core/vocab#>
+PREFIX xsd:  <http://www.w3.org/2001/XMLSchema#>
+
+# 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

src/validation/queries/regional/common/r03_region_titv.rq

@@ -0,0 +1,18 @@
+PREFIX vcfc: <https://w3id.org/vcf-core/vocab#>
+PREFIX xsd:  <http://www.w3.org/2001/XMLSchema#>
+
+# 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)
+}

src/validation/queries/regional/common/r04_region_filter_distribution.rq

@@ -0,0 +1,16 @@
+PREFIX vcfc: <https://w3id.org/vcf-core/vocab#>
+PREFIX xsd:  <http://www.w3.org/2001/XMLSchema#>
+
+# 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

src/validation/queries/regional/expanded/r05_region_sample_genotype_counts.rq

@@ -0,0 +1,30 @@
+PREFIX vcfc: <https://w3id.org/vcf-core/vocab#>
+PREFIX xsd:  <http://www.w3.org/2001/XMLSchema#>
+
+# 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