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, + )