Fix small internal bugs, input validation and report without CompletenessScore

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 ._<name> 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'))
         } 
     }
 

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)
 

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)
 
 

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 ._<name> 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, ..

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
 
 

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: