diff --git a/scripts/gtdb_ground.py b/scripts/gtdb_ground.py index fa2f646a..bddae2cf 100644 --- a/scripts/gtdb_ground.py +++ b/scripts/gtdb_ground.py @@ -157,6 +157,58 @@ def resolve_kg_microbe_dir(explicit: str | None) -> Path: ) +# GTDB appends an alphabetic suffix to every clade split off a genus/family that +# does NOT contain the nomenclatural type: `g__Enterococcus_B` is Enterococcus +# minus the lineage holding *E. faecalis*, which keeps the bare `g__Enterococcus`. +# So the suffix is GTDB's own published marker for "not the type-anchored clade". +_NON_TYPE_SUFFIX = re.compile(r"^(?P[A-Z][A-Za-z0-9-]*)_(?P[A-Z]+)$") + + +def non_type_clade(name: str) -> str | None: + """The unsuffixed sibling of a GTDB name, when the name is a non-type clade. + + Returns the base name (`Enterococcus` for `Enterococcus_B`), else None. + + Used to *flag* rather than to *resolve* (#374). Grounding by majority-of- + genomes and GTDB's type-species rule disagree whenever a non-type clade is + more heavily sequenced. Verified against GTDB R226, that is one taxon of the + five where the two denominators differ — `Enterococcus` grounds to + `g__Enterococcus_B` while *E. faecalis* sits in `g__Enterococcus`. The other + four agree because the majority happens to be the type clade, which is luck + and not construction; nothing stops them flipping as sequencing depth shifts. + + Deliberately not used to override the answer. Preferring the type clade + would change what every existing grounding means, on the evidence of a + single live case — and the majority answer is not obviously wrong, since a + user asking for "NCBI genus X" may well want the clade most of X's genomes + are in. Emitting both and letting a curator see the conflict costs nothing + and forecloses nothing. + """ + match = _NON_TYPE_SUFFIX.match(name or "") + return match.group("base") if match else None + + +def _non_type_warning(chosen: str, weights: dict) -> str: + """Human-readable note for a grounding that is not the type-anchored clade. + + Says whether the unsuffixed sibling was even a contender, because the two + cases call for different curator action: if it drew genomes here, this is a + close call worth a look; if it drew none, the NCBI name's genomes simply do + not sit in the type clade at all, which is a stronger statement. + """ + base = non_type_clade(chosen) + sibling = weights.get(base) + if sibling: + return ( + f"{chosen} is not the type-anchored clade; GTDB reserves {base} for the " + f"lineage holding the nomenclatural type, which drew {int(sibling)} genomes here (#374)" + ) + return ( + f"{chosen} is not the type-anchored clade; GTDB reserves {base} for the " + f"lineage holding the nomenclatural type, which drew no genomes here (#374)" + ) + + def _curie(name: str, prefix: str) -> str: return f"GTDB:{prefix}__" + name.replace(" ", "_") @@ -584,6 +636,15 @@ def resolve_higher( "is_reclassified": top != _clean_label(label), "via": f"ncbi_rank_{prefix}", "n_alt": len(weights), + # The majority answer is a non-type clade, so it disagrees with + # GTDB's own naming rule (#374). Reported, not resolved: the + # curator decides. Key omitted entirely when there is no + # conflict, so its presence always means something. + **( + {"non_type_clade_warning": _non_type_warning(top, weights)} + if non_type_clade(top) + else {} + ), } # Same tie-break as `top` above: ties here were still row-ordered, so the # AMBIGUOUS option list a curator reads was not reproducible (#382 review). @@ -1913,6 +1974,12 @@ def main(argv: list[str] | None = None) -> int: # record's *nitrite* oxidizer. The number alone did not say "look here". near = " ⚠ NEAR-TIE" if _is_near_tie(g.get("majority_fraction")) else "" print(f" majority : {g['majority_fraction']}{via}{of}{thin}{near}") + # Printed, not just stored (#374). The whole value of flagging rather + # than resolving is that a curator sees the conflict at the moment they + # ground the taxon; a key that only ever appears in emitted YAML would + # be found later or not at all. + if g.get("non_type_clade_warning"): + print(f" ⚠ NON-TYPE : {g['non_type_clade_warning']}") if args.emit_yaml: print(" --- gtdb_classification block ---") for line in emit_block(g, mapping_source).splitlines(): diff --git a/tests/test_gtdb_non_type_clade_flag.py b/tests/test_gtdb_non_type_clade_flag.py new file mode 100644 index 00000000..0f6a1bc0 --- /dev/null +++ b/tests/test_gtdb_non_type_clade_flag.py @@ -0,0 +1,156 @@ +"""Grounding to a non-type clade is reported, not resolved (#374). + +`gtdb_ground.py` grounds a higher-rank NCBI taxon to whichever GTDB taxon holds +the most genomes. GTDB names by **nomenclatural type**: the lineage containing +the type species keeps the unsuffixed name and the rest take alphabetic +suffixes. The two rules disagree whenever a non-type clade is more heavily +sequenced. + +Verified against GTDB R226 (the mapping this repo reads), `Enterococcus` grounds +to `g__Enterococcus_B` at 0.598 while *E. faecalis* — the type species — sits in +`g__Enterococcus`, which drew 9904 of the 28462 genomes. The other four taxa +where the denominators differ agree only because the majority happens to be the +type clade, which is luck rather than construction. + +**Reported, not resolved, by decision.** Preferring the type clade would change +what every existing grounding means on the evidence of one live case, and the +majority answer is not obviously wrong — someone asking for "NCBI genus X" may +well want the clade most of X's genomes are in. So the tool emits both and the +curator decides. + +The suffix is not a heuristic: it is GTDB's own published marker for "this clade +does not contain the type". That is why the check can be a regex and still be +correct. +""" + +from __future__ import annotations + +import importlib.util +import pathlib + +import pytest + +REPO = pathlib.Path(__file__).parent.parent +SCRIPT = REPO / "scripts/gtdb_ground.py" + + +@pytest.fixture(scope="module") +def gtdb(): + spec = importlib.util.spec_from_file_location("gtdb_ground_under_test", SCRIPT) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +@pytest.mark.parametrize( + ("name", "expected_base"), + [ + # Every non-type clade actually present in the KB today. + ("Enterococcus_B", "Enterococcus"), + ("Cetobacterium_A", "Cetobacterium"), + ("Methanobrevibacter_A", "Methanobrevibacter"), + ("Bacillota_A", "Bacillota"), + ("Clostridium_AM", "Clostridium"), + ("Acidithiobacillus_A", "Acidithiobacillus"), + ("Pseudomonas_E", "Pseudomonas"), + ("Ruminococcus_B", "Ruminococcus"), + ], +) +def test_a_suffixed_name_is_a_non_type_clade(gtdb, name, expected_base): + assert gtdb.non_type_clade(name) == expected_base + + +@pytest.mark.parametrize( + "name", + [ + # The type-anchored names themselves. + "Enterococcus", + "Pseudomonas", + "Bacillus", + "Leptospirillum", + # A binomial is not a genus name; the species path does not get this + # warning, and passing one here must not produce a false positive. + "Ruminococcus_B gnavus", + "Clostridium_AM drakei", + # Genuine names that merely contain characters the pattern touches. + "Candidatus Nitrosocosmicus", + "CAG-267", + "", + ], +) +def test_these_are_not_flagged(gtdb, name): + assert gtdb.non_type_clade(name) is None + + +def test_the_warning_names_the_type_clade_and_its_support(gtdb): + """The message has to tell a curator what to compare against. + + Numbers from the real Enterococcus case: the type clade is not a fringe + option, it drew 9904 genomes against the winner's 17011, which is exactly + the situation where a curator's judgement is worth having. + """ + warning = gtdb._non_type_warning( + "Enterococcus_B", {"Enterococcus_B": 17011, "Enterococcus": 9904} + ) + + assert "Enterococcus_B is not the type-anchored clade" in warning + assert "reserves Enterococcus for" in warning + assert "9904 genomes" in warning + assert "#374" in warning + + +def test_the_warning_distinguishes_a_type_clade_with_no_genomes(gtdb): + """Different curator action, so a different sentence. + + If the type clade drew nothing, the NCBI name's genomes simply do not sit in + it — a stronger statement than a close call, and it should not read as one. + """ + warning = gtdb._non_type_warning("Cetobacterium_A", {"Cetobacterium_A": 40}) + + assert "drew no genomes here" in warning + assert "0 genomes" not in warning, "a bare 0 reads like a close call that lost" + + +def _ground(name: str) -> str: + """Run the real CLI against the real GTDB mapping.""" + import subprocess + + result = subprocess.run( + ["uv", "run", "python", str(SCRIPT), "--name", name], + capture_output=True, + text=True, + cwd=REPO, + ) + return result.stdout + + +@pytest.mark.integration +def test_enterococcus_warns_and_still_returns_the_majority_answer(): + """End-to-end on the one taxon in the KB where the rules actually disagree. + + Flagging must not change the answer — that was the decision. If someone + later makes this prefer the type clade, the second assertion goes red and + they have to say so deliberately. + """ + out = _ground("Enterococcus") + if "GTDB taxon" not in out: + pytest.skip("local GTDB mapping unavailable") + + assert "NON-TYPE" in out, "the one real conflict in the KB produced no warning" + assert "g__Enterococcus_B" in out, "the majority answer was silently overridden" + + +@pytest.mark.integration +@pytest.mark.parametrize("name", ["Pseudomonas", "Acetobacter", "Leptospirillum", "Bacillus"]) +def test_the_four_that_agree_do_not_warn(name): + """The other four taxa where the denominators differ. + + They agree with the type rule only because the majority happens to be the + type clade. A warning here would be a false positive on 4 of 5 cases, which + would train a curator to ignore it. + """ + out = _ground(name) + if "GTDB taxon" not in out: + pytest.skip("local GTDB mapping unavailable") + + assert "NON-TYPE" not in out, f"{name} agrees with the type rule but was flagged"