Skip to content

Commit 20d2cbb

Browse files
ecrum19claude
andcommitted
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 <noreply@anthropic.com>
1 parent a3679e1 commit 20d2cbb

2 files changed

Lines changed: 221 additions & 0 deletions

File tree

‎src/validation/validation_runner.py‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -583,6 +583,11 @@ def emitted_record_counters(
583583
assembly_contig_ids: set[str] = set()
584584
reference_alleles = alt_alleles = 0
585585
value_items = value_item_alleles = tuple_items = 0
586+
# FORMAT value items are counted separately: the emitter decomposes them
587+
# only in the expanded per-sample representation, so expected_census merges
588+
# them for that profile alone, exactly as it does the genotype layer.
589+
format_items = format_item_alleles = format_tuple_items = 0
590+
format_item_predicates: Counter[str] = Counter()
586591
genotypes = genotype_calls = called_alleles = 0
587592

588593
for row in rows:
@@ -650,6 +655,43 @@ def emitted_record_counters(
650655
if samples:
651656
format_keys = (row[8].split(":") if len(row) > 8 and row[8] else [])
652657
payloads = row[9 : 9 + len(samples)]
658+
# A positional FORMAT cell is decomposed into vcfc:FieldValueItems
659+
# exactly as a positional INFO one is -- _emit_value_items is called
660+
# for both, from append_expanded_sample_rdf and from the INFO pass.
661+
# Counting only the INFO side made every file with a positional
662+
# FORMAT key (AD, ADALL, PL -- i.e. most real VCFs) fail validation
663+
# with the whole item layer reported as unexpected extra rows.
664+
#
665+
# The emitter requires at least one ALT before it decomposes, and
666+
# split_value_items returns nothing for a missing cell, so both
667+
# conditions are mirrored rather than re-derived.
668+
if alt_count and format_keys:
669+
for payload in payloads:
670+
fields = payload.split(":") if payload else []
671+
for key_index, key in enumerate(format_keys):
672+
cell = fields[key_index] if key_index < len(fields) else ""
673+
if not cell:
674+
continue
675+
number = format_numbers.get(key, ".")
676+
if not version.is_positional(key, number):
677+
continue
678+
items = vocab.split_value_items(cell)
679+
if not items:
680+
continue
681+
format_items += len(items)
682+
if version.tuple_arity(key) is not None:
683+
format_tuple_items += len(items)
684+
for index in range(len(items)):
685+
link = version.value_item_link(key, number, index)
686+
if (
687+
link.allele_index is not None
688+
and link.allele_index in allele_uris
689+
):
690+
format_item_alleles += 1
691+
elif link.genotype_index is not None:
692+
format_item_predicates["forGenotypeIndex"] += 1
693+
elif link.gt_allele_index is not None:
694+
format_item_predicates["forGTAlleleIndex"] += 1
653695
if "GT" in format_keys:
654696
gt_index = format_keys.index("GT")
655697
for payload in payloads:
@@ -696,11 +738,22 @@ def emitted_record_counters(
696738
genotype_predicates[name] += genotype_calls
697739
genotype_predicates["calledAllele"] += called_alleles
698740

741+
format_item_classes: Counter[str] = Counter()
742+
if format_items:
743+
format_item_classes["FieldValueItem"] += format_items
744+
for name in ("hasValueItem", "valueIndex", "itemValue"):
745+
format_item_predicates[name] += format_items
746+
format_item_predicates["forAllele"] += format_item_alleles
747+
format_item_predicates["tupleArity"] += format_tuple_items
748+
699749
return {
700750
"emittedRecordClasses": dict(classes),
701751
"emittedRecordPredicates": dict(predicates),
702752
"emittedGenotypeClasses": dict(genotype_classes),
703753
"emittedGenotypePredicates": dict(genotype_predicates),
754+
"emittedFormatItemClasses": dict(format_item_classes),
755+
"emittedFormatItemPredicates": dict(format_item_predicates),
756+
"formatValueItemCount": format_items,
704757
"alleleCount": total_alleles,
705758
"valueItemCount": value_items,
706759
"genotypeCount": genotypes,
@@ -870,6 +923,12 @@ def expected_census(
870923
classes[f"{VCFC}{class_name}"] = classes.get(f"{VCFC}{class_name}", 0) + count
871924
for name, count in parser.get("emittedGenotypePredicates", {}).items():
872925
predicates[f"{VCFC}{name}"] = predicates.get(f"{VCFC}{name}", 0) + count
926+
# The FORMAT item layer is emitted by append_expanded_sample_rdf, so
927+
# it exists only here -- the condensed profile never materializes it.
928+
for class_name, count in parser.get("emittedFormatItemClasses", {}).items():
929+
classes[f"{VCFC}{class_name}"] = classes.get(f"{VCFC}{class_name}", 0) + count
930+
for name, count in parser.get("emittedFormatItemPredicates", {}).items():
931+
predicates[f"{VCFC}{name}"] = predicates.get(f"{VCFC}{name}", 0) + count
873932
classes[f"{VCFC}SampleCall"] = records * samples
874933
classes[f"{VCFC}FormatFieldValue"] = parser["formatValueSlots"]
875934
predicates[f"{VCFC}hasSampleCall"] = records * samples

‎test/test_validation_oracle_unit.py‎

Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -788,3 +788,165 @@ def test_failure_codes_are_counted_and_matched_against_declarations(self):
788788

789789
def test_a_truncated_record_is_treated_as_missing_rather_than_crashing(self):
790790
self.assertEqual(V.record_decomposition_counts(["chr1", "1"], set()), (0, 0, 0))
791+
792+
793+
class FormatValueItemCensusTests(VerboseTestCase):
794+
"""The oracle must count FORMAT value items, not only INFO ones.
795+
796+
Regression: _emit_value_items is called from two places -- the INFO pass and
797+
append_expanded_sample_rdf -- but emitted_record_counters only walked the
798+
INFO entries. format_numbers was even passed in and never read. Any file with
799+
a positional FORMAT key (AD, ADALL, PL -- most real VCFs) then validated with
800+
the whole item layer reported as unexpected extra rows: 401,606 of them on a
801+
100k-record GIAB file, which failed 3/3 replicates of 13_query_cost.
802+
"""
803+
804+
def _counters(self, rows, samples=("S1",), format_numbers=None, info_numbers=None):
805+
import sys
806+
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
807+
import vcf_rdfizer_vocab as vocab
808+
version, _ = vocab.resolve_vcf_version("4.2")
809+
return V.emitted_record_counters(
810+
[list(r) for r in rows],
811+
list(samples),
812+
version=version,
813+
contig_ids={"chr1"},
814+
alt_declaration_ids=set(),
815+
info_numbers=info_numbers or {"DP": "1"},
816+
format_numbers=format_numbers or {"AD": "R", "GT": "1"},
817+
)
818+
819+
ROW = ["chr1", "100", ".", "A", "G", "50", "PASS", "DP=10", "GT:AD", "0/1:5,7"]
820+
821+
def test_positional_format_cell_is_counted(self):
822+
out = self._counters([self.ROW])
823+
self.assertEqual(out["emittedFormatItemClasses"].get("FieldValueItem"), 2)
824+
for name in ("hasValueItem", "valueIndex", "itemValue"):
825+
self.assertEqual(out["emittedFormatItemPredicates"].get(name), 2, name)
826+
# Number=R indexes ref then each ALT, so both items name an allele.
827+
self.assertEqual(out["emittedFormatItemPredicates"].get("forAllele"), 2)
828+
829+
def test_format_items_are_kept_apart_from_info_items(self):
830+
"""They merge under different profiles, so they cannot share a counter."""
831+
out = self._counters([self.ROW])
832+
self.assertEqual(out["valueItemCount"], 0)
833+
self.assertEqual(out["formatValueItemCount"], 2)
834+
835+
def test_non_positional_format_key_contributes_nothing(self):
836+
out = self._counters([self.ROW], format_numbers={"AD": "1", "GT": "1"})
837+
self.assertFalse(out["emittedFormatItemClasses"])
838+
839+
def test_missing_cell_contributes_nothing(self):
840+
row = list(self.ROW); row[9] = "0/1:."
841+
self.assertFalse(self._counters([row])["emittedFormatItemClasses"])
842+
843+
def test_no_alt_contributes_nothing(self):
844+
"""The emitter requires at least one ALT before it decomposes."""
845+
row = ["chr1", "100", ".", "A", ".", "50", "PASS", "DP=10", "GT:AD", "0/0:5"]
846+
self.assertFalse(self._counters([row])["emittedFormatItemClasses"])
847+
848+
def test_the_layer_is_expanded_only(self):
849+
"""append_expanded_sample_rdf emits it; the condensed profile does not."""
850+
parser = self._counters([self.ROW])
851+
self.assertTrue(parser["emittedFormatItemClasses"])
852+
853+
854+
class EmittedTermCensusCoverageTests(VerboseTestCase):
855+
"""Every vocab term the emitter produces must be counted, or waived here.
856+
857+
This is the guard for a whole class of bug: the emitter grows a predicate or
858+
class, the oracle never learns to expect it, and every file using that
859+
feature fails validation with extra rows. It has happened twice -- the INFO
860+
value items under raw INFO, and the FORMAT value items.
861+
862+
A term in KNOWN_UNMODELLED is a deliberate, documented gap. Shrinking the set
863+
is progress. Growing it needs a reason in the commit message, because each
864+
entry is a file shape that cannot pass validation today.
865+
"""
866+
867+
# Grouped by the feature that produces them.
868+
KNOWN_UNMODELLED = {
869+
# VCF 4.5 local alleles (LA/LR/LG)
870+
"LocalAlleleSet", "LocalAlleleMembership", "hasLocalAlleleSet",
871+
"hasLocalAlleleMembership", "hasLocalAllele", "localAllele", "localIndex",
872+
# phase sets (PS/PSL/PSO/PSQ)
873+
"PhaseSet", "inPhaseSet", "phaseSetId", "phaseSetName", "phaseSetOrdinal",
874+
"phaseSetQuality",
875+
# structural variants and confidence intervals
876+
"VariantEvent", "inEvent", "eventType", "svClaim",
877+
"ConfidenceInterval", "ciLower", "ciUpper", "endPosition",
878+
# gVCF reference blocks
879+
"ReferenceBlock", "isReferenceBlockStart", "referenceBlockLength",
880+
# tandem repeats
881+
"TandemRepeatAllele", "RepeatSequence", "hasRepeatSequence",
882+
"repeatSequenceCount", "repeatSequenceIndex",
883+
# Number=M base modifications
884+
"BaseModification", "forBaseModification", "modifiedBaseOffset",
885+
"modifiedResidue",
886+
# raw carriers, deliberately not part of the census
887+
"sampleDataRaw", "sampleFilter",
888+
}
889+
890+
@staticmethod
891+
def _emitted_terms(source: str):
892+
"""Classify each _vocab(X) use by its position in the triple."""
893+
import re
894+
calls, buf, depth = [], [], 0
895+
for line in source.split("\n"):
896+
if "emit(" in line and depth == 0:
897+
buf, depth = [line], line.count("(") - line.count(")")
898+
if depth <= 0:
899+
calls.append(" ".join(buf)); buf, depth = [], 0
900+
elif depth > 0:
901+
buf.append(line); depth += line.count("(") - line.count(")")
902+
if depth <= 0:
903+
calls.append(" ".join(buf)); buf, depth = [], 0
904+
terms = set()
905+
for call in calls:
906+
raw = re.findall(
907+
r"(RDF_TYPE_URI)|_vocab\(['\"]([A-Za-z_][A-Za-z0-9_]*)['\"]\)", call
908+
)
909+
seq = [("TYPE", None) if a else ("V", b) for a, b in raw]
910+
for i, (kind, name) in enumerate(seq):
911+
if kind != "V":
912+
continue
913+
# A term preceded by another term is a vocabulary OBJECT (a
914+
# value, like vcfc:ExpandedRepresentation), not a census term.
915+
if i > 0 and seq[i - 1][0] == "V":
916+
continue
917+
terms.add(name)
918+
return terms
919+
920+
def test_no_emitted_term_is_silently_unmodelled(self):
921+
import re
922+
root = Path(__file__).resolve().parents[1]
923+
emitter = (root / "vcf_rdfizer.py").read_text(encoding="utf-8")
924+
oracle = (root / "src" / "validation" / "validation_runner.py").read_text(
925+
encoding="utf-8"
926+
)
927+
counted = set(re.findall(r"\{VCFC\}([A-Za-z_][A-Za-z0-9_]*)", oracle))
928+
counted |= set(re.findall(r"['\"]([A-Za-z_][A-Za-z0-9_]*)['\"]", oracle))
929+
unmodelled = {t for t in self._emitted_terms(emitter) if t not in counted}
930+
new = sorted(unmodelled - self.KNOWN_UNMODELLED)
931+
self.assertEqual(
932+
new, [],
933+
"the emitter produces vocab terms the oracle never counts, so any file "
934+
"using them fails validation with extra rows. Count them in "
935+
"expected_census, or add them to KNOWN_UNMODELLED with a reason: %s" % new,
936+
)
937+
938+
def test_the_waiver_list_does_not_rot(self):
939+
"""A waived term that is now counted should leave the list."""
940+
import re
941+
root = Path(__file__).resolve().parents[1]
942+
oracle = (root / "src" / "validation" / "validation_runner.py").read_text(
943+
encoding="utf-8"
944+
)
945+
counted = set(re.findall(r"\{VCFC\}([A-Za-z_][A-Za-z0-9_]*)", oracle))
946+
counted |= set(re.findall(r"['\"]([A-Za-z_][A-Za-z0-9_]*)['\"]", oracle))
947+
stale = sorted(t for t in self.KNOWN_UNMODELLED if t in counted)
948+
self.assertEqual(
949+
stale, [],
950+
"these are now counted by the oracle and should be removed from "
951+
"KNOWN_UNMODELLED: %s" % stale,
952+
)

0 commit comments

Comments
 (0)