From 20d2cbb6f5c7fd0408d88b910abf20c27df7f995 Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Wed, 16 Sep 2026 10:43:07 +0200 Subject: [PATCH] Count FORMAT value items in the oracle, not only INFO ones A run over any VCF with a positional FORMAT key failed validation. The oracle reported the entire value-item layer as unexpected -- 401,606 extra hasValueItem/forAllele/FieldValueItem rows on a 100k-record GIAB file -- and 13_query_cost failed 3/3 replicates at its large scale. _emit_value_items is called from two places: the INFO pass, and append_expanded_sample_rdf for FORMAT. emitted_record_counters only ever walked parse_info_entries, so the FORMAT half was never expected. format_numbers was already a parameter of that function and was never read once -- the counting had simply never been written. The layer is expanded-only, because append_expanded_sample_rdf is the only emitter of it, so it is returned separately and merged by expected_census under representation == "expanded", exactly as the genotype layer already is. The emitter's two preconditions are mirrored rather than re-derived: at least one ALT, and split_value_items yielding items (a missing cell yields none). Why the fixtures never caught it: test-1k.vcf and test-10k.vcf declare no positional FORMAT key at all, while HG005 declares AD and ADALL as Number=R. The bug needed real data to appear. Adds EmittedTermCensusCoverageTests, which diffs every vocab term the emitter produces against everything the oracle counts, classifying each use by its position in the triple so vocabulary objects are not mistaken for census terms. This is the guard for the whole bug class -- it has now bitten twice, here and in the INFO items under raw INFO. It ships with a waiver list of 35 terms that are genuinely unmodelled today, grouped by feature: local alleles, phase sets, SV/confidence intervals, gVCF reference blocks, tandem repeats and Number=M base modifications. Each waived entry is a file shape that cannot pass validation, so the list is a to-do, not a licence. A second test fails if a waived term later becomes counted, so the list cannot rot. Co-Authored-By: Claude Opus 5 --- src/validation/validation_runner.py | 59 ++++++++++ test/test_validation_oracle_unit.py | 162 ++++++++++++++++++++++++++++ 2 files changed, 221 insertions(+) diff --git a/src/validation/validation_runner.py b/src/validation/validation_runner.py index c3962de..86977bc 100644 --- a/src/validation/validation_runner.py +++ b/src/validation/validation_runner.py @@ -583,6 +583,11 @@ def emitted_record_counters( assembly_contig_ids: set[str] = set() reference_alleles = alt_alleles = 0 value_items = value_item_alleles = tuple_items = 0 + # FORMAT value items are counted separately: the emitter decomposes them + # only in the expanded per-sample representation, so expected_census merges + # them for that profile alone, exactly as it does the genotype layer. + format_items = format_item_alleles = format_tuple_items = 0 + format_item_predicates: Counter[str] = Counter() genotypes = genotype_calls = called_alleles = 0 for row in rows: @@ -650,6 +655,43 @@ def emitted_record_counters( if samples: format_keys = (row[8].split(":") if len(row) > 8 and row[8] else []) payloads = row[9 : 9 + len(samples)] + # A positional FORMAT cell is decomposed into vcfc:FieldValueItems + # exactly as a positional INFO one is -- _emit_value_items is called + # for both, from append_expanded_sample_rdf and from the INFO pass. + # Counting only the INFO side made every file with a positional + # FORMAT key (AD, ADALL, PL -- i.e. most real VCFs) fail validation + # with the whole item layer reported as unexpected extra rows. + # + # The emitter requires at least one ALT before it decomposes, and + # split_value_items returns nothing for a missing cell, so both + # conditions are mirrored rather than re-derived. + if alt_count and format_keys: + for payload in payloads: + fields = payload.split(":") if payload else [] + for key_index, key in enumerate(format_keys): + cell = fields[key_index] if key_index < len(fields) else "" + if not cell: + continue + number = format_numbers.get(key, ".") + if not version.is_positional(key, number): + continue + items = vocab.split_value_items(cell) + if not items: + continue + format_items += len(items) + if version.tuple_arity(key) is not None: + format_tuple_items += len(items) + for index in range(len(items)): + link = version.value_item_link(key, number, index) + if ( + link.allele_index is not None + and link.allele_index in allele_uris + ): + format_item_alleles += 1 + elif link.genotype_index is not None: + format_item_predicates["forGenotypeIndex"] += 1 + elif link.gt_allele_index is not None: + format_item_predicates["forGTAlleleIndex"] += 1 if "GT" in format_keys: gt_index = format_keys.index("GT") for payload in payloads: @@ -696,11 +738,22 @@ def emitted_record_counters( genotype_predicates[name] += genotype_calls genotype_predicates["calledAllele"] += called_alleles + format_item_classes: Counter[str] = Counter() + if format_items: + format_item_classes["FieldValueItem"] += format_items + for name in ("hasValueItem", "valueIndex", "itemValue"): + format_item_predicates[name] += format_items + format_item_predicates["forAllele"] += format_item_alleles + format_item_predicates["tupleArity"] += format_tuple_items + return { "emittedRecordClasses": dict(classes), "emittedRecordPredicates": dict(predicates), "emittedGenotypeClasses": dict(genotype_classes), "emittedGenotypePredicates": dict(genotype_predicates), + "emittedFormatItemClasses": dict(format_item_classes), + "emittedFormatItemPredicates": dict(format_item_predicates), + "formatValueItemCount": format_items, "alleleCount": total_alleles, "valueItemCount": value_items, "genotypeCount": genotypes, @@ -870,6 +923,12 @@ def expected_census( classes[f"{VCFC}{class_name}"] = classes.get(f"{VCFC}{class_name}", 0) + count for name, count in parser.get("emittedGenotypePredicates", {}).items(): predicates[f"{VCFC}{name}"] = predicates.get(f"{VCFC}{name}", 0) + count + # The FORMAT item layer is emitted by append_expanded_sample_rdf, so + # it exists only here -- the condensed profile never materializes it. + for class_name, count in parser.get("emittedFormatItemClasses", {}).items(): + classes[f"{VCFC}{class_name}"] = classes.get(f"{VCFC}{class_name}", 0) + count + for name, count in parser.get("emittedFormatItemPredicates", {}).items(): + predicates[f"{VCFC}{name}"] = predicates.get(f"{VCFC}{name}", 0) + count classes[f"{VCFC}SampleCall"] = records * samples classes[f"{VCFC}FormatFieldValue"] = parser["formatValueSlots"] predicates[f"{VCFC}hasSampleCall"] = records * samples diff --git a/test/test_validation_oracle_unit.py b/test/test_validation_oracle_unit.py index 374a927..3d00aac 100644 --- a/test/test_validation_oracle_unit.py +++ b/test/test_validation_oracle_unit.py @@ -788,3 +788,165 @@ def test_failure_codes_are_counted_and_matched_against_declarations(self): def test_a_truncated_record_is_treated_as_missing_rather_than_crashing(self): self.assertEqual(V.record_decomposition_counts(["chr1", "1"], set()), (0, 0, 0)) + + +class FormatValueItemCensusTests(VerboseTestCase): + """The oracle must count FORMAT value items, not only INFO ones. + + Regression: _emit_value_items is called from two places -- the INFO pass and + append_expanded_sample_rdf -- but emitted_record_counters only walked the + INFO entries. format_numbers was even passed in and never read. Any file with + a positional FORMAT key (AD, ADALL, PL -- most real VCFs) then validated with + the whole item layer reported as unexpected extra rows: 401,606 of them on a + 100k-record GIAB file, which failed 3/3 replicates of 13_query_cost. + """ + + def _counters(self, rows, samples=("S1",), format_numbers=None, info_numbers=None): + import sys + sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + import vcf_rdfizer_vocab as vocab + version, _ = vocab.resolve_vcf_version("4.2") + return V.emitted_record_counters( + [list(r) for r in rows], + list(samples), + version=version, + contig_ids={"chr1"}, + alt_declaration_ids=set(), + info_numbers=info_numbers or {"DP": "1"}, + format_numbers=format_numbers or {"AD": "R", "GT": "1"}, + ) + + ROW = ["chr1", "100", ".", "A", "G", "50", "PASS", "DP=10", "GT:AD", "0/1:5,7"] + + def test_positional_format_cell_is_counted(self): + out = self._counters([self.ROW]) + self.assertEqual(out["emittedFormatItemClasses"].get("FieldValueItem"), 2) + for name in ("hasValueItem", "valueIndex", "itemValue"): + self.assertEqual(out["emittedFormatItemPredicates"].get(name), 2, name) + # Number=R indexes ref then each ALT, so both items name an allele. + self.assertEqual(out["emittedFormatItemPredicates"].get("forAllele"), 2) + + def test_format_items_are_kept_apart_from_info_items(self): + """They merge under different profiles, so they cannot share a counter.""" + out = self._counters([self.ROW]) + self.assertEqual(out["valueItemCount"], 0) + self.assertEqual(out["formatValueItemCount"], 2) + + def test_non_positional_format_key_contributes_nothing(self): + out = self._counters([self.ROW], format_numbers={"AD": "1", "GT": "1"}) + self.assertFalse(out["emittedFormatItemClasses"]) + + def test_missing_cell_contributes_nothing(self): + row = list(self.ROW); row[9] = "0/1:." + self.assertFalse(self._counters([row])["emittedFormatItemClasses"]) + + def test_no_alt_contributes_nothing(self): + """The emitter requires at least one ALT before it decomposes.""" + row = ["chr1", "100", ".", "A", ".", "50", "PASS", "DP=10", "GT:AD", "0/0:5"] + self.assertFalse(self._counters([row])["emittedFormatItemClasses"]) + + def test_the_layer_is_expanded_only(self): + """append_expanded_sample_rdf emits it; the condensed profile does not.""" + parser = self._counters([self.ROW]) + self.assertTrue(parser["emittedFormatItemClasses"]) + + +class EmittedTermCensusCoverageTests(VerboseTestCase): + """Every vocab term the emitter produces must be counted, or waived here. + + This is the guard for a whole class of bug: the emitter grows a predicate or + class, the oracle never learns to expect it, and every file using that + feature fails validation with extra rows. It has happened twice -- the INFO + value items under raw INFO, and the FORMAT value items. + + A term in KNOWN_UNMODELLED is a deliberate, documented gap. Shrinking the set + is progress. Growing it needs a reason in the commit message, because each + entry is a file shape that cannot pass validation today. + """ + + # Grouped by the feature that produces them. + KNOWN_UNMODELLED = { + # VCF 4.5 local alleles (LA/LR/LG) + "LocalAlleleSet", "LocalAlleleMembership", "hasLocalAlleleSet", + "hasLocalAlleleMembership", "hasLocalAllele", "localAllele", "localIndex", + # phase sets (PS/PSL/PSO/PSQ) + "PhaseSet", "inPhaseSet", "phaseSetId", "phaseSetName", "phaseSetOrdinal", + "phaseSetQuality", + # structural variants and confidence intervals + "VariantEvent", "inEvent", "eventType", "svClaim", + "ConfidenceInterval", "ciLower", "ciUpper", "endPosition", + # gVCF reference blocks + "ReferenceBlock", "isReferenceBlockStart", "referenceBlockLength", + # tandem repeats + "TandemRepeatAllele", "RepeatSequence", "hasRepeatSequence", + "repeatSequenceCount", "repeatSequenceIndex", + # Number=M base modifications + "BaseModification", "forBaseModification", "modifiedBaseOffset", + "modifiedResidue", + # raw carriers, deliberately not part of the census + "sampleDataRaw", "sampleFilter", + } + + @staticmethod + def _emitted_terms(source: str): + """Classify each _vocab(X) use by its position in the triple.""" + import re + calls, buf, depth = [], [], 0 + for line in source.split("\n"): + if "emit(" in line and depth == 0: + buf, depth = [line], line.count("(") - line.count(")") + if depth <= 0: + calls.append(" ".join(buf)); buf, depth = [], 0 + elif depth > 0: + buf.append(line); depth += line.count("(") - line.count(")") + if depth <= 0: + calls.append(" ".join(buf)); buf, depth = [], 0 + terms = set() + for call in calls: + raw = re.findall( + r"(RDF_TYPE_URI)|_vocab\(['\"]([A-Za-z_][A-Za-z0-9_]*)['\"]\)", call + ) + seq = [("TYPE", None) if a else ("V", b) for a, b in raw] + for i, (kind, name) in enumerate(seq): + if kind != "V": + continue + # A term preceded by another term is a vocabulary OBJECT (a + # value, like vcfc:ExpandedRepresentation), not a census term. + if i > 0 and seq[i - 1][0] == "V": + continue + terms.add(name) + return terms + + def test_no_emitted_term_is_silently_unmodelled(self): + import re + root = Path(__file__).resolve().parents[1] + emitter = (root / "vcf_rdfizer.py").read_text(encoding="utf-8") + oracle = (root / "src" / "validation" / "validation_runner.py").read_text( + encoding="utf-8" + ) + counted = set(re.findall(r"\{VCFC\}([A-Za-z_][A-Za-z0-9_]*)", oracle)) + counted |= set(re.findall(r"['\"]([A-Za-z_][A-Za-z0-9_]*)['\"]", oracle)) + unmodelled = {t for t in self._emitted_terms(emitter) if t not in counted} + new = sorted(unmodelled - self.KNOWN_UNMODELLED) + self.assertEqual( + new, [], + "the emitter produces vocab terms the oracle never counts, so any file " + "using them fails validation with extra rows. Count them in " + "expected_census, or add them to KNOWN_UNMODELLED with a reason: %s" % new, + ) + + def test_the_waiver_list_does_not_rot(self): + """A waived term that is now counted should leave the list.""" + import re + root = Path(__file__).resolve().parents[1] + oracle = (root / "src" / "validation" / "validation_runner.py").read_text( + encoding="utf-8" + ) + counted = set(re.findall(r"\{VCFC\}([A-Za-z_][A-Za-z0-9_]*)", oracle)) + counted |= set(re.findall(r"['\"]([A-Za-z_][A-Za-z0-9_]*)['\"]", oracle)) + stale = sorted(t for t in self.KNOWN_UNMODELLED if t in counted) + self.assertEqual( + stale, [], + "these are now counted by the oracle and should be removed from " + "KNOWN_UNMODELLED: %s" % stale, + )