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/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..6beff51 100644 --- a/FastOMA/_utils_roothog.py +++ b/FastOMA/_utils_roothog.py @@ -24,128 +24,24 @@ Gene = collections.namedtuple("Gene", ("numeric_id", "prot_id", "main_isoform"), defaults=(None,) ) -"""UnionFind.py +FASTA_EXTENSIONS = ("fa", "fasta") -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. -""" +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(".") -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: +def list_proteome_files(folder): + """(species name, extension) of the proteome FASTA files in folder. - - 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 + 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 @@ -160,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 2fa4f07..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 @@ -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/check_input.py b/FastOMA/check_input.py index e7d6c2b..259dc0b 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,40 @@ 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() + # 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 and candidate not in reserved: + reserved.add(candidate) + 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: + # 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_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 +214,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 +253,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/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/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/nf-tests/helpers/input_archive_macos.tar.gz b/nf-tests/helpers/input_archive_macos.tar.gz new file mode 100644 index 0000000..5c5a7f0 Binary files /dev/null and b/nf-tests/helpers/input_archive_macos.tar.gz differ 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 f7b8f43..3cf9e85 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 }, @@ -75,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 { 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_check_input.py b/tests/test_check_input.py index 05ddc0c..a775502 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,25 @@ 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_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")) + 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")) + 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 +144,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_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_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"]) 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 b35b7f0..2868517 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" @@ -61,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): 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