From 27d23aa56e3e21e1fef7bbdab139feb11b4ed7ab Mon Sep 17 00:00:00 2001 From: Travis Date: Thu, 13 Aug 2026 20:29:55 -0700 Subject: [PATCH] Make --taxonomy select the priors it names lineage_paths() required One_Hop_Refilled_*.json and Two_Hop_Refilled_*.json under ONE_HOP_NEIGHBOR_DATA and TWO_HOP_NEIGHBOR_DATA. Neither directory is present, and it looks like neither was ever part of the code, so the check failed for every taxonomy including the default and every run loaded domain priors. The two paths were returned and never read, so nothing failed. Four runs of the 0.1.5 container differing only in -t produced byte identical output. Check only the counts table and the neighbour graph, the two files the code reads. The accepted names now come from the installed files rather than a hardcoded list. The lists they replace had drifted three ways: they named "cyanobacteriota" while the file is "cyanobacteriota_melainabacteria_group", omitted "fusobacteriota", and filed "mycoplasmatota" as a kingdom though the -ota suffix makes it a phylum. Two prior sets stayed unreachable even after that check was corrected. A name may be shortened where exactly one canonical name extends it at an underscore, so "cyanobacteriota" and "fcb" resolve while "cyano" does not. An unrecognized name causes an error that lists the valid names by rank, or the near misses when the input is close to one. It used to select domain priors in silence. The domain-level taxon is called bacteria, since every other entry in the list is a taxon rather than a rank. Logging reports the priors in use instead of the string the user typed. Those disagreed: asking for cyanobacteriota logged cyanobacteriota while loading the domain file. tests/taxonomy_test.py runs one input under four taxonomies and requires the results to differ. It fails against the current release. Implemented with assistance from Claude (Opus 5) --- .github/workflows/test.yml | 6 + BLIMMP_Scripts/module_detection.py | 170 +++++++++++++++++++++----- tests/taxonomy_test.py | 184 +++++++++++++++++++++++++++++ 3 files changed, 329 insertions(+), 31 deletions(-) create mode 100644 tests/taxonomy_test.py diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index aa6d1cb..987226d 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -42,3 +42,9 @@ jobs: # 9-genome run. - name: Concurrent extraction run: docker run --rm blimmp:ci python /app/tests/concurrency_test.py + + # --taxonomy was a no-op: every run fell back to domain priors, so four + # runs differing only in -t produced byte identical output. This asserts + # the results differ. + - name: Taxonomy priors + run: docker run --rm blimmp:ci python /app/tests/taxonomy_test.py diff --git a/BLIMMP_Scripts/module_detection.py b/BLIMMP_Scripts/module_detection.py index f423075..96b2748 100644 --- a/BLIMMP_Scripts/module_detection.py +++ b/BLIMMP_Scripts/module_detection.py @@ -27,18 +27,110 @@ ## CONSTANTS -KINGDOM = {"bacillati", "fusobacteriati", "mycoplasmatota", "pseudomonadati", "thermotogati"} -PHYLUM = {"bacillota", "acidobacteriota", "actinomycetota", "campylobacterota", "cyanobacteriota", - "deinococcota", "fcb_group", "mycoplasmatota", "myxococcota", "pseudomonadota", - "pvc_group", "spirochaetota", "thermodesulfobacteriota", "thermotogota"} KO_RE = re.compile(r'^K\d{5}$') +COUNTS_PREFIX = "ko_freq_ko_matrix_sampleids_" +GRAPH_PREFIX = "Module_AllHop_Refilled_" +RANKS = ("domain", "kingdom", "phylum") + + +def available_priors(counts_dir: Path, module_neighbor_dir: Path): + """Map canonical taxon name -> (tag, rank, counts file, neighbour graph file). + + Built from the files present, not from a hardcoded list. The lists this + replaced had drifted from the data three ways: they named "cyanobacteriota" + while the file is "cyanobacteriota_melainabacteria_group", they omitted + "fusobacteriota", and they filed "mycoplasmatota" as a kingdom when the + -ota suffix makes it a phylum. Discovery keeps the accepted names in step + with what is installed. + + Match the filename pattern rather than listing the directory. A README.md + sits in both directories and would otherwise be offered as a taxon. + """ + found = {} + for counts_file in sorted(counts_dir.glob(f"{COUNTS_PREFIX}*.tsv")): + tag = counts_file.name[len(COUNTS_PREFIX):-len(".tsv")] + graph_file = module_neighbor_dir / f"{GRAPH_PREFIX}{tag}.json" + if not graph_file.exists(): + continue + if tag == "domain_level_priors": + # The rank is domain, the taxon is Bacteria. Naming it "domain" + # would put a rank label among taxon names. + name, rank = "bacteria", "domain" + else: + match = re.search(r"_(phylum|kingdom)_level_priors$", tag) + if not match: + continue + rank = match.group(1) + name = tag[:match.start()] + found[name] = (tag, rank, counts_file, graph_file) + return found + + +def describe_priors(available: dict, indent: str = " ") -> str: + """List the taxa grouped by rank. + + Rank separates the confusable pairs: bacillati is a kingdom and bacillota a + phylum, and the same one-letter difference splits fusobacteriati from + fusobacteriota, pseudomonadati from pseudomonadota, and thermotogati from + thermotogota. A flat list sets them side by side with no way to tell which + is which. + """ + lines = [] + for rank in RANKS: + names = sorted(n for n, (_tag, r, _c, _g) in available.items() if r == rank) + if names: + lines.append(f"{indent}{rank + ':':9}{', '.join(names)}") + return "\n".join(lines) + + +def resolve_taxonomy(taxonomy: str, available: dict) -> str: + """Turn what the user typed into one of the canonical names. + + Accepts a canonical name, or a taxon name that exactly one canonical name + extends with a clade qualifier: "cyanobacteriota" reaches + "cyanobacteriota_melainabacteria_group", and "fcb" reaches "fcb_group". + The match must break on an underscore, so "cyano" and "bacillat" are + refused. + + An unrecognized name causes an error. It used to select domain-level + priors silently, so a typo produced plausible output from the wrong priors + with nothing in the log to say so. + """ + val = (taxonomy or "").strip().lower() + if val in ("", "bacteria", "domain"): + return "bacteria" + if val in available: + return val + + hits = sorted(name for name in available if name.startswith(val + "_")) + if len(hits) == 1: + return hits[0] + if len(hits) > 1: + detail = " and ".join(f"{n} ({available[n][1]})" for n in hits) + raise SystemExit( + f"[taxonomy] '{taxonomy}' is ambiguous between {detail}.\n" + f" Give the full name." + ) + # A truncation such as "fusobacteri" is not a taxon, but it is one letter + # from two of them that differ only by rank. Point at those rather than + # reprinting all twenty names. + near = sorted(name for name in available if name.startswith(val)) + if near: + detail = ", ".join(f"{n} ({available[n][1]})" for n in near) + raise SystemExit( + f"[taxonomy] '{taxonomy}' is an unrecognized taxonomic group.\n" + f" Closest: {detail}" + ) + raise SystemExit( + f"[taxonomy] '{taxonomy}' is an unrecognized taxonomic group.\n" + f"{describe_priors(available)}" + ) + ## Configurations @dataclass(frozen=True) class Paths: counts_dir: Path - onehop_dir: Path - twohop_dir: Path module_neighbor_dir: Path module_eq_json: Path module_json_dir: Path @@ -308,24 +400,27 @@ def read_ko_occurrence(kooccpath): @staticmethod def lineage_paths(taxonomy: str, paths: Paths): - val = (taxonomy or "").strip().lower() - if val in PHYLUM: level, name = "phylum", val - elif val in KINGDOM: level, name = "kingdom", val - else: level, name = "domain", "bacteria" - tag = "domain_level_priors" if level == "domain" else f"{name}_{level}_level_priors" - want_counts = paths.counts_dir / f"ko_freq_ko_matrix_sampleids_{tag}.tsv" - want_one = paths.onehop_dir / f"One_Hop_Refilled_{tag}.json" - want_two = paths.twohop_dir / f"Two_Hop_Refilled_{tag}.json" - want_all = paths.module_neighbor_dir / f"Module_AllHop_Refilled_{tag}.json" - if not (want_counts.exists() and want_one.exists() and want_two.exists()): - if tag != "domain_level_priors": - print(f"[taxonomy] Using domain-level fallbacks for '{tag}'.", file=sys.stderr) - tag = "domain_level_priors" - want_counts = paths.counts_dir / f"ko_freq_ko_matrix_sampleids_{tag}.tsv" - want_one = paths.onehop_dir / f"One_Hop_Refilled_{tag}.json" - want_two = paths.twohop_dir / f"Two_Hop_Refilled_{tag}.json" - want_all = paths.module_neighbor_dir / f"Module_AllHop_Refilled_{tag}.json" - return want_counts, want_one, want_two, want_all, tag + """Resolve --taxonomy to the counts table and neighbour graph to use. + + This used to require One_Hop_Refilled_*.json and + Two_Hop_Refilled_*.json under ONE_HOP_NEIGHBOR_DATA and + TWO_HOP_NEIGHBOR_DATA. Neither directory is present, and it looks like + neither was ever part of the code, so the check failed for every + taxonomy including the domain-level one and every run fell back to + domain priors. The two paths were returned and never read, so nothing + failed and --taxonomy did nothing. Check only the files the code + reads. + """ + available = available_priors(paths.counts_dir, paths.module_neighbor_dir) + if not available: + raise SystemExit( + "[taxonomy] No prior sets found. The installed package is missing " + "Data_Dependencies/ATB_Taxonomy_Frequency or " + "Graph_Dependencies/MODULE_ALL_NEIGHBOR_DATA." + ) + name = resolve_taxonomy(taxonomy, available) + tag, _rank, counts_file, graph_file = available[name] + return counts_file, graph_file, tag @staticmethod def modules_to_kos(module_json_dir): @@ -2011,8 +2106,11 @@ def run(self): raise ValueError("--sigma must be between 0 and 1") # taxonomy-driven paths - counts_tsv, onehop_json, twohop_json, all_neighbor_json, tag = File_Helpers.lineage_paths(self.cfg.taxonomy, self.paths) - logging.info(f"Taxonomic level chosen: {self.cfg.taxonomy}") + counts_tsv, all_neighbor_json, tag = File_Helpers.lineage_paths(self.cfg.taxonomy, self.paths) + # Report the priors in use, not the string the user typed. These used + # to disagree: asking for cyanobacteriota logged "cyanobacteriota" + # while loading the domain file. + logging.info(f"Taxonomic priors in use: {tag}") ko_occ = File_Helpers.read_ko_occurrence(str(counts_tsv)) if self.cfg.verbose: @@ -2322,6 +2420,14 @@ def validate_paths(paths: "Paths") -> None: def main(): + # Discovered before the parser is built so --help lists what is installed. + # Two globs, no extraction. + _here = os.path.dirname(os.path.abspath(__file__)) + _priors = available_priors( + Path(_here) / "Data_Dependencies" / "ATB_Taxonomy_Frequency", + Path(_here) / "Graph_Dependencies" / "MODULE_ALL_NEIGHBOR_DATA", + ) + p = argparse.ArgumentParser( description='BLIMMP: Bayesian Likelihood Inference of Metabolic Module Presence. ' 'Evaluates KEGG module completeness from HMM search results.', @@ -2334,9 +2440,13 @@ def main(): p.add_argument('file',help='Path to the HMMER .tblout or .domtblout file') p.add_argument('-f', '--format',choices=['domtblout'], required=True,help='Input file format: "domtblout" for --domtblout') p.add_argument('-s', '--sigma',type=float, required=True,help='Genome completeness estimate (0.0-1.0). Use 1.0 if unknown or for complete genomes') - p.add_argument('-t', '--taxonomy',default="bacteria", metavar="NAME", - help='Taxonomic group for priors (default: bacteria). ' - 'Options include phylum names like "cyanobacteriota" or kingdom names like "pseudomonadati"') + p.add_argument('-t', '--taxonomy', default="bacteria", metavar="NAME", + help='Taxon whose priors to use (default: bacteria, the whole domain). ' + 'Names may be shortened where that is unambiguous, so ' + '"cyanobacteriota" selects cyanobacteriota_melainabacteria_group. ' + 'Available:\n' + + describe_priors(_priors, indent=' ') + + '\nRank suffixes: -ati is a kingdom, -ota a phylum.') p.add_argument('-o', '--output', required=True, metavar="PREFIX", help='Output prefix for result files (e.g., "results/Genomename_Result")') @@ -2391,8 +2501,6 @@ def main(): paths = Paths( counts_dir = Path(DD)/ "ATB_Taxonomy_Frequency", - onehop_dir = Path(GD)/ "ONE_HOP_NEIGHBOR_DATA", - twohop_dir = Path(GD)/ "TWO_HOP_NEIGHBOR_DATA", module_neighbor_dir = Path(GD)/ "MODULE_ALL_NEIGHBOR_DATA", module_eq_json = Path(GD)/ "KEGG_Module_Equations_Jan26.json", module_json_dir = ensure_module_graphs(GD), diff --git a/tests/taxonomy_test.py b/tests/taxonomy_test.py new file mode 100644 index 0000000..07e3510 --- /dev/null +++ b/tests/taxonomy_test.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python +"""--taxonomy must change which priors are used. + +It did not. lineage_paths() required One_Hop_Refilled_*.json and +Two_Hop_Refilled_*.json under two directories that appear never to have been +part of the code, so the check failed for every taxonomy, every run fell back +to domain priors, and four runs differing only in -t produced byte identical +output. Nothing crashed, because the two paths were returned and never read. + +The integration section below would have caught that: run one input under +several taxonomies and require the results to differ. The unit section covers +the name resolution that picks the priors. + +Run inside the container, against the installed package: + docker run --rm blimmp:ci python /app/tests/taxonomy_test.py +""" + +import hashlib +import os +import subprocess +import sys +import tempfile +from pathlib import Path + +from BLIMMP_Scripts.module_detection import ( + available_priors, + describe_priors, + resolve_taxonomy, +) +import BLIMMP_Scripts.module_detection as md + +PKG = Path(md.__file__).parent +COUNTS_DIR = PKG / "Data_Dependencies" / "ATB_Taxonomy_Frequency" +GRAPH_DIR = PKG / "Graph_Dependencies" / "MODULE_ALL_NEIGHBOR_DATA" + +failures = [] + + +def check(condition, message): + if condition: + print(f" ok {message}") + else: + print(f" FAIL {message}") + failures.append(message) + + +def expect_exit(taxonomy, available, why): + try: + got = resolve_taxonomy(taxonomy, available) + except SystemExit: + print(f" ok {taxonomy!r} rejected ({why})") + return + print(f" FAIL {taxonomy!r} resolved to {got!r}, expected rejection ({why})") + failures.append(f"{taxonomy!r} not rejected") + + +print("== discovery ==") +priors = available_priors(COUNTS_DIR, GRAPH_DIR) +check(len(priors) >= 20, f"found {len(priors)} prior sets in the installed package") +check("README.md" not in priors, "README.md is not offered as a taxon") +check("README" not in " ".join(priors), "no README-derived name leaked in") + +by_rank = {} +for name, (_tag, rank, _c, _g) in priors.items(): + by_rank.setdefault(rank, []).append(name) +check(len(by_rank.get("domain", [])) == 1, "exactly one domain-level set") +check(by_rank.get("domain") == ["bacteria"], "the domain-level taxon is named bacteria, not domain") +check(len(by_rank.get("kingdom", [])) == 4, f"4 kingdoms, got {len(by_rank.get('kingdom', []))}") +check(len(by_rank.get("phylum", [])) >= 14, f"{len(by_rank.get('phylum', []))} phyla") + +# Every -ati name is a kingdom and every -ota name a phylum. The list this +# replaced filed mycoplasmatota as a kingdom, which the suffix contradicts. +for name, (_tag, rank, _c, _g) in priors.items(): + if name.endswith("ati"): + check(rank == "kingdom", f"{name} is a kingdom") + elif name.endswith("ota"): + check(rank == "phylum", f"{name} is a phylum") + +print("== every discovered name resolves to itself ==") +bad = [n for n in priors if resolve_taxonomy(n, priors) != n] +check(not bad, f"all {len(priors)} names round-trip" if not bad else f"these did not: {bad}") + +print("== the two sets that used to be unreachable ==") +check("cyanobacteriota_melainabacteria_group" in priors, "cyanobacteria priors are reachable") +check("fusobacteriota" in priors, "fusobacteriota priors are reachable") +check( + resolve_taxonomy("cyanobacteriota", priors) == "cyanobacteriota_melainabacteria_group", + "'cyanobacteriota' reaches the melainabacteria group set", +) + +print("== shortening, only where unambiguous ==") +check(resolve_taxonomy("fcb", priors) == "fcb_group", "'fcb' reaches fcb_group") +check(resolve_taxonomy("pvc", priors) == "pvc_group", "'pvc' reaches pvc_group") +for name in ("", "bacteria", "domain", "BACTERIA"): + check(resolve_taxonomy(name, priors) == "bacteria", f"{name!r} means the whole domain") + +expect_exit("cyano", priors, "an arbitrary truncation") +expect_exit("bacillat", priors, "an arbitrary truncation") +expect_exit("nonsense", priors, "nothing like a taxon name") + +print("== a near miss points at the candidates ==") +# fusobacteri sits one letter from a kingdom and a phylum. Guessing between +# them would pick priors built from a different set of genomes, so it is +# refused, but the message should name the two rather than reprint all twenty. +for stem, expected in ( + ("fusobacteri", ("fusobacteriati", "fusobacteriota")), + ("pseudomonad", ("pseudomonadati", "pseudomonadota")), + ("thermotog", ("thermotogati", "thermotogota")), + ("cyano", ("cyanobacteriota_melainabacteria_group",)), + ("bacillat", ("bacillati",)), +): + try: + resolve_taxonomy(stem, priors) + check(False, f"{stem!r} should have been refused") + except SystemExit as exc: + message = str(exc) + check(all(name in message for name in expected), + f"{stem!r} is refused and named {', '.join(expected)}") + check(any(f"({rank})" in message for rank in ("kingdom", "phylum")), + f"{stem!r} suggestion carries the rank") + +print("== the listing names the rank ==") +described = describe_priors(priors) +check("kingdom:" in described and "phylum:" in described, "output is grouped by rank") +check("bacillati" in described and "bacillota" in described, "both confusable names are listed") + +print("== end to end: different taxonomy, different result ==") +example = Path("/app/Examples/example.domtblout") +if not example.exists(): + example = PKG.parent / "Examples" / "example.domtblout" + +if not example.exists(): + print(f" SKIP no example domtblout at {example}") +else: + digests = {} + with tempfile.TemporaryDirectory() as tmp: + for taxonomy in ("bacteria", "cyanobacteriota", "actinomycetota", "bacillati"): + prefix = os.path.join(tmp, taxonomy) + result = subprocess.run( + ["BLIMMP", str(example), "-f", "domtblout", "--sigma", "1.0", + "-t", taxonomy, "-o", prefix], + capture_output=True, text=True, + ) + if result.returncode != 0: + print(f" FAIL -t {taxonomy} exited {result.returncode}") + print(result.stderr[-2000:]) + failures.append(f"-t {taxonomy} failed") + continue + out = Path(f"{prefix}_BLIMMP_module_probabilities.csv") + if not out.exists(): + print(f" FAIL -t {taxonomy} wrote no probabilities file") + failures.append(f"-t {taxonomy} produced nothing") + continue + digests[taxonomy] = hashlib.md5(out.read_bytes()).hexdigest() + + for taxonomy, digest in digests.items(): + print(f" {taxonomy:20s} {digest}") + check( + len(set(digests.values())) == len(digests), + f"all {len(digests)} taxonomies gave distinct results " + f"(got {len(set(digests.values()))} distinct)", + ) + +print("== an unknown taxonomy stops the run ==") +if example.exists(): + with tempfile.TemporaryDirectory() as tmp: + result = subprocess.run( + ["BLIMMP", str(example), "-f", "domtblout", "--sigma", "1.0", + "-t", "not_a_real_taxon", "-o", os.path.join(tmp, "x")], + capture_output=True, text=True, + ) + check(result.returncode != 0, "an invalid --taxonomy exits non-zero") + combined = result.stdout + result.stderr + check("unrecognized taxonomic group" in combined, + "the error calls the name unrecognized") + check("kingdom:" in combined, "the error lists the valid names by rank") + +if failures: + print(f"\nFAILED: {len(failures)}") + for item in failures: + print(f" - {item}") + sys.exit(1) + +print("\nAll taxonomy checks passed.")