From 42eb0e06fd96c7cbd4c896e324c6b0c8a72de7b2 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Wed, 4 Nov 2026 16:21:25 +0100 Subject: [PATCH 1/5] fix: small internal bugs in HOG class, UnionFind, batching and helpers - UnionFind: take over the fixed zoo version. remove() of a set's root or of an element others still point to no longer breaks find(), and union() counts the weight of a shared root only once. Drop the identical copy of the class in _utils_roothog.py, which shadowed the import. - Representative(rep, elements) copies rep's sets instead of extending them in place; find_representative() uses get_subelements(); HOG() without representatives for a merged HOG raises a clear ValueError (the default was a typing object). - batch_roothogs rejects rootHOG files with the same name in different subfolders instead of letting them overwrite each other. - pairwise relations: an empty orthologGroup no longer crashes set.union(). - prepare_species_tree returns the species leaf for a single-species rootHOG (no output change: such rootHOGs still produce no group). - copy_tree: drop the dead feature-collecting loop (its set.union() result was discarded) and document what the newick copy keeps; same output as before. Fixes #122 Co-Authored-By: Claude Opus 5.5 --- FastOMA/_hog_class.py | 11 ++- FastOMA/_infer_subhog.py | 17 +---- FastOMA/_utils_roothog.py | 124 ------------------------------- FastOMA/_utils_subhog.py | 8 +- FastOMA/batch_roothogs.py | 7 ++ FastOMA/zoo/hog/transform.py | 2 +- FastOMA/zoo/unionfind.py | 24 +++++- tests/test_batch_roothogs.py | 10 +-- tests/test_hog_class.py | 16 ++-- tests/test_pairwise_relations.py | 2 - tests/test_species_tree_prep.py | 2 - tests/test_unionfind.py | 61 +++++++++++++-- 12 files changed, 110 insertions(+), 174 deletions(-) diff --git a/FastOMA/_hog_class.py b/FastOMA/_hog_class.py index 9e80586..ced8d8c 100644 --- a/FastOMA/_hog_class.py +++ b/FastOMA/_hog_class.py @@ -246,8 +246,9 @@ def __init__(self, the_one, elements:Optional[List]=None): self.enabled = True if isinstance(the_one, Representative): self._seq = the_one.get_record() - self._subelements = the_one.get_subelements() - self._species = the_one.get_species() + # copies: extending them below must not change the_one + self._subelements = set(the_one.get_subelements()) + self._species = set(the_one.get_species()) elif isinstance(the_one, SeqRecord): self._seq = the_one self._subelements = set([the_one.id]) @@ -277,7 +278,7 @@ class HOG: _hogid_iter = 10000 def __init__(self, input_instantiate, taxnomic_range:TreeNode, rhogid:str, msa:Optional[MultipleSeqAlignment] = None, - representatives=Optional[List[Representative]], conf_infer_subhhogs=None): + representatives: Optional[List[Representative]] = None, conf_infer_subhhogs=None): """HOG class constructor. A HOG can be initialized with either a single SeqRecord or from a list of HOGs @@ -316,6 +317,8 @@ def __init__(self, input_instantiate, taxnomic_range:TreeNode, rhogid:str, msa:O hog_members |= sub_hog.get_members() self._members = hog_members # set.union(*tup) # this also include dubious_members self._subhogs = list(input_instantiate) # full members of subhog, children + if representatives is None: + raise ValueError("a HOG merged from sub-hogs needs its representatives") self._representatives = [r for r in representatives] representatives_id = set(r.get_id() for r in representatives) assert len(representatives_id) == len(self._representatives), "Non-unique ids among representatives" @@ -380,7 +383,7 @@ def find_representative(self, seq_id): if r.get_id() == seq_id: return r for r in self._representatives: - if seq_id in r.get_members(): + if seq_id in r.get_subelements(): return r raise KeyError("No representative found for " + seq_id) diff --git a/FastOMA/_infer_subhog.py b/FastOMA/_infer_subhog.py index 8a2460b..8af9b82 100644 --- a/FastOMA/_infer_subhog.py +++ b/FastOMA/_infer_subhog.py @@ -773,18 +773,9 @@ def infer_hogs_this_level(node_species_tree, rhogid, pickles_subhog_folder_all, def copy_tree(t:TreeNode): - """This function creates a deep copy of an ete3 tree while - preserving the names - - standard code in ete3: - new_node = self.__class__(self.write(features=[])) - """ - features = set([]) - for n in t.traverse(): - features.union(n.features) - features -= {"name"} - new_node = t.__class__(t.write(features=features, quoted_node_names=True, format_root_node=True), - quoted_node_names=True) - return new_node + """Copy of an ete3 tree via newick: topology, (quoted) leaf names, branch lengths and support + values. Internal node names and other node features are not copied (newick format 0 writes the + support value in place of an internal node's name).""" + return t.__class__(t.write(quoted_node_names=True, format_root_node=True), quoted_node_names=True) diff --git a/FastOMA/_utils_roothog.py b/FastOMA/_utils_roothog.py index f11dc65..e87571d 100644 --- a/FastOMA/_utils_roothog.py +++ b/FastOMA/_utils_roothog.py @@ -24,130 +24,6 @@ Gene = collections.namedtuple("Gene", ("numeric_id", "prot_id", "main_isoform"), defaults=(None,) ) -"""UnionFind.py - -Union-find data structure. Based on Josiah Carlson's code, -http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/215912 -with significant additional changes by D. Eppstein and -Adrian Altenhoff. -""" - - -class UnionFind(object): - """Union-find data structure. - - Each unionFind instance X maintains a family of disjoint sets of - hashable objects, supporting the following two methods: - - - X[item] returns a name for the set containing the given item. - Each set is named by an arbitrarily-chosen one of its members; as - long as the set remains unchanged it will keep the same name. If - the item is not yet part of a set in X, a new singleton set is - created for it. - - - X.union(item1, item2, ...) merges the sets containing each item - into a single larger set. If any item is not yet part of a set - in X, it is added to X as one of the members of the merged set. - """ - - def __init__(self, elements=None): - """Create a new union-find structure. - - If elements is not None, the structure gets initialized - with each element as a singleton component. - - :param elements: an iterable to initialize the structure. - """ - - self.weights = {} - self.parents = {} - if elements is not None: - for elem in iter(elements): - self.parents[elem] = elem - self.weights[elem] = 1 - - def __getitem__(self, obj): - """return the name of set which contains obj. - - :param obj: the query object - - :SeeAlso: :meth:`find`""" - return self.find(obj) - - def find(self, obj): - """Find and return the name of the set containing the obj. - - If the object is not found in any set, a new singleton set - is created that holds only this object until it is further merged.""" - - # check for previously unknown obj. If unknown, add it - # as a new cluster - if obj not in self.parents: - self.parents[obj] = obj - self.weights[obj] = 1 - return obj - - # find path of objects leading to the root - path = [obj] - root = self.parents[obj] - while root != path[-1]: - path.append(root) - root = self.parents[root] - - # compress the path and return - for ancestor in path: - self.parents[ancestor] = root - return root - - def remove(self, obj): - """Remove an object from the sets. - - Removes an object entirly from the datastructure. The - containing set will shrink by this one element. - - :Note: If one tries to accessed it afterwards using - :meth:`find`, it will be created newly and put as a - singleton. - """ - if obj not in self.parents: - return - comp = self.find(obj) - self.weights[comp] -= 1 - self.parents.pop(obj) - - def __iter__(self): - """Iterate through all items ever found or unioned by this structure.""" - return iter(self.parents) - - def union(self, *objects): - """Find the sets containing the objects and merge them. - - any number of objects can be passed to this method and - all of them will be merged into one set containing at - least these objects. - - :param objects: the objects to be merged. they have to be all - hashable. If they haven't been initialy added to the UnionFind - datastructre at instantiation time, they are added at this point - in time. - """ - roots = [self[x] for x in objects] - heaviest = max([(self.weights[r], r) for r in roots], key=lambda x: x[0])[1] - for r in roots: - if r != heaviest: - self.weights[heaviest] += self.weights[r] - self.parents[r] = heaviest - - def get_components(self): - """return a list of sets corresponding to the connected - components of the structure.""" - comp_dict = collections.defaultdict(set) - for elem in iter(self): - comp_dict[self[elem]].add(elem) - comp = list(comp_dict.values()) - return comp - - def parse_proteomes(folder=None, min_sequence_length=0): # list_oma_species """ parsing fasta files of proteins located in /proteome/ diff --git a/FastOMA/_utils_subhog.py b/FastOMA/_utils_subhog.py index 2fa4f07..cd48bcd 100644 --- a/FastOMA/_utils_subhog.py +++ b/FastOMA/_utils_subhog.py @@ -280,8 +280,12 @@ def prepare_species_tree(rhog_i: List[SeqRecord], species_tree: Tree, rhogid: st full_species_tree = species_tree.copy() - mrca = species_tree.get_common_ancestor(species_names_uniqe) - mrca.prune(species_names_uniqe, preserve_branch_length=True) + if len(species_names_uniqe) == 1: + # get_common_ancestor of a single name returns the tree root, not the species itself + mrca = species_tree.search_nodes(name=next(iter(species_names_uniqe)))[0] + else: + mrca = species_tree.get_common_ancestor(species_names_uniqe) + mrca.prune(species_names_uniqe, preserve_branch_length=True) return mrca, species_names_rhog, prot_names_rhog, full_species_tree diff --git a/FastOMA/batch_roothogs.py b/FastOMA/batch_roothogs.py index 332446b..e380ebd 100644 --- a/FastOMA/batch_roothogs.py +++ b/FastOMA/batch_roothogs.py @@ -1,4 +1,5 @@ +import collections import shutil from pathlib import Path from . import __version__ as fastoma_version, logger @@ -52,6 +53,12 @@ def _flush(self): def folder_1h_rhog(roothog_path: Path, output_folder_big: Path, output_folder_rest: Path): # create a list of hogs in descending filesize order hog_size_tuples = sorted([(f, f.stat().st_size) for f in roothog_path.rglob("*.fa")], key=lambda x: -x[1]) + # the file name is the rootHOG id downstream, and batches are flat folders + names = collections.Counter(f.name for f, _ in hog_size_tuples) + duplicates = sorted(name for name, cnt in names.items() if cnt > 1) + if duplicates: + raise ValueError(f"rootHOG files with the same name in different subfolders of {roothog_path}: " + f"{', '.join(duplicates)}") with BatchBuilder(output_folder_big, 1) as big_hogs, \ BatchBuilder(output_folder_rest, sum_list_rhogs_filesize_thresh) as rest_hogs: for hog, fsize in hog_size_tuples: diff --git a/FastOMA/zoo/hog/transform.py b/FastOMA/zoo/hog/transform.py index c7a1c73..ba7db05 100644 --- a/FastOMA/zoo/hog/transform.py +++ b/FastOMA/zoo/hog/transform.py @@ -49,7 +49,7 @@ def _extract_pw(self, node): for child1, child2 in itertools.combinations(nodes_of_children, 2): for gId1, gId2 in itertools.product(child1, child2): self.relations.append((gId1, gId2)) - nodes = set.union(*nodes_of_children) + nodes = set().union(*nodes_of_children) # a group can be empty return nodes else: return set([]) diff --git a/FastOMA/zoo/unionfind.py b/FastOMA/zoo/unionfind.py index 83036e1..6cdad30 100644 --- a/FastOMA/zoo/unionfind.py +++ b/FastOMA/zoo/unionfind.py @@ -87,9 +87,23 @@ def remove(self, obj): """ if obj not in self.parents: return - comp = self.find(obj) - self.weights[comp] -= 1 - self.parents.pop(obj) + root = self.find(obj) + # other elements may still point to obj (it can be the root or an + # internal node of an uncompressed path); they must be reparented, + # otherwise later find() calls hit a dangling parent (KeyError). + children = [x for x, p in self.parents.items() if p == obj and x != obj] + if root == obj: + weight = self.weights.pop(obj) + if children: + new_root = children[0] + for child in children: + self.parents[child] = new_root + self.weights[new_root] = weight - 1 + else: + self.weights[root] -= 1 + for child in children: + self.parents[child] = root + del self.parents[obj] def __iter__(self): """Iterate through all items ever found or unioned by this structure.""" @@ -107,7 +121,9 @@ def union(self, *objects): datastructre at instantiation time, they are added at this point in time. """ - roots = [self[x] for x in objects] + # deduplicate: several objects may already share a root, whose + # weight must then only be added once + roots = {self[x] for x in objects} heaviest = max([(self.weights[r], r) for r in roots], key=lambda x: x[0])[1] for r in roots: if r != heaviest: diff --git a/tests/test_batch_roothogs.py b/tests/test_batch_roothogs.py index 1575819..eb293fa 100644 --- a/tests/test_batch_roothogs.py +++ b/tests/test_batch_roothogs.py @@ -50,9 +50,9 @@ def test_rerun_into_existing_output_fails(tmp_path): run(tmp_path) -@pytest.mark.xfail(strict=True, reason="BUG: files with the same name in different subfolders overwrite " - "each other in a batch (shutil.copy at batch_roothogs.py:45)") -def test_same_basename_in_subfolders_kept(tmp_path): +def test_same_basename_in_subfolders_rejected(tmp_path): + # the file name is the rootHOG id downstream and batches are flat, so two files with the same + # name would silently overwrite each other make_files(tmp_path / "in", {"a/HOG_1.fa": 10, "b/HOG_1.fa": 20}) - _, rest = run(tmp_path) - assert sum(len(v) for v in rest.values()) == 2 + with pytest.raises(ValueError, match="HOG_1.fa"): + run(tmp_path) diff --git a/tests/test_hog_class.py b/tests/test_hog_class.py index 1b1c837..54ee0bf 100644 --- a/tests/test_hog_class.py +++ b/tests/test_hog_class.py @@ -69,13 +69,6 @@ def test_merged_hog_with_empty_msa_is_rejected(sp): HOG([leaf_hog(A1, sp)], sp, "r", msa=msa(), representatives=[Representative(rec(A1))]) -def test_merged_hog_without_representatives_argument_fails(sp): - # characterises current behaviour: the default of `representatives` is a typing object - # (_hog_class.py:280), so callers always have to pass it - with pytest.raises(TypeError): - HOG([leaf_hog(A1, sp)], sp, "r", msa=msa(rec(A1))) - - @pytest.mark.xfail(strict=True, reason="BUG: HOG.__init__ column filter can drop all representative rows " "(_hog_class.py:333 -> _utils_subhog.py:485); the HOG is then lost at the next level") def test_merged_hog_column_filter_keeps_representatives(sp, subhog_conf): @@ -125,8 +118,6 @@ def test_representative_with_elements_collects_subelements_and_species(): assert r.get_species() == {"A", "B"} -@pytest.mark.xfail(strict=True, reason="BUG: Representative(rep, elements) mutates rep's subelement/species sets in place " - "(_hog_class.py:249-259 alias the sets)") def test_representative_with_elements_does_not_mutate_source(): a, b = Representative(rec(A1)), Representative(rec(B1)) Representative(a, [a, b]) @@ -141,8 +132,11 @@ def test_find_representative_by_id(sp): assert h.find_representative(B1).get_id() == B1 -@pytest.mark.xfail(strict=True, raises=AttributeError, - reason="BUG: find_representative calls the non-existent Representative.get_members (_hog_class.py:383)") +def test_merged_hog_needs_representatives(sp): + with pytest.raises(ValueError, match="representatives"): + HOG([leaf_hog(A1, sp), leaf_hog(B1, sp)], sp, "r", msa=msa(rec(A1), rec(B1))) + + def test_find_representative_of_non_representative_member(sp): a1, b1 = leaf_hog(A1, sp), leaf_hog(B1, sp) rep = Representative(rec(A1), [Representative(rec(B1))]) diff --git a/tests/test_pairwise_relations.py b/tests/test_pairwise_relations.py index 6317fba..d473385 100644 --- a/tests/test_pairwise_relations.py +++ b/tests/test_pairwise_relations.py @@ -104,8 +104,6 @@ def test_gene_without_attribute_gives_none(self): genes='') assert _pairs(doc, "ortholog", id_attribute="protId") == [("a1", None)] - @pytest.mark.xfail(strict=True, reason="BUG: an orthologGroup without children raises TypeError in " - "set.union() (zoo/hog/transform.py:52)") def test_empty_group_has_no_pairs(self): doc = _xml('') assert _pairs(doc, "ortholog", id_attribute="protId") == [("a1", "b1")] diff --git a/tests/test_species_tree_prep.py b/tests/test_species_tree_prep.py index b35b7f0..ffa9d3b 100644 --- a/tests/test_species_tree_prep.py +++ b/tests/test_species_tree_prep.py @@ -46,8 +46,6 @@ def test_prepare_species_tree_empty_roothog(): prepare_species_tree([], Tree(TREE, format=1), "r") -@pytest.mark.xfail(strict=True, reason="BUG: for a single-species rootHOG get_common_ancestor({'A'}) returns the root, " - "giving a unary root instead of the species leaf (_utils_subhog.py:283)") def test_prepare_species_tree_single_species(): pruned, _, _, _ = prepare_species_tree(recs("A", "A"), Tree(TREE, format=1), "r") assert pruned.is_leaf() and pruned.name == "A" diff --git a/tests/test_unionfind.py b/tests/test_unionfind.py index f67fe5f..53ef3cb 100644 --- a/tests/test_unionfind.py +++ b/tests/test_unionfind.py @@ -1,13 +1,17 @@ -"""Tests for the UnionFind data structure (zoo/unionfind.py and its copy in _utils_roothog.py).""" +"""Tests for the UnionFind data structure (zoo/unionfind.py, also used by _utils_roothog.py).""" import pytest from FastOMA import _utils_roothog from FastOMA.zoo import unionfind -@pytest.fixture(params=[unionfind.UnionFind, _utils_roothog.UnionFind], ids=["zoo", "utils_roothog"]) -def UF(request): - return request.param +@pytest.fixture +def UF(): + return unionfind.UnionFind + + +def test_roothog_utils_use_the_zoo_unionfind(): + assert _utils_roothog.UnionFind is unionfind.UnionFind def components(uf): @@ -50,11 +54,56 @@ def test_remove_non_root(UF): assert components(uf) == {frozenset(set("abc") - {member})} -@pytest.mark.xfail(strict=True, reason="BUG: removing the representative of a set breaks find() for the " - "other members (KeyError, zoo/unionfind.py:78)") def test_remove_root(UF): uf = UF() uf.union("a", "b", "c") root = uf["a"] uf.remove(root) assert components(uf) == {frozenset(set("abc") - {root})} + + +def test_remove_element_others_point_through(UF): + # a -> b -> c chain (no path compression yet): removing b must not break find(a) + uf = UF() + uf.parents.update({"a": "b", "b": "c", "c": "c"}) + uf.weights["c"] = 3 + uf.remove("b") + assert components(uf) == {frozenset("ac")} + assert uf.weights[uf["a"]] == 2 + + +def test_remove_singleton(UF): + uf = UF() + uf.union("a", "b") + uf["x"] + uf.remove("x") + assert components(uf) == {frozenset("ab")} + + +def test_union_with_shared_roots_counts_weight_once(UF): + uf = UF() + uf.union("a", "x") + uf.union("b", "c") + uf.union("a", "b", "c") + assert components(uf) == {frozenset("abcx")} + assert uf.weights[uf["a"]] == 4 + + +def test_union_same_element_twice(UF): + uf = UF("ab") + uf.union("a", "a", "b") + assert uf.weights[uf["a"]] == 2 + + +def test_remove_root_of_chained_set(UF): + # parents {'a': 'a', 'b': 'a', 'c': 'a', 'd': 'c'}: d still points to c, not to the root + uf = UF() + uf.union("a", "b") + uf.union("c", "d") + uf.union("a", "c") + root = uf["a"] + uf.remove(root) + others = set("abcd") - {root} + assert components(uf) == {frozenset(others)} + assert uf.weights[uf[next(iter(others))]] == 3 + assert root not in uf.weights From 7ec0b80a5ff675bf72d862403c467bf7ecaa6754 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Wed, 4 Nov 2026 16:33:14 +0100 Subject: [PATCH 2/5] fix: input validation in check-input and infer-roothogs Input folders - hidden files (the ._ AppleDouble files of macOS tar/zip archives, .DS_Store) are no longer taken as proteomes or hogmaps; check-input and parse_proteomes share list_proteome_files(), so both agree on what a proteome is - files without an extension no longer crash parse_proteomes - species names derived from hogmap files keep their dots (Homo.sapiens.fa.hogmap -> Homo.sapiens) - an empty proteome is reported as invalid instead of raising ZeroDivisionError Species tree - duplicate internal node names are reported, and trees without internal nodes no longer raise UnboundLocalError - internal node names that are empty or shorter than 3 characters after removing special characters, or that equal a leaf name, are replaced by internal_; generated names never clash with existing ones. The checked tree of testdata is unchanged. - an unreadable species tree exits with status 1 instead of 0 (read_species_tree), so Nextflow no longer treats the task as successful Splice files - a species whose splice file has a single line is no longer ignored by check_splice Fixes #118 Co-Authored-By: Claude Opus 5.5 --- FastOMA/_utils_roothog.py | 33 +++++++++++++---- FastOMA/_utils_subhog.py | 8 ++-- FastOMA/check_input.py | 66 ++++++++++++++++++++------------- tests/test_check_input.py | 29 +++++---------- tests/test_roothog_io.py | 4 -- tests/test_species_tree_prep.py | 7 +++- 6 files changed, 84 insertions(+), 63 deletions(-) diff --git a/FastOMA/_utils_roothog.py b/FastOMA/_utils_roothog.py index e87571d..6beff51 100644 --- a/FastOMA/_utils_roothog.py +++ b/FastOMA/_utils_roothog.py @@ -24,6 +24,26 @@ Gene = collections.namedtuple("Gene", ("numeric_id", "prot_id", "main_isoform"), defaults=(None,) ) +FASTA_EXTENSIONS = ("fa", "fasta") + + +def is_hidden_file(name): + """hidden files, e.g. the ._ AppleDouble files that macOS tar/zip add to archives, or .DS_Store""" + return name.startswith(".") + + +def list_proteome_files(folder): + """(species name, extension) of the proteome FASTA files in folder. + + Hidden files and files with other (or without) extensions are ignored.""" + proteomes = [] + for name in sorted(listdir(folder)): + species_name, dot, ext = name.rpartition(".") + if not is_hidden_file(name) and dot and species_name and ext in FASTA_EXTENSIONS: + proteomes.append((species_name, ext)) + return proteomes + + def parse_proteomes(folder=None, min_sequence_length=0): # list_oma_species """ parsing fasta files of proteins located in /proteome/ @@ -36,15 +56,12 @@ def parse_proteomes(folder=None, min_sequence_length=0): # list_oma_species folder = "./proteome" fasta_format_keep = "" - project_files = listdir(folder) - logger.debug("using '%s' as proteome folder, found %d files", folder, len(project_files)) + proteome_files = list_proteome_files(folder) + logger.debug("using '%s' as proteome folder, found %d proteomes", folder, len(proteome_files)) species_names = [] # query/input species name based on the file name - for file in project_files: - species_name, ext = file.rsplit('.', 1) - logger.debug("%s: %s | %s", file, species_name, ext) - if ext in ("fa", "fasta"): - species_names.append(species_name) - fasta_format_keep = ext # last one is stored either fa or fasta + for species_name, ext in proteome_files: + species_names.append(species_name) + fasta_format_keep = ext # last one is stored either fa or fasta (check_input rejects mixed extensions) species_names.sort() # sort the species names alphabetically -> reproducible results # todo accept all fasta formats in the input prtoeome folder, fasta, fa, fna, .. diff --git a/FastOMA/_utils_subhog.py b/FastOMA/_utils_subhog.py index cd48bcd..564be4e 100644 --- a/FastOMA/_utils_subhog.py +++ b/FastOMA/_utils_subhog.py @@ -70,15 +70,15 @@ def read_species_tree(species_tree_address): elif format_tree == "nwk": try: species_tree = Tree(species_tree_address, format=1) - except: + except Exception: try: species_tree = Tree(species_tree_address) - except: + except Exception: logger.error("Format of species tree is not known or the file doesn't exist "+species_tree_address ) - sys.exit() + sys.exit(1) else: logger.error("For now we accept phyloxml or nwk format for input species tree.or the file doesn't exist "+species_tree_address) - sys.exit() + sys.exit(1) return species_tree diff --git a/FastOMA/check_input.py b/FastOMA/check_input.py index e7d6c2b..3659329 100644 --- a/FastOMA/check_input.py +++ b/FastOMA/check_input.py @@ -21,12 +21,16 @@ def check_proteome_files(proteome_folder): proteome_files = os.listdir(proteome_folder) logger.info("There are %d files in the proteome folder.", len(proteome_files)) + proteomes = _utils_roothog.list_proteome_files(proteome_folder) - not_fa = [i for i in proteome_files if not (i.endswith(".fa") or i.endswith(".fasta"))] + hidden = [i for i in proteome_files if _utils_roothog.is_hidden_file(i)] + if hidden: + logger.info("Ignoring hidden files in the proteome folder: %s", hidden) + not_fa = [i for i in proteome_files if not _utils_roothog.is_hidden_file(i) and not (i.endswith(".fa") or i.endswith(".fasta"))] if not_fa: logger.warning("We expect that only fa/fasta files are in the proteome folder. Better to remove these %s", not_fa) - fa_fasta = [i.split(".")[-1] for i in proteome_files if i.endswith(".fa") or i.endswith(".fasta")] + fa_fasta = [ext for _, ext in proteomes] if len(fa_fasta) <= 2: logger.error("There are not enough proteomes in the folder ") logger.error("Check input failed. FastOMA halted!") @@ -36,7 +40,7 @@ def check_proteome_files(proteome_folder): logger.warning("We expect that all fasta files are with the same format either fa or fasta, we won't include some of them %s.", proteome_files) return False - species_names = [os.path.splitext(i)[0] for i in proteome_files if i.endswith(".fa") or i.endswith(".fasta")] + species_names = [species_name for species_name, _ in proteomes] return species_names @@ -60,6 +64,10 @@ def check_proteome(species_names, prot_recs_lists, proteome_folder): # checking only first record in each species # for ii in range(len(prot_recs_list)): # check that at least 20% of all residue are not DNA/RNA characters (as a check for protein sequences) residue_cnts = collections.Counter(itertools.chain.from_iterable(rec.seq.upper() for rec in prot_recs_list)) + if sum(residue_cnts.values()) == 0: + logger.error("The proteome %s/%s.fa contains no sequences.", proteome_folder, species_name) + logger.error("Check input failed. FastOMA halted!") + return False aa_chars = set(extended_protein_letters) - set(unambiguous_dna_letters) - set(unambiguous_rna_letters) - set('N') total_aa = sum([residue_cnts[aa] for aa in aa_chars]) if total_aa / sum(residue_cnts.values()) < 0.2: @@ -83,21 +91,23 @@ def check_proteome(species_names, prot_recs_lists, proteome_folder): def check_hogmap_files(hogmap_folder): hogmap_files = os.listdir(hogmap_folder) logger.info("There are %d files in the hogmap folder.", len(hogmap_files)) + hogmap_files = [i for i in hogmap_files if not _utils_roothog.is_hidden_file(i)] not_hogmap = [i for i in hogmap_files if not i.endswith(".hogmap")] if not_hogmap: logger.warning("We expect that only hogmap files are in the %s folder. Better to remove these %s", hogmap_folder, not_hogmap) - species_hogmaps = ["".join(i.split(".")[:-2]) for i in hogmap_files if i.endswith(".hogmap")] + # ..hogmap; the species name itself may contain dots + species_hogmaps = [i[:-len(".hogmap")].rpartition(".")[0] for i in hogmap_files if i.endswith(".hogmap")] logger.info("There are "+str(len(species_hogmaps))+" hogmaps.") return species_hogmaps def check_speciestree_internalnode(species_tree): # All the internal node of the input species tree should have a name + internal_node_name = [] for node in species_tree.traverse(strategy="postorder"): if not node.is_leaf(): - internal_node_name = [] if (not node.name) or len(node.name) < 1: logger.warning("one of the internal node in species tree doesn't have a name. we'll update the species tree.") else: @@ -153,28 +163,32 @@ def check_omamer_db(omamerdb_adress=None): def add_internal_node_prune(species_tree, species_names, out_fname): # add name for the internal, if no name is provided, removintg special chars counter_internal = 0 - node_names = set() + # leaf names count as taken: an internal node named like a species would make the names + # ambiguous (pruning by name fails, and per-level pickle files are named after the nodes) + node_names = set(species_tree.get_leaf_names()) + + def numbered(prefix): + nonlocal counter_internal + while True: + candidate = prefix + str(counter_internal) + counter_internal += 1 + if candidate not in node_names: + return candidate + for node in species_tree.traverse(strategy="postorder"): if not node.is_leaf(): node_name = node.name - - if len(node_name) < 3 or node_name in node_names: - if not node.is_leaf(): - node.name = "internal_" + str(counter_internal) - counter_internal += 1 - node_names.add(node.name) - logger.debug("The internal node name was too small or repeated "+node_name+" which is changed to "+node.name) - elif not all(e.isalnum() or e in ("_", "-") for e in node_name): #any(not c.isalnum() for c in node_name): - node_name_new = ''.join(e for e in node_name if e.isalnum() or e in ("_", "-")) # removing special chars - if node_name_new in node_names: - node_name_new += "_" + str(counter_internal) - counter_internal += 1 + cleaned = ''.join(e for e in node_name if e.isalnum() or e in ("_", "-")) # removing special chars + + if len(node_name) < 3 or node_name in node_names or len(cleaned) < 3: + node.name = numbered("internal_") + logger.debug("The internal node name '%s' was too short, repeated or made of special chars, " + "it is changed to '%s'", node_name, node.name) + elif cleaned != node_name: + node.name = numbered(cleaned + "_") if cleaned in node_names else cleaned logger.debug("The internal node name has special chars or repeated '%s' which is changed to '%s'", - node.name, node_name_new) - node.name = node_name_new - node_names.add(node.name) - else: - node_names.add(node_name) + node_name, node.name) + node_names.add(node.name) logger.info("The species tree has %d leaves", len(species_tree)) species_tree.prune(species_names) # , preserve_branch_length=True) @@ -192,7 +206,7 @@ def check_splice(isoform_by_gene_all): total_isoforms = 0 spliced_species_num = 0 for sp, isoforms in isoform_by_gene_all.items(): - if len(isoforms)>1: + if isoforms: spliced_species_num+=1 for isoform in isoforms: @@ -231,11 +245,11 @@ def fastoma_check_input(): try: species_tree = Tree(conf.species_tree, format=1) - except: + except Exception: try: species_tree = Tree(conf.species_tree) # todo add check fro Phyloxml - except: + except Exception: logger.error("We have problem with parsing species tree %s using ete3 Tree. Maybe there are some special chars. FastOMA check input failed! Exiting...", conf.species_tree) sys.exit(1) diff --git a/tests/test_check_input.py b/tests/test_check_input.py index 05ddc0c..6eac191 100644 --- a/tests/test_check_input.py +++ b/tests/test_check_input.py @@ -29,9 +29,6 @@ def test_needs_three_proteomes(self, tmp_path): def test_mixed_extensions_rejected(self, tmp_path): assert C.check_proteome_files(make_files(tmp_path / "p", ["A.fa", "B.fa", "C.fasta"])) is False - @pytest.mark.xfail(strict=True, reason="BUG: hidden files are taken as proteomes, e.g. the ._A.fa AppleDouble " - "files of archives made with macOS tar (check_input.py:39); check-input " - "then reports '._A' as missing from the species tree") def test_hidden_files_ignored(self, tmp_path): folder = make_files(tmp_path / "p", ["A.fa", "B.fa", "C.fa", "._A.fa", ".DS_Store"]) assert sorted(C.check_proteome_files(folder)) == ["A", "B", "C"] @@ -54,8 +51,6 @@ def test_ok(self): def test_rejected(self, recs): assert self.check(recs) is False - @pytest.mark.xfail(strict=True, reason="BUG: empty proteome raises ZeroDivisionError (check_input.py:65) " - "instead of being reported as invalid") def test_empty_proteome_rejected(self): assert self.check([]) is False @@ -65,13 +60,9 @@ def test_species_names(self, tmp_path): folder = make_files(tmp_path / "h", ["A.fa.hogmap", "B.fasta.hogmap", "notes.txt"]) assert sorted(C.check_hogmap_files(folder)) == ["A", "B"] - @pytest.mark.xfail(strict=True, reason="BUG: hidden ._A.fa.hogmap files (macOS tar) count as hogmaps " - "(check_input.py:91), so the hogmap set never equals the proteomes") def test_hidden_files_ignored(self, tmp_path): assert C.check_hogmap_files(make_files(tmp_path / "h", ["A.fa.hogmap", "._A.fa.hogmap"])) == ["A"] - @pytest.mark.xfail(strict=True, reason="BUG: dots are dropped from species names (check_input.py:91), " - "so they never match the proteome names") def test_dotted_species_name(self, tmp_path): assert C.check_hogmap_files(make_files(tmp_path / "h", ["Homo.sapiens.fa.hogmap"])) == ["Homo.sapiens"] @@ -86,15 +77,11 @@ def test_leaves_missing_species(self): def test_leaves_duplicated(self): assert C.check_speciestree_leaves(Tree("((A,B),A);"), ["A", "B"]) is False - @pytest.mark.xfail(strict=True, reason="BUG: name list is reset inside the loop (check_input.py:100), " - "duplicate internal names are never reported") def test_internalnode_duplicates_reported(self, caplog): with caplog.at_level(logging.WARNING, logger="FastOMA"): C.check_speciestree_internalnode(Tree("((A,B)x,(C,D)x)r;", format=1)) assert any("unique" in r.message for r in caplog.records) - @pytest.mark.xfail(strict=True, reason="BUG: UnboundLocalError for trees without internal nodes " - "(check_input.py:106)") def test_internalnode_single_leaf(self): C.check_speciestree_internalnode(Tree("A;")) @@ -116,14 +103,20 @@ def prune(tmp_path, nwk, species): def test_rename_and_prune(self, tmp_path, nwk, species, expected): assert self.prune(tmp_path, nwk, list(species)) == expected - @pytest.mark.xfail(strict=True, reason="BUG: a name made only of special chars becomes empty " - "(check_input.py:168)") + def test_generated_names_do_not_clash_with_existing(self, tmp_path): + # internal_0 is already used further up the tree, so the unnamed node gets internal_1 + res = self.prune(tmp_path, "(((A,B),C)internal_0,D)root;", list("ABCD")) + names = [n.name for n in Tree(res, format=1).traverse() if not n.is_leaf()] + assert sorted(names) == ["internal_0", "internal_1", "root"] + + def test_name_too_short_after_removing_special_chars(self, tmp_path): + res = self.prune(tmp_path, "((A,B)'a.b',(C,D)cde)root;", list("ABCD")) + assert res == "((A:1,B:1)internal_0:1,(C:1,D:1)cde:1)root:0;" + def test_special_char_only_name(self, tmp_path): tree = Tree(self.prune(tmp_path, "((A,B)'!!!',(C,D)cde)root;", list("ABCD")), format=1) assert all(n.name for n in tree.traverse()) - @pytest.mark.xfail(strict=True, reason="BUG: internal node named like a leaf makes ete3 prune raise " - "'Ambiguous node name' (check_input.py:180)") def test_internal_node_named_like_leaf(self, tmp_path): res = self.prune(tmp_path, "((AQUAE,CHLTR)MYCGE,MYCGE)root;", ["AQUAE", "CHLTR", "MYCGE"]) assert {n.name for n in Tree(res, format=1).get_leaves()} == {"AQUAE", "CHLTR", "MYCGE"} @@ -146,7 +139,5 @@ def test_no_isoforms_exits(self): with pytest.raises(SystemExit): C.check_splice({"A": [["a1"], ["b1"]]}) - @pytest.mark.xfail(strict=True, reason="BUG: species with a single splice line are ignored " - "(len(isoforms)>1 at check_input.py:195)") def test_single_gene_with_isoforms(self): assert C.check_splice({"A": [["a1", "a2"]]}) is None diff --git a/tests/test_roothog_io.py b/tests/test_roothog_io.py index b93ee03..d03cf5a 100644 --- a/tests/test_roothog_io.py +++ b/tests/test_roothog_io.py @@ -39,15 +39,11 @@ def test_fasta_extension(self, tmp_path): write_fasta(tmp_path / "A.fasta", [("a1", SEQ)]) assert U.parse_proteomes(str(tmp_path))[2] == "fasta" - @pytest.mark.xfail(strict=True, reason="BUG: files without an extension crash parse_proteomes " - "(file.rsplit at _utils_roothog.py:166); check_input only warns") def test_non_fasta_files_ignored(self, tmp_path): write_fasta(tmp_path / "A.fa", [("a1", SEQ)]) (tmp_path / "README").write_text("x") assert U.parse_proteomes(str(tmp_path))[0] == ["A"] - @pytest.mark.xfail(strict=True, reason="BUG: hidden files such as the ._A.fa AppleDouble files of macOS tar " - "archives are parsed as proteomes (_utils_roothog.py:165-170)") def test_hidden_files_ignored(self, tmp_path): write_fasta(tmp_path / "A.fa", [("a1", SEQ)]) (tmp_path / "._A.fa").write_bytes(b"\x00\x05\x16\x07\x00\x02\x00\x00Mac OS X ") diff --git a/tests/test_species_tree_prep.py b/tests/test_species_tree_prep.py index ffa9d3b..2868517 100644 --- a/tests/test_species_tree_prep.py +++ b/tests/test_species_tree_prep.py @@ -59,15 +59,18 @@ def test_read_species_tree_newick_with_internal_names(tmp_path): def test_read_species_tree_missing_file(tmp_path): - with pytest.raises(SystemExit): + # a non-zero exit status, otherwise Nextflow takes the failed task for a success + with pytest.raises(SystemExit) as e: read_species_tree(str(tmp_path / "missing.nwk")) + assert e.value.code == 1 def test_read_species_tree_unknown_extension(tmp_path): p = tmp_path / "sp.tree" p.write_text(TREE + "\n") - with pytest.raises(SystemExit): + with pytest.raises(SystemExit) as e: read_species_tree(str(p)) + assert e.value.code == 1 def test_list_rhog_fastas(tmp_path): From 54463b1fd242bd96c4ee56d7b7bb2e8a5ddd3ae3 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Wed, 4 Nov 2026 16:41:21 +0100 Subject: [PATCH 3/5] fix: report works when CompletenessScore is disabled The four report cells that plot CompletenessScore now check that the column exists and print a note otherwise, so --report together with --disable_completeness_score no longer fails. The options nf-test now disables all three scores (it kept CompletenessScore before because the report failed without it) and checks that report.html is still produced. Fixes #123 Co-Authored-By: Claude Opus 5.5 --- FastOMA/fastoma_notebook_stat.ipynb | 68 +++++++++++++++++------------ nf-tests/options.nf.test | 6 ++- tests/test_report_contract.py | 3 -- 3 files changed, 44 insertions(+), 33 deletions(-) diff --git a/FastOMA/fastoma_notebook_stat.ipynb b/FastOMA/fastoma_notebook_stat.ipynb index e2e83b6..4d7caab 100644 --- a/FastOMA/fastoma_notebook_stat.ipynb +++ b/FastOMA/fastoma_notebook_stat.ipynb @@ -398,8 +398,11 @@ "id": "c90ff7ef-c88f-4c74-bfbb-0cf2fa0d5901", "metadata": {}, "source": [ - "sns.displot(roothog_df, x=\"CompletenessScore\", kind=\"ecdf\")\n", - "plt.title(\"Cumulative distribution of CompletenessScore for all RootHOGs\", fontsize=16);\n" + "if \"CompletenessScore\" in roothog_df.columns:\n", + " sns.displot(roothog_df, x=\"CompletenessScore\", kind=\"ecdf\")\n", + " plt.title(\"Cumulative distribution of CompletenessScore for all RootHOGs\", fontsize=16);\n", + "else:\n", + " print(\"CompletenessScore was not computed for this run (--disable_completeness_score).\")" ], "outputs": [], "execution_count": null @@ -409,18 +412,21 @@ "id": "073a3976-c942-42dc-815b-7fc1655d8837", "metadata": {}, "source": [ - "g = sns.jointplot(data=roothog_df, \n", - " kind=\"scatter\",\n", - " x='nr_members', \n", - " y='CompletenessScore', \n", - " marginal_kws=dict(bins=20, element=\"step\"), \n", - " marginal_ticks=True,\n", - " height=11)\n", - "g.fig.suptitle(\"Distribution of HOG sizes and CompletenessScores for all RootHOGs\", fontsize=16)\n", - "g.ax_joint.set_xlabel(\"HOG size (number of member genes)\", fontsize=14) \n", - "g.ax_joint.set_ylabel(\"Completeness Score of HOG\", fontsize=14)\n", - "g.fig.tight_layout()\n", - "g.fig.subplots_adjust(top=0.95)" + "if \"CompletenessScore\" in roothog_df.columns:\n", + " g = sns.jointplot(data=roothog_df, \n", + " kind=\"scatter\",\n", + " x='nr_members', \n", + " y='CompletenessScore', \n", + " marginal_kws=dict(bins=20, element=\"step\"), \n", + " marginal_ticks=True,\n", + " height=11)\n", + " g.fig.suptitle(\"Distribution of HOG sizes and CompletenessScores for all RootHOGs\", fontsize=16)\n", + " g.ax_joint.set_xlabel(\"HOG size (number of member genes)\", fontsize=14) \n", + " g.ax_joint.set_ylabel(\"Completeness Score of HOG\", fontsize=14)\n", + " g.fig.tight_layout()\n", + " g.fig.subplots_adjust(top=0.95)\n", + "else:\n", + " print(\"CompletenessScore was not computed for this run (--disable_completeness_score).\")" ], "outputs": [], "execution_count": null @@ -481,8 +487,11 @@ "id": "ee844afe-fa0a-4b3b-8fe3-4e9fec48bf36", "metadata": {}, "source": [ - "sns.displot(level_df, x=\"CompletenessScore\", kind=\"ecdf\")\n", - "plt.title(f\"Cumulative distribution on CompletenessScore for HOGs at {level}\");\n" + "if \"CompletenessScore\" in level_df.columns:\n", + " sns.displot(level_df, x=\"CompletenessScore\", kind=\"ecdf\")\n", + " plt.title(f\"Cumulative distribution on CompletenessScore for HOGs at {level}\");\n", + "else:\n", + " print(\"CompletenessScore was not computed for this run (--disable_completeness_score).\")" ], "outputs": [], "execution_count": null @@ -492,18 +501,21 @@ "id": "3f572007-d324-4310-bc55-b210e19034f8", "metadata": {}, "source": [ - "g = sns.jointplot(data=level_df, \n", - " kind=\"scatter\",\n", - " x='nr_members', \n", - " y='CompletenessScore', \n", - " marginal_kws=dict(bins=20, element=\"step\"), \n", - " marginal_ticks=True,\n", - " height=11)\n", - "g.fig.suptitle(f\"Distribution of HOG sizes and CompletenessScores for HOGs at {level}\", fontsize=16)\n", - "g.ax_joint.set_xlabel(\"HOG size (number of member genes)\", fontsize=14) \n", - "g.ax_joint.set_ylabel(\"Completeness Score of HOG\", fontsize=14)\n", - "g.fig.tight_layout()\n", - "g.fig.subplots_adjust(top=0.95) " + "if \"CompletenessScore\" in level_df.columns:\n", + " g = sns.jointplot(data=level_df, \n", + " kind=\"scatter\",\n", + " x='nr_members', \n", + " y='CompletenessScore', \n", + " marginal_kws=dict(bins=20, element=\"step\"), \n", + " marginal_ticks=True,\n", + " height=11)\n", + " g.fig.suptitle(f\"Distribution of HOG sizes and CompletenessScores for HOGs at {level}\", fontsize=16)\n", + " g.ax_joint.set_xlabel(\"HOG size (number of member genes)\", fontsize=14) \n", + " g.ax_joint.set_ylabel(\"Completeness Score of HOG\", fontsize=14)\n", + " g.fig.tight_layout()\n", + " g.fig.subplots_adjust(top=0.95) \n", + "else:\n", + " print(\"CompletenessScore was not computed for this run (--disable_completeness_score).\")" ], "outputs": [], "execution_count": null diff --git a/nf-tests/options.nf.test b/nf-tests/options.nf.test index f7b8f43..398e290 100644 --- a/nf-tests/options.nf.test +++ b/nf-tests/options.nf.test @@ -16,7 +16,8 @@ nextflow_pipeline { fasta_header_id_transformer = "UniProt" write_msas = true write_genetrees = true - // the report needs the CompletenessScore, so only the other two are disabled here + // all scores disabled; the report (on in the test profile) must still be produced + disable_completeness_score = true disable_implied_losses_score = true disable_tcs_score = true } @@ -33,7 +34,8 @@ nextflow_pipeline { // UniProt transformer: "sp|P13927|EFTU_MYCGE" -> "P13927" in all outputs { assert roothogs.flatten().every { it ==~ /[A-Z0-9]{6,10}/ } }, { assert FastomaOutput.canonicalHogs(orthoxml).every { !it.contains('|') } }, - { assert orthoxml_text.contains('CompletenessScore') }, + { assert !orthoxml_text.contains('CompletenessScore') }, + { assert out.resolve('report.html').toFile().exists() }, { assert !orthoxml_text.contains('ImpliedLosses') }, { assert !orthoxml_text.contains('TCSScore') }, { assert file("${params.msa_folder}").list().size() > 0 }, diff --git a/tests/test_report_contract.py b/tests/test_report_contract.py index bb57adc..021fd17 100644 --- a/tests/test_report_contract.py +++ b/tests/test_report_contract.py @@ -88,9 +88,6 @@ def test_reads_these_outputs(self): assert '"species_tree_checked.nwk"' in source assert '"phylostratigraphy.html"' in source - @pytest.mark.xfail(strict=True, reason="BUG: the report plots CompletenessScore unconditionally, so " - "--disable_completeness_score with --report fails " - "(fastoma_notebook_stat.ipynb, roothog_df/level_df plots)") def test_completeness_plots_are_guarded(self): for cell in _code_cells(): source = "".join(cell["source"]) From 85af6ed3f6964ee6e4f915d39e45bc0afc3ad016 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Wed, 4 Nov 2026 20:49:14 +0100 Subject: [PATCH 4/5] fix: skip hidden files in the Nextflow proteomes channel check-input and parse_proteomes ignore hidden files, but FastOMA.nf selected proteomes by their .fa/.fasta suffix only, so the ._.fa AppleDouble files of macOS archives still reached omamer. New nf-test runs the pipeline on the original macOS-made test archive (with ._ files for every entry). Refs #118 (review comment on #125) Co-Authored-By: Claude Opus 5.5 --- FastOMA.nf | 3 ++- nf-tests/helpers/input_archive_macos.tar.gz | Bin 0 -> 9973 bytes nf-tests/lib/FastomaOutput.groovy | 7 +++++ nf-tests/options.nf.test | 27 ++++++++++++++++++++ 4 files changed, 36 insertions(+), 1 deletion(-) create mode 100644 nf-tests/helpers/input_archive_macos.tar.gz diff --git a/FastOMA.nf b/FastOMA.nf index b14c23c..52ea4e6 100644 --- a/FastOMA.nf +++ b/FastOMA.nf @@ -495,8 +495,9 @@ workflow { proteomes = proteome_folder.flatMap { dir -> // Use file() to list files, don't create new channels def proteomeDir = file(dir.toString()) + // skip hidden files, e.g. the ._ AppleDouble files of archives made on macOS return proteomeDir.listFiles().findAll { - it.name.endsWith('.fa') || it.name.endsWith('.fasta') + !it.name.startsWith('.') && (it.name.endsWith('.fa') || it.name.endsWith('.fasta')) } } diff --git a/nf-tests/helpers/input_archive_macos.tar.gz b/nf-tests/helpers/input_archive_macos.tar.gz new file mode 100644 index 0000000000000000000000000000000000000000..5c5a7f0f09597c49350830315400b6870671b395 GIT binary patch literal 9973 zcmVji-?bazaCG`?pM2ol*5Or8y}!t@8-Xn#kq97j zG64+Q>5GCRwLFRyNexMjXCw4M?6W_Ezl_iR3cmQ{-;-6HBQ4I(|R?ynHfWFE(rbtrpuS zv-xtf+Uq5~N&3o_Twy(bWXP$H1%9FeA|G$iX@o)a+ zCx89PC+*WOp8P63`HxRdLtgvjimwGQT+Yw^QX_hxc%g>{`!;pfB)P6^WXl*|NLJ+!X_Wb`$qrt56$h< zFK=J|^26c459-pU>A&d9TlRmoMgP|ya7}mPXRiNu{J;KO4Zl|A`m3iu z`W(<>vs(Qz{@;4GUOoDMA7woHe?RMSr~mhk;E#?2KdArj_5(f|9&k30SUVAlVK^nbosZ{~05|7x>ZKidC~G5-G5x8MJ2v)#;h-v`sHqZ3%Q zX7F8!U43v)HAOepryP`yA($~+<&ufUM>B+=bI^s~qIJ3Wc2HG`Z9RJBwe9lgLld0N zz2zx&8cj)21*Hzgst~xM(qNrYLoTYTZ+xYS8FVPl^~H7G>Ch*oi;l(n!5LG64Z1f* zxuBBr#yZ|@T;PS$Rnhf*@2&0!-v_UPt92Sj-H*|#qUsWP!J6WGtG#POo3izWHDyz5 z#Tt6wW?yR`lCrwk(J7Uj*ClEfoz>apM)|JK-t?un_0TKdw`Ekea$~QI&UrK@M6H7H zqjH^AENe8qNvT)CwziJN)@&g(t#M^AF+jO$w8cA@v+BFzM%k*bvo^&USkSk+?scbZ zaAqvV^mPd(b9>LyS$6N-5Do*}unwlRd4xb+E+d>WZ6A_x4$?R5*e92wtyLQJXliH2 zD)vbS=R+<5YWu#h^i;djz-l%4TE$qD@gbV*`Z$Ezl*SE3d24-DoI(s7q>id9S^2|g zAzO~DA35m;*yDX{{ZLz9hR_d2)ghJOi>Xd>H(M^gkM;a4cEKApvd91r*wPr~n+eJH zP-8&O{K&JDH{8e-*O@Sy$ZKk-rWX;35PVA7YZe`pL1elZvN0`hHV)`mr*}mL`@VmL5?` z$wt@Yloc{v0xQ=hYX|vznckUPN0ns4VMJzgHd!H=xp95Au9V?0XcN;=51k(nVviM7 zK6pK{RkQ$cXP42Np;R_oZZ+P;Ja%3QWKLk8T6?6IY5Yvwhs<@ODe^zBUBurONcx% z2HOE<^$pxHRXxT!r!yN@i~aYe&c8nMkxAG|1g@F_C9<=3P@b7O9|}lpkp(D?_SPQs za@s586m1efcO7zB01~u2Hlvt=15mE?=ztoBB~U}Ft!%L*G{^RBKjbPhAykb6)Al1W zUALV-=(cE$3QXb`z&3tN%sB$b241?zhLIqJW(PBLski>77!a5nJJo74q9Q$VRJ8Bw z(DpE+)PS0e(6(e7D=ZQ67pEtJuoe(7m;y90!@9PP9+5^0xS@008MxhgEf&mH1N0br zn~GPBk7(&o#v*6~fnt;Ux-Q)q`=K7|n6XM(X7y+|l2^5*m4{9{Gr_68pcsY))XGJ9 zKrZM!t3){T8UEI*?f2=|2Ytse2-$gnsel8}A1Sn{#^?+O(B9^u6rO=K1%O-TK~+i( z*aLhO|G|&$kT83tSG^z^iYg+nkVQCjaCv%%*KOPK{|Mf8W8W!XA-sJDYYTS8f}+Jb zw$r#8v*NL_nC)Z{n-Pp0hmidU7!F_@H&&p24_t)YL|uf+Y@tV@*-h)LaTawqg03*o z_=$U^2IpB;u2j!%VvA5AMy^E^EyAXJE5L%1s*y)_SX!P08%4KyL7?2kDkrVZwuGo< z^L_4JIukQS*r2)_kMIeJSj6bMTnAql(_oYVq~H#Hj|znm!umE0ASx_H2#4;_jAB() z4a}kVJy-D)4p?0TtHJa|=wnUNu)6_tiaG?-LXX(Ha3~g}HG+?}&9$sKqtihJCNQSS z*kkP!&Oy0RK!eS_%DvIz0F?$Ph#e+*2!Kfk80Cl+u!v>{$fMVF>l>7{!<9g>RT8Ht z)fMC^Bi~QrwOj1I=c)g2Ctlbtk1Z}kakei7u?OiOO#rg+*H9D}L3gbI{hM0ctr{^? zLj#iyCXQ>WwSDH+5-sGIV7kY+^{A0}XSi4$q#Syu&A=RDk;m$$NA@{Yv(Ez>&AfQ=meL*Gt{Zm;Ra%RF}-B4(%NN>X5A!Ls|K~oEHfi!N{kjBpwm4jG6QOulwRC3IE{a)AfktZ1Fp*Cw>YtS zt#~%t0J!+B5myTvgMT#ibyMMW4aglj=ttxcN9mqGAOBxZWUSo1`cx&xI6r7F{u_NMFV+-)(yy!pvACJVrRR9eX9qZ1- zAr2~gK+1y%peN2yfRPz*fMo+TC>OxWt7X&(h|RSYpA1k%cyW?F&L`e^gD9#Fxr?sV zH^EiL)LuZ59Y?Q$Y#1QxNM=x6z|#d~GFlY+8e1oBb2sz^EQ1FTGb{<>0CtBU?n->r+oDvKz zek?0wMa6W0vL{|*hqy!{iXEy)FoVzB`RIutk#!Ixg?dCCP;cBj;tt}Z(S!~tz|cSo z2|R$e2H!JB9)!El8)7u)Tg!-C;(9XfKB2l0}AGmgvo46Bn%*; z0=7VGTvQ@b@Y*GAKCtFPQ2hbrM445Mq~npdB@Mp0uWL+BMA>mROP{ki zl|UkL)+p)h9MWfQOrlon-@K#2|QLUqP>7Jm9S3x)6OK zK7Er^Nf{5_4wRwVSQ*p3D`@jMmb0%8tT#gRnS&>jgGaX^(5 za~CjS5+!zJvK@F8*oh)uh6Irk10!=11&SuBjL>x`2)-5&kKfk{^vcXI$e2`lVj1)k zIP$g~YV7+MQa{#&i)ce@2nBGmMkmU{Bm_?4bmPZJlv(io5PWE(jF?h9GnMftyT%=( zy#p&#JwiA{fEJek&qgej9iTXo4Jw5VOVCHRxQ&7)$Q{y>GX}{@#>U_rDFPh6#ewiH zWlw&H}{rPkEsM$5{G5(oAR6!b5Z|BVdp6D2WFP_N_^hZiJ06 zzvFHoI6-iv5*y?tLq{w^MA#L)b6nB^)isj6>4?1}-dG3HhD7LsjRqaqVT;XzZvn&3 zW&jP|fHPPlDFOY;1N`Vv{wNV@g@{Pv0-FrQ2F$h_d2@s6;0f{?mKAMAQAWi}crxIj zAWuRRD#^O=c$CRcAz1D4?n?icR_M5>|F`48{>9Y)H;?uIk24hzm$ALek z{y%@L|9_P6p7Wpc)nYSWJ=Xt!^5ahbKbZCZLH%E^R@-;#|C`05{r?z4f|P!>UC(Ay z9d?p}s@1jjL#XRcLIqsmY$PkBL+AT0q*|ghDbEvwHGb#=u}3Y5O^@4Lb&7bOJ3`l` zBz4+mN*swVj0c`VpCozjd>sfdi#A;+XAPJIpS@DbT5Zz7)gg%|-;TZ;4>=J{6Uc_n zI7e6^DUZ6OT9wXbVg-rAYg65qVOg*u5HC`OEypq}8<0#|cAeC`Bk|05F2gpZi&Zx1 zScVP9GORib#`-kWsf$`7hgRyR;fCN(3bP^QK}lXQ4aG9EiRx5_9iP?9Z`HH>N-fSKnm=M zFg1?)T(~2Mww3HMDJk4ZCH#jAcUh+*g}Z1PR~-v?l2mL6O;zQmi6+7SCuE6|rM=tk^k8D%#}hbWyR(ayld&i%clWX1PC>R8DM?Y-LA4Z+&U*sw%Ed zCLtqzETcXUf(|T+5VI!S%tIis)~>|rz*20ViPcr?Bz@G4UDG?Mg+f>mJ=*>s;pB|Oi7{BxeQWFkrIM$reeyV>UtC{8;mMaOo^cz3@DiWne$Ss zoeC-^&L^ci8!iee2PvqC-1KG2rtxrOu1`)1B%MH4kphX2L)F4-IL8&oOs&nS^9aDu zDMR3F1q)r)rMzoa_PrArOt#5P(0nO-=lxLBZfrm=2(YLA{s(7s8fP;=DDiAC?2cp*j zEg3w*FZFLPr4y+Sh3tDt4CUDMu1e^I9MGB0!LUddPPLP^S7kCCsVzNw?7NuUsdnO~ z+6l`psJ$*;BX?mco@DgEJjIz?kjM}8xPQjY&$2%!My79bespMGacY|6OC!Hk7I}3k&@}4E<|oQCHv1N zfrmz>ki`BeeBOWc&xM`i3Q*GXCYwj6b%hhqPcuXWJC8gQC9;5zAp~mMd zKk~(Y23sV@{-Aip)1?N~`T%T<$4q}pBS0NWo>BrpwdCIw5GFms+k|Nv{(i`VaY|7HqE?$x9ww#aBqjuRdr31&hpKf@1k^en30QX-?{&)MB z|NS`QG5`Cs9uMSye{>x9Q}Vxy$NcX{8Slyeu6E1$eEpdJ{mG9z{r_Or{|EcOX6yCl z-Td$N(f)soF&$Fbt!A6;vHxolgQ~@yCdTh1iISuR%pLm%+_42LWqve-AMLRPEVVkE zTELvIRMVH-=FyBq#aRyHz1;zp?mEC+At)Yf8PvdwU8%;PcodJ54zO_S0OJngLU-%{ ztBDhJPK4U(kd)7Mu&$R>4zao;9+nQUz8_3B=MJz&I>6*WJMj>aY_bx~5j7F?w<>#^ zlB7Wh+N`Pr6Xnnnsf9>cLA7SRr)sjl{wPWLC0BzAu%e)7;oQO&D71 z9K4o3nJ^F%B^88Nrb|r9A6hM^#16q{-3C20kRm|s&`Z+1u&R!LmOzpX6C&4o zjKNlVgc8ch=`0sqCuf1Rs?&iWUXtDfMm8pQFkViqI><`s?1h*QEIo+sUR zu1<2iP{QcsO{1dpwfRmpCWQhMy$ge`#4b*Gry@tL z2BQ+MCWUyW3PF}Gytf*_tTcdiK}(2iiJx0(FAOC&Fcv1-#+a7IjeQPslGEib+gw|t4(A52>2#AcfO%;EtHxn!02`D`QadW4MQjIKr&8;| zW$8$1ot7+T-$`cd)Bsjc{!p=5rMeu3*oQGnX8~b<0}T&-+e$lxwQZwvXTwlW4PddC zv{aI6k93{Nk^Z_37=j#My>9?hsgu)sI)osH__DHf?xWEzRm`pr@H^FGmj~_h zC?}Ffr6IAC&4IsqS{uQLv z4PJ)K*976MOOOmQ79EOdr~a=X&6T6Khu$3z_;OM|a!GKq`uo=}Mb~zV{ps78QDJI> z(~>aD$DS`GJzqHwhoXS$4oo`te5Ik5L(^i_07+H(0TDk$m0Icf>ZRuk`Kk?&GxdDA zbI+IaSOqz6jBWD?K@YAvcYJyNJ4llk`x#vmdg;U|7ad>aW(awR+%hB{2>C!u?265) zmTV-t2(qOCH7w_jFX^2GQ1wlYLz5CV(kV=-6$@ASi;k~^qFFIIevmd!>63u3MgV~) zp$^h8+{ZIkEtiX9lh+xnvWC}9gjIdMk)(H%OwYsWQrlDSm*B%7eP1jW`|L28{otJB zJDqfYL4S1^khY`~v;~xgoExKbIyZl*zQv?FiwA=`OY>JXb$>Nx=x|e%5#&;?NwyWY zw1aGJDK>2jm$L-b%l3eZH#G?D!to4>r$ zovLu3u=l+#$L=pl)Td;m+U49ZIPLmVO?kUp%#S6jGbR|xC3pUTqSerQ>E}+y=2Ov1 z)d(Eh00Dcc0YIPJ^#*cEHSg2}lrYg^i8JcXCQv150s`!d8GIP0juWe@!G+f21_=Ct z>qd&I7iFw8M0qAM8>Nh;rCAE6?pVgEj&&?!HMlO14#yTCCoMou&Pvylh3-x*KpF=O zsdn8Beo!SXK!cI(Ld>o`wg3S#Qsa{%7N#dj3(#;+vP6=NCS|Rqu@ZBjzuh0q`jw>S}pF8EpHaFu!~k2fKD|mH#GpAYFKUV0mWkt zD-RyUo{CuTd*H!f*{K1jof?4HcC4j-dTIb_q=r?5CbEG_8h~((wKNdmpnC^KDro@n zMjJ%rSi=epD2I<0r6UTGcL=;~D_6@&1|{vG@;3u|L2cG9uDSC9oQ!G-~0|i`D}u4Gy^!DO_{W z85FcMU9b&2YAKvx6p3KEs#pxX2Nn1?NGiX7~Sc;MXI6VwEZR>Vr? z3j2wn0E~u&2q9rvT85|YAYi1&QGkst+KuYCRu`}|LTrVXYpYTEouo&o>xwkF$2WU~ zaBN#d<;J8s#OBohjSFDr$8zC`aM6{EIw8>NaH_i`2MmsdLg_of7aDQuM-3M-JF3T` z)RtBDn4EXT(J5+*F2$g%`ntj z`pDkw3?gJe546+?rA@lcQYKWn0xy#^CG{zKeG@UGfbh^w=YEcbLX0>@``8+kZL4Z& z4dVHUFicv5P)=MDF0(HJ0``r)D{s8j|I z%s0}9hDUYZ6XY*@g4oF@$3&t8S-8wS@%DXBko3|F4)ZKlc4$y*T+h3nAjL*HaYz-G zU(S_Cd9E079&Byo%mqSi+pg}Wo}e_=N?<+@Q%?}qUD_T5dXtq_m4JC0hEfq5PCY?R zdV<_VPY|JqmA`Xx*At{;^b(3`Ij12#L7@_NFi2h2bOhI2l)pcK{Xu6ZKUfxI%@kRv zc#yvUAfZ38_0$vO-|7h}=bj)Z4Z+^yjca&Hq~K0HK|y+gh<{ua58$gbeFQ9e?!ZGF zbAl@gvB30WSCEqaRU^(x%~k<-2x%-cEBV_lqXgK}5fqSW_{-H`(y=2*LW?|g1c4AZ zxF|HxgwiD_@!_ri59|MIwRq?J&+^g! zf0Xeze>1hE|LqjX{JXE8zrKAr|J&&ZA>*pVA zj@cjo`9IsozyJD?#-soLvmW>U|A)hYAJqTZ+x?&G^=$s=|9_0J*gxR^@3xE8V*U8> zKR<(Ur}ZE31Au-=|7Y)h{?BH%eEj^Mk1_7~fBykM*T)xsv^oCx{O@M9eDwc5&Up0y ze%9k&|KI%ppdSkd{^&*^`v0M<{}1Z_a<*E( zZT}aG<)i)o7~>Z&o;`hb`|{J@-@bhH^^519Jy~4OuNT!Xe)rw&%YXXx*_Z3jo_zcA z7WMu1yVtj0e)`qde|Y`f%iB+%J$?C)w=b(-eEIbC?WZrlllRRQn@{JnPZzU1Uw*dO zf3{d%&vE-#bME@hr{6w({oBu;yuN+)`uex)>6NdafA!+iS1-PM`Ni#LPkz;^Cs6qL zmruWW@%;A5FQ+?SeE0g7_0w0cx%heWEBosD#mm>fees89-`v0d<%<`uKmX?T_qX4C z_T-qazQ4xWp8n?PtJ}|i`}7raf4Tk5Zuy(r-Ttffd^i8{%gubf`{D~QZ?oPmmR~HY z-@Us1-Pd1UA^ykzr^jA>_4L`--~7|(-=2T?;>CB*U!Q+@{^I4c;|s68`_1XKS6{q% zd3*8u@wzq39n1g=%@?~F+gq%sOZ+~4aZB~U2k=}>vOi9; zxY{ij*X#LWJ(D(pbzTkbcz)CF(T(O<$%91G-``yJp-uw`Yizyc8DYmQG^=h%(&Q{AAG@0`R z2zHud4F{I;xeu#_#N5**7IL-N%r5rv=I27(PcomMiN#|#*Q?EJznm|Za~`(3V)^G; z=4$|Y4~>w(-Eu!&T5aX~_F|SdKY`$?`=xsC*Jimq|6D!vSK=?|Fux=9YJ0t0BSmm! zv7awiSBu@|j@G*^eB2`XEEM+0rNtgv?Dm%{y!>P87i3u63IBX~JzpW=yUlXG*shT5 z<(=@)=kxvM?*HSZ<^3(yLw`T~Vv^aN^e;Ax>y7wc%f)_+6mG88t9!{@Ef%}Q7D1lx zWfyX3y^`-2lKIvjnZ6*`_MW4=)%9kFuENtbT)kQy1#g;UyPFB(tTt;AN4bK+U!joG z9Di;j^*BD*xXk~YTY&!(@;~dx{LjZ3kNKaU^>`@%^JC+{pOF7q%pUVUA7wn4|H0?q zE^urf^FKfNaj*X$%KHDH{u7(5-pT*09`ipRWjyA8{xb7FZ~sY__vC*z*Fg1>plvhX zZa2%T^?v(s{%18^+7Jp~=6~M$>l+W|e-`WO)qcGt$|aQE5aO;^?@#}%u9lnG8|k07 z{=CF{(mw>iMA$5IK@7TBU#;f*hx0$%>5}ArF7iL`{P}|i^FMP*%`6t1{hIi9N0hzX z?H|nlOc$2>^L4#R?EfpVz%Ec_r@Q3y^&^UD|BsSDVH5Qnudcj%HsVw9LeA8E^|0(%-$h4bgc|I7IQ{1Fg;3Gx5q`QMK;9^?O?^>`@$|FLo4Pl*5L+sF9-ql|~+|NRnc zxp|ELfBNHI|38%V|3UpHN__k8e-i&ce*V)(86Q0U-_9Q6|35zdfBPf3K0svt&Sy^j z$U*Z%pO^HWqXzS9o{h?H_q*+qKzFshJN@zQQG@w(%;4gv!8@PL^I$wn)Oby}xhD2r zuLNzb*8ADR@$7QCwBBznMyUDr z$B;Z&$X_gG*Q+H=-)v^f#gf&|m+vX$%a!$f|5(U>9FOC1JTU$*4!98G0Hy!{E9~bw literal 0 HcmV?d00001 diff --git a/nf-tests/lib/FastomaOutput.groovy b/nf-tests/lib/FastomaOutput.groovy index 8eccc4a..f09d43a 100644 --- a/nf-tests/lib/FastomaOutput.groovy +++ b/nf-tests/lib/FastomaOutput.groovy @@ -110,6 +110,13 @@ class FastomaOutput { return counts } + // Sorted names of the species listed in the orthoXML header. + static List speciesNames(Path orthoxml) { + def root = parseOrthoxml(orthoxml) + def species = root.getElementsByTagNameNS('*', 'species') + return (0.. elementCounts(Path orthoxml) { def root = parseOrthoxml(orthoxml) return ['orthologGroup', 'paralogGroup', 'gene', 'geneRef'].collectEntries { name -> diff --git a/nf-tests/options.nf.test b/nf-tests/options.nf.test index 398e290..3cf9e85 100644 --- a/nf-tests/options.nf.test +++ b/nf-tests/options.nf.test @@ -77,6 +77,33 @@ nextflow_pipeline { } } + test("archive made with macOS tar (AppleDouble ._ files)") { + + when { + params { + outdir = "$outputDir" + output_folder = "$outputDir" + hogmap_folder = "$outputDir/hogmap" + // same content as input_archive.tar.gz, plus a ._ metadata file for every entry + input = "${projectDir}/nf-tests/helpers/input_archive_macos.tar.gz" + } + } + + then { + def out = file(params.output_folder).toPath() + def hogmap_files = file("${params.hogmap_folder}").listFiles() + def species = FastomaOutput.speciesNames(out.resolve('FastOMA_HOGs.orthoxml')) + + assertAll( + { assert workflow.success }, + // the hidden files are neither taken as proteomes nor passed to omamer + { assert hogmap_files.findAll { it.name.startsWith('.') }.isEmpty() }, + { assert hogmap_files.findAll { it.name.endsWith('.hogmap') }.size() == 3 }, + { assert species == ['AQUAE', 'CHLTR', 'MYCGE'] }, + ) + } + } + test("splice files select one isoform per gene") { when { From 6da5ed82b8ecffcca79c50702bc14096e8cfb331 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Wed, 4 Nov 2026 20:49:14 +0100 Subject: [PATCH 5/5] fix: renaming internal nodes must not move supplied labels to other clades Generated names (internal_, _) only avoided leaf names and names of nodes visited before, so in (((A,B),C)internal_0,D)root; the unnamed (A,B) clade took internal_0 and the clade that was called internal_0 was renamed. Original names of all internal nodes are now reserved up front, also for names cleaned of special characters. The tests assert the clade-to-name mapping. Refs #118 (review comment on #125) Co-Authored-By: Claude Opus 5.5 --- FastOMA/check_input.py | 16 ++++++++++++---- tests/test_check_input.py | 13 +++++++++---- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/FastOMA/check_input.py b/FastOMA/check_input.py index 3659329..259dc0b 100644 --- a/FastOMA/check_input.py +++ b/FastOMA/check_input.py @@ -163,16 +163,21 @@ def check_omamer_db(omamerdb_adress=None): def add_internal_node_prune(species_tree, species_names, out_fname): # add name for the internal, if no name is provided, removintg special chars counter_internal = 0 - # leaf names count as taken: an internal node named like a species would make the names - # ambiguous (pruning by name fails, and per-level pickle files are named after the nodes) + # names already given in the final tree; leaf names count as taken: an internal node named like a + # species would make the names ambiguous (pruning by name fails, and per-level pickle files are named + # after the nodes) node_names = set(species_tree.get_leaf_names()) + # generated names must also avoid the original names of nodes visited later, otherwise a supplied + # label would end up on a different clade + reserved = node_names | {n.name for n in species_tree.traverse() if n.name} def numbered(prefix): nonlocal counter_internal while True: candidate = prefix + str(counter_internal) counter_internal += 1 - if candidate not in node_names: + if candidate not in node_names and candidate not in reserved: + reserved.add(candidate) return candidate for node in species_tree.traverse(strategy="postorder"): @@ -185,7 +190,10 @@ def numbered(prefix): logger.debug("The internal node name '%s' was too short, repeated or made of special chars, " "it is changed to '%s'", node_name, node.name) elif cleaned != node_name: - node.name = numbered(cleaned + "_") if cleaned in node_names else cleaned + # the cleaned name must not take another node's original name either (this node's own + # original name has special chars, so it can't be the one found in `reserved`) + node.name = numbered(cleaned + "_") if cleaned in node_names or cleaned in reserved else cleaned + reserved.add(node.name) logger.debug("The internal node name has special chars or repeated '%s' which is changed to '%s'", node_name, node.name) node_names.add(node.name) diff --git a/tests/test_check_input.py b/tests/test_check_input.py index 6eac191..a775502 100644 --- a/tests/test_check_input.py +++ b/tests/test_check_input.py @@ -103,11 +103,16 @@ def prune(tmp_path, nwk, species): def test_rename_and_prune(self, tmp_path, nwk, species, expected): assert self.prune(tmp_path, nwk, list(species)) == expected - def test_generated_names_do_not_clash_with_existing(self, tmp_path): - # internal_0 is already used further up the tree, so the unnamed node gets internal_1 + def test_generated_names_do_not_take_existing_labels(self, tmp_path): + # internal_0 is the supplied name of the (A,B,C) clade, visited after (A,B): the unnamed (A,B) + # clade must get another name, and (A,B,C) must keep its label res = self.prune(tmp_path, "(((A,B),C)internal_0,D)root;", list("ABCD")) - names = [n.name for n in Tree(res, format=1).traverse() if not n.is_leaf()] - assert sorted(names) == ["internal_0", "internal_1", "root"] + assert res == "(((A:1,B:1)internal_1:1,C:1)internal_0:1,D:1)root:0;" + + def test_cleaned_names_do_not_take_existing_labels(self, tmp_path): + # 'a.bc' cleans to abc, which is the supplied name of a clade visited later + res = self.prune(tmp_path, "((A,B)'a.bc',(C,D)abc)root;", list("ABCD")) + assert res == "((A:1,B:1)abc_0:1,(C:1,D:1)abc:1)root:0;" def test_name_too_short_after_removing_special_chars(self, tmp_path): res = self.prune(tmp_path, "((A,B)'a.b',(C,D)cde)root;", list("ABCD"))