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.")