Tcs and implied loss scores

FastOMA/_hog_class.py

@@ -23,6 +23,178 @@
 from . import _utils_subhog, logger
 
 
+def _member_species(members) -> set:
+    return set(m.split("||")[1] for m in members)
+
+
+def _species_name_index(tree: TreeNode) -> dict:
+    """species name -> TreeNode lookup for `tree`, built once and cached on its root.
+
+    ete3's search_nodes()/get_common_ancestor(names) each do a full linear scan of the tree;
+    at scale (e.g. a rootHOG with 100k+ genes, each needing a per-species lookup) calling those
+    once per gene dominates runtime. Building this index once per species tree and doing O(1)
+    dict lookups instead is ~1000x faster in practice and is safe to cache on the tree: a fresh
+    species tree object is built per rootHOG (see prepare_species_tree), so the cache can't go
+    stale across rootHOGs, and it's built lazily on first use so callers never pass it explicitly."""
+    root = tree.get_tree_root()
+    index = getattr(root, "_species_name_index", None)
+    if index is None:
+        index = {n.name: n for n in root.traverse()}
+        root.add_feature("_species_name_index", index)
+    return index
+
+
+def _species_mrca(tree: TreeNode, species_of_members: set) -> TreeNode:
+    """Common ancestor of `species_of_members` within `tree`.
+
+    `tree` should be the real, unpruned species tree whenever one is available: the working
+    species (sub)tree used during inference is pruned to the species present in the rootHOG
+    (see `_utils_subhog.prepare_species_tree`), so a species entirely absent from the rootHOG
+    no longer exists in it, and computing losses against that tree would silently ignore it.
+    A node that genuinely branches for `species_of_members` is never affected by that pruning
+    (pruning only ever removes nodes that become unary), so this mrca has the same identity/name
+    whether computed against the pruned or the unpruned tree -- only the topology below it, which
+    ImpliedLosses/TCSScore need, differs."""
+    index = _species_name_index(tree)
+    return tree.get_common_ancestor(*[index[x] for x in species_of_members])
+
+
+def _count_implied_losses(mrca: TreeNode, species_with_members: set) -> int:
+    """Dollo-parsimony style count of minimal loss events between `mrca` and the species
+    that actually have a member gene: a whole clade lacking any member counts as one loss,
+    regardless of how many species it spans."""
+    def recurse(node):
+        leaves_under = set(n.name for n in node.iter_leaves())
+        if leaves_under.isdisjoint(species_with_members):
+            return 1
+        if node.is_leaf():
+            return 0
+        return sum(recurse(c) for c in node.children)
+    return sum(recurse(c) for c in mrca.children)
+
+
+def _count_untouched_clades(tree_node: TreeNode, reached_nodes: list) -> int:
+    """Counts, below `tree_node`, the maximal clades that contain none of `reached_nodes` --
+    i.e. clades no sub-hog lineage touches at all -- as one loss each. Does not descend into a
+    clade that a lineage already reaches: that clade's internal losses (if any) are counted
+    separately, per lineage, by _hog_implied_losses -- re-scanning them here would double-count."""
+    if tree_node in reached_nodes:
+        return 0
+    if not any(tree_node in r.get_ancestors() for r in reached_nodes):
+        return 1
+    return sum(_count_untouched_clades(c, reached_nodes) for c in tree_node.children)
+
+
+def _hog_implied_losses(hog: "HOG", tree_node: TreeNode) -> int:
+    """Recursively counts implied gene losses for `hog` within `tree_node`'s subtree.
+
+    Mirrors _tax_overlap's _tax_now-based grouping of _subhogs, but sums losses per group member
+    instead of intersecting/folding lineages: a group of sub-hogs sharing the same _tax_now
+    represents a duplication, and each such paralogous copy's losses are counted independently
+    so that one copy's presence in a species can't mask a sibling copy's loss there."""
+    subhogs = getattr(hog, "_subhogs", None)
+    if not subhogs:
+        return 0
+
+    losses = 0
+    clade_groups = {}
+    for sh in subhogs:
+        clade_groups.setdefault(sh._tax_now.name, []).append(sh)
+    for subhogs_of_clade in clade_groups.values():
+        clade_node = subhogs_of_clade[0]._tax_now
+        for sh in subhogs_of_clade:
+            losses += _hog_implied_losses(sh, clade_node)
+
+    reached_nodes = [subhogs_of_clade[0]._tax_now for subhogs_of_clade in clade_groups.values()]
+    losses += _count_untouched_clades(tree_node, reached_nodes)
+    return losses
+
+
+def _species_lineage_index(tree: TreeNode) -> dict:
+    """species name -> frozenset of ancestor node names (inclusive of the species itself and the
+    tree's root), cached on the tree's root like _species_name_index.
+
+    Unlike _species_mrca (which needs a lineage relative to some already-known mrca), the TCS
+    taxonomy-overlap score below wants each species' *absolute* lineage: it discovers the HOG's
+    own mrca as a side effect of intersecting members' lineages, rather than needing it supplied
+    up front (see attach_scores)."""
+    root = tree.get_tree_root()
+    index = getattr(root, "_species_lineage_index", None)
+    if index is None:
+        index = {}
+        for leaf in root.iter_leaves():
+            ancestors = set()
+            n = leaf
+            while True:
+                ancestors.add(n.name)
+                if n is root:
+                    break
+                n = n.up
+            index[leaf.name] = frozenset(ancestors)
+        root.add_feature("_species_lineage_index", index)
+    return index
+
+
+def _combine_tax_overlap(parts):
+    """Folds sibling (nset, leaf_size, leaf_acc, tax_score) tuples into their parent's, per
+    Moi/Kim's taxonomy-overlap algorithm: a "match" (nonempty lineage shared by every
+    contributing part) earns len(nset) points per gene (leaf_size) below this node, on top of
+    whatever each part already scored deeper down. A part with an empty nset (no informative
+    match anywhere in its own subtree) is excluded from the intersection rather than zeroing it
+    out -- one already-incongruent branch shouldn't poison a genuinely-matching sibling."""
+    leaf_size = sum(p[1] for p in parts)
+    leaf_acc = sum(p[2] for p in parts) + leaf_size
+    tax_score = sum(p[3] for p in parts)
+    nonempty = [p[0] for p in parts if p[0]]
+    nset = frozenset.intersection(*nonempty) if nonempty else frozenset()
+    tax_score += len(nset) * leaf_size
+    return nset, leaf_size, leaf_acc, tax_score
+
+
+def _tax_overlap(hog: "HOG", species_lineage_index: dict):
+    """Recursively computes (nset, leaf_size, leaf_acc, tax_score) for `hog`, mirroring the same
+    _tax_now-based grouping of _subhogs used by HOG.to_orthoxml() (a group of sub-hogs sharing a
+    _tax_now is a duplication -- its members are folded together as siblings, same as a group of
+    literal gene-tree children would be). `leaf_size` is read off len(hog).
+    See attach_scores for how the four returned values become TCSScore."""
+    if not hog._subhogs:
+        species = next(iter(_member_species(hog.get_members())))
+        return species_lineage_index[species], len(hog), 0, 0
+
+    groups = []
+    for _, subhogs_of_clade in itertools.groupby(
+            sorted(hog._subhogs, key=lambda h: h._tax_now.name), key=lambda h: h._tax_now.name):
+        subhogs_of_clade = list(subhogs_of_clade)
+        if len(subhogs_of_clade) == 1:
+            groups.append(_tax_overlap(subhogs_of_clade[0], species_lineage_index))
+        else:
+            groups.append(_combine_tax_overlap([_tax_overlap(sh, species_lineage_index) for sh in subhogs_of_clade]))
+    return groups[0] if len(groups) == 1 else _combine_tax_overlap(groups)
+
+
+def attach_scores(hog_element: ET.Element, hog: "HOG", mrca: TreeNode, species_of_members: set) -> None:
+    """Computes and attaches CompletenessScore, TCSScore and ImpliedLosses as <score>
+    sub-elements of `hog_element`.
+
+    TCSScore follows Moi et al. 2025 / Kim et al. 2026's taxonomy-overlap score: _tax_overlap()
+    computes it using each species' absolute lineage (not one relative to `mrca`), and the HOG's
+    own mrca-relative "ideal" contribution is discovered and subtracted algebraically at the end
+    (tax_score - leaf_acc * len(nset)) rather than needing `mrca` supplied up front, then
+    normalized by leaf_size -- the gene count of `hog` itself."""
+    completeness_score = round(len(species_of_members) / mrca.size, 4)
+    ET.SubElement(hog_element, "score", attrib={"id": "CompletenessScore", "value": str(completeness_score)})
+
+    if getattr(hog, "_subhogs", None):
+        implied_losses = _hog_implied_losses(hog, mrca)
+    else:
+        implied_losses = _count_implied_losses(mrca, species_of_members)
+    ET.SubElement(hog_element, "score", attrib={"id": "ImpliedLosses", "value": str(implied_losses)})
+
+    nset, leaf_size, leaf_acc, tax_score = _tax_overlap(hog, _species_lineage_index(mrca))
+    tcs_score = (tax_score - leaf_acc * len(nset)) / leaf_size
+    ET.SubElement(hog_element, "score", attrib={"id": "TCSScore", "value": str(round(tcs_score, 4))})
+
+
 # from .infer_subhogs import conf_infer_subhhogs #fastoma_infer_subhogs #
 
 
@@ -257,7 +429,7 @@ def merge_prots_msa(self, merged_fragment_name, merged_msa_new):   # merged_frag
     #     self._msa = MultipleSeqAlignment(msa_new)
     #     return 1
 
-    def to_orthoxml(self):
+    def to_orthoxml(self, full_species_tree: Optional[TreeNode] = None):
         if len(self._subhogs) == 0:
             list_member = list(self._members)
             if len(list_member) == 1:
@@ -313,7 +485,7 @@ def _sorter_key(sh):
                 # the following line could be improved, instead of tax_now we can use the least common ancestor of all members
                 # property_element = ET.SubElement(paralog_element, "property",attrib={"name": "TaxRange", "value": str(sub_clade)}) # self._tax_now
                 for sh in list_of_subhogs_of_same_clade:
-                    element_p = sh.to_orthoxml()
+                    element_p = sh.to_orthoxml(full_species_tree)
                     if str(element_p):
                         paralog_element.append(element_p)  # ,**gene_id_name  indent+2
                     else:
@@ -323,7 +495,7 @@ def _sorter_key(sh):
             elif len(list_of_subhogs_of_same_clade) == 1:
                 subhog = list_of_subhogs_of_same_clade[0]
                 if len(subhog._members):
-                    element = subhog.to_orthoxml()
+                    element = subhog.to_orthoxml(full_species_tree)
                     if str(element):  # element could be  <Element 'geneRef' at 0x7f7f9bacb450>
                         element_list.append(element)  # indent+2
                     else:
@@ -337,16 +509,13 @@ def _sorter_key(sh):
         elif len(element_list) > 1:
             #hog_elemnt = ET.Element('orthologGroup', attrib={"id": str(self._hogid)})
             hog_elemnt = ET.Element('orthologGroup', attrib={"id": str(self._hogid)}, )
-            species_of_members = set([i.split("||")[1] for i in self._members])  #  'tr|H2MU14|H2MU14_ORYLA||ORYLA||1056022282'
-            num_species_tax_hog = len(species_of_members)
-            mrca = self.taxlevel.get_common_ancestor(
-                *[self.taxlevel.search_nodes(name=x)[0] for x in species_of_members])
-            if mrca != self.taxlevel:
+            species_of_members = _member_species(self._members)  #  'tr|H2MU14|H2MU14_ORYLA||ORYLA||1056022282'
+            mrca = _species_mrca(full_species_tree if full_species_tree is not None else self.taxlevel, species_of_members)
+            if mrca.name != self.taxlevel.name:
                 logger.info(f"mrca ({mrca.name}) != self.taxlevel ({self.taxlevel.name})")
                 logger.info(f"<{hog_elemnt.tag} {hog_elemnt.attrib}>")
 
-            completeness_score = round(num_species_tax_hog/mrca.size, 4)
-            property_element = ET.SubElement(hog_elemnt, "score", attrib={"id": "CompletenessScore", "value": str(completeness_score)})
+            attach_scores(hog_elemnt, self, mrca, species_of_members)
             property_element = ET.SubElement(hog_elemnt, "property", attrib={"name": "TaxRange", "value": str(mrca.name)})
 
             for element in element_list:

FastOMA/_infer_subhog.py

@@ -26,7 +26,7 @@
 from . import _wrappers, logger
 from . import _utils_subhog
 from . import _utils_frag_SO_detection
-from ._hog_class import HOG, Representative, split_hog
+from ._hog_class import HOG, Representative, split_hog, attach_scores, _member_species, _species_name_index
 from ._utils_subhog import MSAFilter, MSAFilterElbow, MSAFilterTrimAL
 
 from .zoo.utils import unique
@@ -73,7 +73,7 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder,  pi
     # the file "species_tree_checked.nwk" is created by the check_input.py
     (species_tree) = _utils_subhog.read_species_tree(conf_infer_subhhogs.species_tree)
 
-    (species_tree, species_names_rhog, prot_names_rhog) = _utils_subhog.prepare_species_tree(rhog_i, species_tree, rhogid)
+    (species_tree, species_names_rhog, prot_names_rhog, full_species_tree) = _utils_subhog.prepare_species_tree(rhog_i, species_tree, rhogid)
     species_names_rhog = list(set(species_names_rhog))
     logger.info("Number of unique species in rHOG " + rhogid + " is " + str(len(species_names_rhog)) + ".")
 
@@ -102,14 +102,15 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder,  pi
         tot_genes += len(hog_i)
         if len(hog_i) >= inferhog_min_hog_size_xml:
             # could be improved   # hogs_a_rhog_xml = hog_i.to_orthoxml(**gene_id_name)
-            hogs_a_rhog_xml_raw = hog_i.to_orthoxml()    # <generef  >      <paralg object >
+            hogs_a_rhog_xml_raw = hog_i.to_orthoxml(full_species_tree)    # <generef  >      <paralg object >
             if orthoxml_v03 and 'paralogGroup' in str(hogs_a_rhog_xml_raw) :
                 # in version v0.3 of orthoxml, there shouldn't be any paralogGroup at root level. Let's put them inside an orthogroup should be in
                 hog_elemnt = ET.Element('orthologGroup', attrib={"id": str(hog_i.hogid)})
-                num_species_tax_hog = len(set([i.split("||")[1] for i in hog_i.get_members()]))
-                completeness_score = round(num_species_tax_hog / hog_i.taxlevel.size, 4)
-                ET.SubElement(hog_elemnt, "score",
-                              attrib={"id": "CompletenessScore", "value": str(completeness_score)})
+                species_of_members = _member_species(hog_i.get_members())
+                scoring_node = hog_i.taxlevel
+                if full_species_tree is not None:
+                    scoring_node = _species_name_index(full_species_tree)[hog_i.taxlevel.name]
+                attach_scores(hog_elemnt, hog_i, scoring_node, species_of_members)
                 ET.SubElement(hog_elemnt, "property", attrib={"name": "TaxRange", "value": str(hog_i.taxname)})
                 hog_elemnt.append(hogs_a_rhog_xml_raw)
                 hogs_a_rhog_xml = hog_elemnt
@@ -163,6 +164,10 @@ def build_xml_from_rhog(rhogid:str, seqs:List[SeqRecord], hogs:List[ET.Element])
     scores = ET.SubElement(root, "scores")
     ET.SubElement(scores, "scoreDef",
                   {"id": "CompletenessScore","desc": "Fraction of expected species with genes in the (Sub)HOG"})
+    ET.SubElement(scores, "scoreDef",
+                  {"id": "TCSScore", "desc": "Taxonomic Congruence Score: how well the (Sub)HOG structure matches the species tree topology"})
+    ET.SubElement(scores, "scoreDef",
+                  {"id": "ImpliedLosses", "desc": "Number of implied gene loss events (Dollo parsimony) within the (Sub)HOG's taxonomic range"})
     groups = ET.SubElement(root, 'groups')
     for hog in hogs:
         groups.append(hog)

FastOMA/_utils_subhog.py

@@ -253,7 +253,10 @@ def prepare_species_tree(rhog_i: List[SeqRecord], species_tree: Tree, rhogid: st
     orthoxml_to_newick.py function for extracting orthoxml_to_newick.py subtree from the input species tree  orthoxml_to_newick.py.k.orthoxml_to_newick.py pruning,
     based on the names of species in the rootHOG.
 
-    output: species_tree (pruned), species_names_rhog, prot_names_rhog
+    output: species_tree (pruned), species_names_rhog, prot_names_rhog, full_species_tree (unpruned,
+    with the `size` feature set) -- kept around so that scores computed later (e.g. ImpliedLosses,
+    TCSScore) can be evaluated against the real species tree topology, not just the species present
+    in this rootHOG.
     """
     assert len(rhog_i) > 0, 'input hog_i is empty, probably previous step find_rhog has issue, rhogs/HOG_B0'+rhogid+'is empty?'
     species_names_rhog = []
@@ -275,9 +278,11 @@ def prepare_species_tree(rhog_i: List[SeqRecord], species_tree: Tree, rhogid: st
     for n in species_tree.traverse():
         n.add_feature("size", len(n))
 
+    full_species_tree = species_tree.copy()
+
     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
+    return mrca, species_names_rhog, prot_names_rhog, full_species_tree
 
 
 def label_sd_internal_nodes(tree_out, threshold_dubious_sd):

FastOMA/collect_subhogs.py

@@ -102,7 +102,14 @@ def _annotateGroupR(node: ET.ElementTree, og: str, idx: int = 0):
 
     omamer_roothog_id = ":".join(hog.get('id').split("_")[0:2])
     fam_elem = ET.Element("property", {"name": "OMAmerRootHOG", "value": omamer_roothog_id})
-    hog.insert(1, fam_elem)
+    # orthoxml requires all <score> children before any <property> children, so insert
+    # the new property right after the trailing run of <score> elements, not at a fixed index.
+    insert_idx = 0
+    for child in hog:
+        if child.tag != "score":
+            break
+        insert_idx += 1
+    hog.insert(insert_idx, fam_elem)
     _annotateGroupR(hog, "HOG:{:07d}".format(fam))
     return hog
 
@@ -189,6 +196,10 @@ def write_hog_orthoxml(pickle_folder, output_xml_name, gene_id_pickle_file, id_t
     scores = ET.SubElement(orthoxml_file, "scores")
     ET.SubElement(scores, "scoreDef", {"id": "CompletenessScore",
                                        "desc": "Fraction of expected species with genes in the (Sub)HOG"})
+    ET.SubElement(scores, "scoreDef", {"id": "TCSScore",
+                                       "desc": "Taxonomic Congruence Score: how well the (Sub)HOG structure matches the species tree topology"})
+    ET.SubElement(scores, "scoreDef", {"id": "ImpliedLosses",
+                                       "desc": "Number of implied gene loss events (Dollo parsimony) within the (Sub)HOG's taxonomic range"})
 
     #  #### create the groups of orthoxml   ####
     groups_xml = ET.SubElement(orthoxml_file, "groups")

FastOMA/zoo/wrappers/treebuilders/fasttree.py

@@ -94,7 +94,7 @@ def _call(self, filename, *args, **kwargs):
             self.stderr = self.cli.get_stderr()
             last_error_line = self.stderr.split('\n')[-1].strip()
             logger.error('FastTree returned non-zero exit status: {}'.format(self.returncode))
-            logger.error('Output of FastTree:\n\n%s\nstdout=\n%s\n{}\n\n%s\nstderr=\n%s\n{}\n\n',
+            logger.error('Output of FastTree:\n\n%s\nstdout=\n%s\n%s\n\n%s\nstderr=\n%s\n%s\n\n',
                          "=" * 30, "=" * 30, summarize_long_message(self.stdout),
                          "=" * 30, "=" * 30, summarize_long_message(self.stderr))
             if self.returncode < 0:

nf-tests/default.nf.test.snap

@@ -50,7 +50,7 @@
             [
                 {
                     "name": "FastOMA_HOGs.orthoxml",
-                    "lineCount": 162
+                    "lineCount": 194
                 },
                 {
                     "name": "OrthologousGroups.tsv",
@@ -80,8 +80,8 @@
         ],
         "meta": {
             "nf-test": "0.9.3",
-            "nextflow": "24.10.5"
+            "nextflow": "26.04.6"
         },
-        "timestamp": "2025-10-23T12:09:08.301966"
+        "timestamp": "2026-08-25T18:07:32.287389"
     }
 }

tests/data/tcs_taxonomy.tsv

@@ -0,0 +1,5 @@
+query	lineage
+SP1	1 (class), 2 (order), 3 (family), 4 (genus), 5 (species)
+SP2	1 (class), 2 (order), 3 (family), 4 (genus), 6 (species)
+SP3	1 (class), 7 (order), 8 (family), 9 (genus), 10 (species)
+SP4	1 (class), 7 (order), 11 (family), 12 (subfamily), 13 (genus), 14 (species), 15 (subspecies)