From eb40beed7d8fbcf83e9a983396671524796bfd49 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Tue, 18 Aug 2026 13:23:45 +0200 Subject: [PATCH 1/4] [WIP] adding initial implementation of add scores adds TCS score according to reference implementation. adds implied loss score computing the nr of implied loss of a (sub)HOG. --- FastOMA/_hog_class.py | 86 ++++++++++++++++++++++++-- FastOMA/_infer_subhog.py | 12 ++-- FastOMA/collect_subhogs.py | 4 ++ tests/test_hog_scores.py | 123 +++++++++++++++++++++++++++++++++++++ 4 files changed, 216 insertions(+), 9 deletions(-) create mode 100644 tests/test_hog_scores.py diff --git a/FastOMA/_hog_class.py b/FastOMA/_hog_class.py index f1712be..fa263d3 100644 --- a/FastOMA/_hog_class.py +++ b/FastOMA/_hog_class.py @@ -23,6 +23,86 @@ from . import _utils_subhog, logger +def _member_species(members) -> set: + return set(m.split("||")[1] for m in 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 _induced_congruent_score(node: TreeNode, species_with_members: set, depth: int = 0): + """Best achievable TCS-style score for `species_with_members` given the real species + tree topology below `node`, i.e. the score a perfectly congruent tree would obtain. + + Returns (has_member, score) where `has_member` indicates whether any leaf under `node` + is in `species_with_members`.""" + if node.is_leaf(): + return node.name in species_with_members, 0.0 + results = [_induced_congruent_score(c, species_with_members, depth + 1) for c in node.children] + total = sum(s for _, s in results) + n_contributing = sum(1 for has, _ in results if has) + if n_contributing >= 2: + total += depth + return n_contributing > 0, total + + +def _combine_tcs_parts(parts): + lineages = [lin for lin, _ in parts if lin is not None] + score = sum(s for _, s in parts) + lineage = frozenset.intersection(*lineages) if lineages else None + return lineage, score + (len(lineage) if lineage is not None else 0) + + +def _tcs_partial(hog: "HOG", mrca: TreeNode): + """Recursively computes (lineage relative to `mrca`, raw TCS score) for `hog`, mirroring + the same _tax_now-based grouping of _subhogs used by HOG.to_orthoxml().""" + if not hog._subhogs: + species = next(iter(_member_species(hog.get_members()))) + lineage = set() + n = mrca.search_nodes(name=species)[0] + while n is not mrca: + lineage.add(n.name) + n = n.up + return frozenset(lineage), 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(_tcs_partial(subhogs_of_clade[0], mrca)) + else: + groups.append(_combine_tcs_parts([_tcs_partial(sh, mrca) for sh in subhogs_of_clade])) + return _combine_tcs_parts(groups) if len(groups) > 1 else groups[0] + + +def attach_scores(hog_element: ET.Element, hog: "HOG", mrca: TreeNode, species_of_members: set) -> None: + """Computes and attaches CompletenessScore, TCSScore (if defined) and ImpliedLosses + as sub-elements of `hog_element`.""" + completeness_score = round(len(species_of_members) / mrca.size, 4) + ET.SubElement(hog_element, "score", attrib={"id": "CompletenessScore", "value": str(completeness_score)}) + + implied_losses = _count_implied_losses(mrca, species_of_members) + ET.SubElement(hog_element, "score", attrib={"id": "ImpliedLosses", "value": str(implied_losses)}) + + _, ideal_score = _induced_congruent_score(mrca, species_of_members) + if ideal_score > 0: + _, raw_score = _tcs_partial(hog, mrca) + tcs_score = round(raw_score / ideal_score, 4) + ET.SubElement(hog_element, "score", attrib={"id": "TCSScore", "value": str(tcs_score)}) + + # from .infer_subhogs import conf_infer_subhhogs #fastoma_infer_subhogs # @@ -337,16 +417,14 @@ 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) + species_of_members = _member_species(self._members) # 'tr|H2MU14|H2MU14_ORYLA||ORYLA||1056022282' mrca = self.taxlevel.get_common_ancestor( *[self.taxlevel.search_nodes(name=x)[0] for x in species_of_members]) if mrca != self.taxlevel: 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: diff --git a/FastOMA/_infer_subhog.py b/FastOMA/_infer_subhog.py index d098acb..160d48f 100644 --- a/FastOMA/_infer_subhog.py +++ b/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 from ._utils_subhog import MSAFilter, MSAFilterElbow, MSAFilterTrimAL from .zoo.utils import unique @@ -106,10 +106,8 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi 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()) + attach_scores(hog_elemnt, hog_i, hog_i.taxlevel, 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 +161,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) diff --git a/FastOMA/collect_subhogs.py b/FastOMA/collect_subhogs.py index 26c7fb4..6314f19 100644 --- a/FastOMA/collect_subhogs.py +++ b/FastOMA/collect_subhogs.py @@ -189,6 +189,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") diff --git a/tests/test_hog_scores.py b/tests/test_hog_scores.py new file mode 100644 index 0000000..5296d25 --- /dev/null +++ b/tests/test_hog_scores.py @@ -0,0 +1,123 @@ +import xml.etree.ElementTree as ET +from unittest import TestCase + +from ete3 import Tree + +from FastOMA._hog_class import ( + HOG, + attach_scores, + _member_species, + _count_implied_losses, + _induced_congruent_score, + _tcs_partial, +) + + +def leaf_hog(prot_id, species_node): + hog = HOG.__new__(HOG) + hog._rhogid = "test" + hog._hogid = "HOG_test_" + prot_id + hog._tax_now = species_node + hog._members = {prot_id} + hog._subhogs = [] + hog._dubious_members = set() + return hog + + +def merged_hog(subhogs, tax_node): + hog = HOG.__new__(HOG) + hog._rhogid = "test" + hog._hogid = "HOG_test_" + tax_node.name + hog._tax_now = tax_node + hog._subhogs = list(subhogs) + hog._members = set().union(*(sh._members for sh in subhogs)) + hog._dubious_members = set() + return hog + + +def score_value(elem, score_id): + for score in elem.findall("score"): + if score.get("id") == score_id: + return score.get("value") + return None + + +class BalancedSpeciesTreeTests(TestCase): + """Species tree: M(A(sp1,sp2), B(sp3,sp4))""" + + def setUp(self): + self.sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + for n in self.sptree.traverse(): + n.add_feature("size", len(n)) + self.M = self.sptree + self.A = self.sptree.search_nodes(name="A")[0] + self.B = self.sptree.search_nodes(name="B")[0] + + def _make_full_hog(self): + leaf1 = leaf_hog("p1||sp1", self.sptree.search_nodes(name="sp1")[0]) + leaf2 = leaf_hog("p2||sp2", self.sptree.search_nodes(name="sp2")[0]) + leaf3 = leaf_hog("p3||sp3", self.sptree.search_nodes(name="sp3")[0]) + leaf4 = leaf_hog("p4||sp4", self.sptree.search_nodes(name="sp4")[0]) + hogA = merged_hog([leaf1, leaf2], self.A) + hogB = merged_hog([leaf3, leaf4], self.B) + return merged_hog([hogA, hogB], self.M) + + def test_fully_congruent_hog(self): + hogM = self._make_full_hog() + species_of_members = _member_species(hogM.get_members()) + self.assertEqual(species_of_members, {"sp1", "sp2", "sp3", "sp4"}) + + elem = ET.Element("orthologGroup") + attach_scores(elem, hogM, self.M, species_of_members) + + self.assertEqual(score_value(elem, "CompletenessScore"), "1.0") + self.assertEqual(score_value(elem, "ImpliedLosses"), "0") + self.assertEqual(score_value(elem, "TCSScore"), "1.0") + + def test_clean_loss_does_not_penalize_tcs(self): + leaf1 = leaf_hog("p1||sp1", self.sptree.search_nodes(name="sp1")[0]) + leaf3 = leaf_hog("p3||sp3", self.sptree.search_nodes(name="sp3")[0]) + leaf4 = leaf_hog("p4||sp4", self.sptree.search_nodes(name="sp4")[0]) + hogA = merged_hog([leaf1], self.A) # sp2 lost, single sub-hog passed through + hogB = merged_hog([leaf3, leaf4], self.B) + hogM = merged_hog([hogA, hogB], self.M) + + species_of_members = _member_species(hogM.get_members()) + self.assertEqual(species_of_members, {"sp1", "sp3", "sp4"}) + + elem = ET.Element("orthologGroup") + attach_scores(elem, hogM, self.M, species_of_members) + + self.assertEqual(score_value(elem, "CompletenessScore"), "0.75") + self.assertEqual(score_value(elem, "ImpliedLosses"), "1") + self.assertEqual(score_value(elem, "TCSScore"), "1.0") + + def test_direct_sibling_species_have_no_tcs_score(self): + leaf1 = leaf_hog("p1||sp1", self.sptree.search_nodes(name="sp1")[0]) + leaf2 = leaf_hog("p2||sp2", self.sptree.search_nodes(name="sp2")[0]) + hogA = merged_hog([leaf1, leaf2], self.A) + + species_of_members = _member_species(hogA.get_members()) + elem = ET.Element("orthologGroup") + attach_scores(elem, hogA, self.A, species_of_members) + + self.assertIsNone(score_value(elem, "TCSScore")) + self.assertIsNotNone(score_value(elem, "CompletenessScore")) + self.assertIsNotNone(score_value(elem, "ImpliedLosses")) + + +class ImpliedLossesUnitTests(TestCase): + def test_whole_missing_clade_counts_as_one_loss(self): + sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + losses = _count_implied_losses(sptree, {"sp3", "sp4"}) + self.assertEqual(losses, 1) + + def test_no_losses_when_all_present(self): + sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + losses = _count_implied_losses(sptree, {"sp1", "sp2", "sp3", "sp4"}) + self.assertEqual(losses, 0) + + def test_single_species_loss_within_a_clade(self): + sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + losses = _count_implied_losses(sptree, {"sp1", "sp3", "sp4"}) + self.assertEqual(losses, 1) From 585610f13d234b77a353e236cff06ae70af6ce0a Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Wed, 26 Aug 2026 07:47:47 +0200 Subject: [PATCH 2/4] compute TCSscore and ImpliedLosses as expected. adds unittests about the score calculation and update the nf-test output with the addtional score lines --- FastOMA/_hog_class.py | 189 ++++++--- FastOMA/_infer_subhog.py | 11 +- FastOMA/_utils_subhog.py | 9 +- FastOMA/collect_subhogs.py | 9 +- FastOMA/zoo/wrappers/treebuilders/fasttree.py | 2 +- nf-tests/default.nf.test.snap | 6 +- tests/data/tcs_taxonomy.tsv | 5 + tests/data/tcs_tree1.nwk | 1 + tests/data/tcs_tree2.nwk | 1 + tests/data/tcs_tree3.nwk | 1 + tests/test_hog_scores.py | 366 +++++++++++++++++- 11 files changed, 534 insertions(+), 66 deletions(-) create mode 100644 tests/data/tcs_taxonomy.tsv create mode 100644 tests/data/tcs_tree1.nwk create mode 100644 tests/data/tcs_tree2.nwk create mode 100644 tests/data/tcs_tree3.nwk diff --git a/FastOMA/_hog_class.py b/FastOMA/_hog_class.py index fa263d3..9e57ee0 100644 --- a/FastOMA/_hog_class.py +++ b/FastOMA/_hog_class.py @@ -27,6 +27,38 @@ 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, @@ -41,66 +73,126 @@ def recurse(node): return sum(recurse(c) for c in mrca.children) -def _induced_congruent_score(node: TreeNode, species_with_members: set, depth: int = 0): - """Best achievable TCS-style score for `species_with_members` given the real species - tree topology below `node`, i.e. the score a perfectly congruent tree would obtain. - - Returns (has_member, score) where `has_member` indicates whether any leaf under `node` - is in `species_with_members`.""" - if node.is_leaf(): - return node.name in species_with_members, 0.0 - results = [_induced_congruent_score(c, species_with_members, depth + 1) for c in node.children] - total = sum(s for _, s in results) - n_contributing = sum(1 for has, _ in results if has) - if n_contributing >= 2: - total += depth - return n_contributing > 0, total - - -def _combine_tcs_parts(parts): - lineages = [lin for lin, _ in parts if lin is not None] - score = sum(s for _, s in parts) - lineage = frozenset.intersection(*lineages) if lineages else None - return lineage, score + (len(lineage) if lineage is not None else 0) - - -def _tcs_partial(hog: "HOG", mrca: TreeNode): - """Recursively computes (lineage relative to `mrca`, raw TCS score) for `hog`, mirroring - the same _tax_now-based grouping of _subhogs used by HOG.to_orthoxml().""" +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()))) - lineage = set() - n = mrca.search_nodes(name=species)[0] - while n is not mrca: - lineage.add(n.name) - n = n.up - return frozenset(lineage), 0.0 + 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(_tcs_partial(subhogs_of_clade[0], mrca)) + groups.append(_tax_overlap(subhogs_of_clade[0], species_lineage_index)) else: - groups.append(_combine_tcs_parts([_tcs_partial(sh, mrca) for sh in subhogs_of_clade])) - return _combine_tcs_parts(groups) if len(groups) > 1 else groups[0] + 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 (if defined) and ImpliedLosses - as sub-elements of `hog_element`.""" + """Computes and attaches CompletenessScore, TCSScore and ImpliedLosses as + 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)}) - implied_losses = _count_implied_losses(mrca, species_of_members) + 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)}) - _, ideal_score = _induced_congruent_score(mrca, species_of_members) - if ideal_score > 0: - _, raw_score = _tcs_partial(hog, mrca) - tcs_score = round(raw_score / ideal_score, 4) - ET.SubElement(hog_element, "score", attrib={"id": "TCSScore", "value": str(tcs_score)}) + 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 # @@ -337,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: @@ -393,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: @@ -403,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_list.append(element) # indent+2 else: @@ -418,9 +510,8 @@ def _sorter_key(sh): #hog_elemnt = ET.Element('orthologGroup', attrib={"id": str(self._hogid)}) hog_elemnt = ET.Element('orthologGroup', attrib={"id": str(self._hogid)}, ) species_of_members = _member_species(self._members) # 'tr|H2MU14|H2MU14_ORYLA||ORYLA||1056022282' - mrca = self.taxlevel.get_common_ancestor( - *[self.taxlevel.search_nodes(name=x)[0] for x in species_of_members]) - if mrca != self.taxlevel: + 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}>") diff --git a/FastOMA/_infer_subhog.py b/FastOMA/_infer_subhog.py index 160d48f..93a4455 100644 --- a/FastOMA/_infer_subhog.py +++ b/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, attach_scores, _member_species +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,12 +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() # + hogs_a_rhog_xml_raw = hog_i.to_orthoxml(full_species_tree) # 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)}) species_of_members = _member_species(hog_i.get_members()) - attach_scores(hog_elemnt, hog_i, hog_i.taxlevel, species_of_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 diff --git a/FastOMA/_utils_subhog.py b/FastOMA/_utils_subhog.py index 423c826..2fa4f07 100644 --- a/FastOMA/_utils_subhog.py +++ b/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): diff --git a/FastOMA/collect_subhogs.py b/FastOMA/collect_subhogs.py index 6314f19..e21a6a3 100644 --- a/FastOMA/collect_subhogs.py +++ b/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 children before any children, so insert + # the new property right after the trailing run of 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 diff --git a/FastOMA/zoo/wrappers/treebuilders/fasttree.py b/FastOMA/zoo/wrappers/treebuilders/fasttree.py index 5aa4eee..a0de294 100644 --- a/FastOMA/zoo/wrappers/treebuilders/fasttree.py +++ b/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: diff --git a/nf-tests/default.nf.test.snap b/nf-tests/default.nf.test.snap index 402bf29..5131a83 100644 --- a/nf-tests/default.nf.test.snap +++ b/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" } } \ No newline at end of file diff --git a/tests/data/tcs_taxonomy.tsv b/tests/data/tcs_taxonomy.tsv new file mode 100644 index 0000000..10fa299 --- /dev/null +++ b/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) diff --git a/tests/data/tcs_tree1.nwk b/tests/data/tcs_tree1.nwk new file mode 100644 index 0000000..1de792c --- /dev/null +++ b/tests/data/tcs_tree1.nwk @@ -0,0 +1 @@ +((SP1,SP2),(SP3,SP4)); diff --git a/tests/data/tcs_tree2.nwk b/tests/data/tcs_tree2.nwk new file mode 100644 index 0000000..96f9841 --- /dev/null +++ b/tests/data/tcs_tree2.nwk @@ -0,0 +1 @@ +(((SP1,SP2),SP3),SP4); diff --git a/tests/data/tcs_tree3.nwk b/tests/data/tcs_tree3.nwk new file mode 100644 index 0000000..425b276 --- /dev/null +++ b/tests/data/tcs_tree3.nwk @@ -0,0 +1 @@ +((SP1,SP3),(SP2,SP4)); diff --git a/tests/test_hog_scores.py b/tests/test_hog_scores.py index 5296d25..9a69549 100644 --- a/tests/test_hog_scores.py +++ b/tests/test_hog_scores.py @@ -1,3 +1,4 @@ +import os import xml.etree.ElementTree as ET from unittest import TestCase @@ -6,10 +7,14 @@ from FastOMA._hog_class import ( HOG, attach_scores, + _species_mrca, + _species_name_index, + _species_lineage_index, _member_species, _count_implied_losses, - _induced_congruent_score, - _tcs_partial, + _hog_implied_losses, + _tax_overlap, + _combine_tax_overlap, ) @@ -42,6 +47,133 @@ def score_value(elem, score_id): return None +DATA_DIR = os.path.join(os.path.dirname(__file__), "data") + + +def load_species_tree_from_lineage_tsv(path): + """Builds a species tree from a lineage TSV (columns: species name, comma-separated + " ()" entries from root to species), e.g. tests/data/tcs_taxonomy.tsv.""" + nodes = {} + root = None + with open(path) as fh: + for line in fh: + species, _, lineage = line.rstrip("\n").partition("\t") + if not lineage or species == "query": + continue + parent = None + for entry in lineage.split(","): + node_id = entry.strip().split(" ", 1)[0] + node = nodes.get(node_id) + if node is None: + node = Tree(name=node_id) + nodes[node_id] = node + if parent is not None: + parent.add_child(node) + elif root is None: + root = node + parent = node + parent.name = species # rename the lineage's terminal node to the species name + for n in root.traverse(): + n.add_feature("size", len(n)) + return root + + +def load_hog_from_genetree(nwk_path, species_tree): + """Builds a nested HOG mirroring a gene tree's topology, treating that topology as + already-reconciled sub-HOG grouping: each internal HOG's taxonomic level is the species-tree + MRCA of the species under it.""" + gene_tree = Tree(nwk_path) + + def build(node): + if node.is_leaf(): + species = node.name + return leaf_hog(f"{species}_p||{species}", species_tree.search_nodes(name=species)[0]) + subhogs = [build(child) for child in node.children] + members = set().union(*(sh._members for sh in subhogs)) + tax_now = _species_mrca(species_tree, _member_species(members)) + return merged_hog(subhogs, tax_now) + + return build(gene_tree) + + +class SpeciesNameIndexCachingTests(TestCase): + """_species_name_index()/_species_mrca() replace ete3's search_nodes()/get_common_ancestor- + by-name (each an O(species tree size) linear scan -- costly once you're doing one lookup per + gene on a rootHOG with 100k+ genes) with a name->TreeNode dict, built once and cached on the + tree's root.""" + + def setUp(self): + self.sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + for n in self.sptree.traverse(): + n.add_feature("size", len(n)) + + def test_index_maps_every_name_to_the_right_node(self): + index = _species_name_index(self.sptree) + for name in ("M", "A", "B", "sp1", "sp2", "sp3", "sp4"): + node = index[name] + self.assertEqual(node.name, name) + self.assertIs(node, self.sptree.search_nodes(name=name)[0]) + + def test_index_is_built_once_and_cached_on_the_root(self): + first = _species_name_index(self.sptree) + second = _species_name_index(self.sptree) + self.assertIs(first, second) + # also reachable (and the same cache) starting from a non-root node + deep_node = self.sptree.search_nodes(name="A")[0] + self.assertIs(_species_name_index(deep_node), first) + + def test_mrca_via_index_matches_mrca_via_search_nodes(self): + expected = self.sptree.get_common_ancestor( + self.sptree.search_nodes(name="sp1")[0], self.sptree.search_nodes(name="sp3")[0] + ) + self.assertIs(_species_mrca(self.sptree, {"sp1", "sp3"}), expected) + + +class TaxOverlapUnitTests(TestCase): + """Unit tests for the building blocks of the Moi/Kim taxonomy-overlap TCSScore + (_species_lineage_index, _tax_overlap, _combine_tax_overlap), isolated from the rest of + attach_scores (CompletenessScore/ImpliedLosses).""" + + def setUp(self): + self.sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + for n in self.sptree.traverse(): + n.add_feature("size", len(n)) + + def test_species_lineage_is_absolute_not_relative_to_any_mrca(self): + # unlike the old mrca-relative lineage, this includes the species tree's own root (M) + # and the species' own name, regardless of what node the index happens to be requested + # from. + index = _species_lineage_index(self.sptree) + self.assertEqual(index["sp1"], frozenset({"sp1", "A", "M"})) + self.assertEqual(index["sp3"], frozenset({"sp3", "B", "M"})) + deep_node = self.sptree.search_nodes(name="A")[0] + self.assertEqual(_species_lineage_index(deep_node), index) + + def test_tax_overlap_leaf_size_comes_from_nrmembers_not_a_recount(self): + # a leaf hog whose _members set happens to hold 2 entries (e.g. a dubious merged + # fragment) should report leaf_size=2 straight from NrMembers/len(hog), not "1 per leaf + # node" the way the reference script's minimal tree class assumes. + leaf = leaf_hog("p1||sp1", self.sptree.search_nodes(name="sp1")[0]) + leaf._members = {"p1||sp1", "p1b||sp1"} + self.assertEqual(len(leaf), 2) + nset, leaf_size, leaf_acc, tax_score = _tax_overlap(leaf, _species_lineage_index(self.sptree)) + self.assertEqual(leaf_size, 2) + self.assertEqual((nset, leaf_acc, tax_score), (frozenset({"sp1", "A", "M"}), 0, 0)) + + def test_combine_ignores_an_empty_nset_instead_of_zeroing_the_whole_intersection(self): + # a part with an empty nset (an already fully-incongruent branch) must not poison an + # otherwise-matching sibling by naive-intersecting it down to nothing. + matching_a = (frozenset({"X", "Y"}), 1, 0, 0) + matching_b = (frozenset({"X", "Y", "Z"}), 1, 0, 0) + incongruent = (frozenset(), 1, 0, 0) + + nset, leaf_size, leaf_acc, tax_score = _combine_tax_overlap([matching_a, matching_b, incongruent]) + self.assertEqual(nset, frozenset({"X", "Y"})) # incongruent's empty nset was excluded + self.assertEqual(leaf_size, 3) + self.assertEqual(leaf_acc, 3) # 0+0+0 (children) + leaf_size + self.assertEqual(tax_score, len({"X", "Y"}) * 3) # 0+0+0 (children) + len(nset)*leaf_size + + class BalancedSpeciesTreeTests(TestCase): """Species tree: M(A(sp1,sp2), B(sp3,sp4))""" @@ -74,7 +206,13 @@ def test_fully_congruent_hog(self): self.assertEqual(score_value(elem, "ImpliedLosses"), "0") self.assertEqual(score_value(elem, "TCSScore"), "1.0") - def test_clean_loss_does_not_penalize_tcs(self): + def test_clean_loss_lowers_the_tcs_score(self): + # TCSScore follows the Moi/Kim taxonomy-overlap formula (see attach_scores' docstring): + # losing sp2's whole branch removes a genuine nested match, so unlike CompletenessScore + # (species-level presence) it does drop here relative to test_fully_congruent_hog's + # 1.0: hogA becomes a bare pass-through of leaf1 (nset={sp1,A,M}, tax_score=0), hogB + # unchanged (nset={B,M}, leaf_size=2, tax_score=4); folded: leaf_size=3, leaf_acc=5, + # tax_score=4+len({M})*3=7 -> (7 - 5*1) / 3 = 0.6667. leaf1 = leaf_hog("p1||sp1", self.sptree.search_nodes(name="sp1")[0]) leaf3 = leaf_hog("p3||sp3", self.sptree.search_nodes(name="sp3")[0]) leaf4 = leaf_hog("p4||sp4", self.sptree.search_nodes(name="sp4")[0]) @@ -90,9 +228,13 @@ def test_clean_loss_does_not_penalize_tcs(self): self.assertEqual(score_value(elem, "CompletenessScore"), "0.75") self.assertEqual(score_value(elem, "ImpliedLosses"), "1") - self.assertEqual(score_value(elem, "TCSScore"), "1.0") + self.assertEqual(score_value(elem, "TCSScore"), "0.6667") - def test_direct_sibling_species_have_no_tcs_score(self): + def test_direct_sibling_species_score_zero(self): + # sp1 and sp2's absolute lineages share {A, M} (nset, len 2), giving tax_score = + # len(nset)*leaf_size = 2*2 = 4 -- but leaf_acc*len(nset) = 2*2 = 4 too (a plain 2-leaf, + # single-copy family earns exactly enough "self" bonus to cancel its own match), so + # (4 - 4) / 2 = 0.0. leaf1 = leaf_hog("p1||sp1", self.sptree.search_nodes(name="sp1")[0]) leaf2 = leaf_hog("p2||sp2", self.sptree.search_nodes(name="sp2")[0]) hogA = merged_hog([leaf1, leaf2], self.A) @@ -101,7 +243,7 @@ def test_direct_sibling_species_have_no_tcs_score(self): elem = ET.Element("orthologGroup") attach_scores(elem, hogA, self.A, species_of_members) - self.assertIsNone(score_value(elem, "TCSScore")) + self.assertEqual(score_value(elem, "TCSScore"), "0.0") self.assertIsNotNone(score_value(elem, "CompletenessScore")) self.assertIsNotNone(score_value(elem, "ImpliedLosses")) @@ -121,3 +263,215 @@ def test_single_species_loss_within_a_clade(self): sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) losses = _count_implied_losses(sptree, {"sp1", "sp3", "sp4"}) self.assertEqual(losses, 1) + + +class HierarchicalImpliedLossesUnitTests(TestCase): + """_hog_implied_losses() recurses through _subhogs (unlike the flat _count_implied_losses, + which only sees a single flattened species set) so that a loss hidden behind a sibling + paralog's presence still gets counted -- while not double-counting a plain (non-duplicated) + nested loss that _count_implied_losses would already catch on its own.""" + + def setUp(self): + self.sptree = Tree("((sp1,sp2)A,(sp3,sp4)B)M;", format=1) + for n in self.sptree.traverse(): + n.add_feature("size", len(n)) + self.M, self.A, self.B = self.sptree, *(self.sptree.search_nodes(name=n)[0] for n in ("A", "B")) + + def _sp(self, name): + return self.sptree.search_nodes(name=name)[0] + + def test_matches_flat_count_without_duplication(self): + # sp2 lost under A, no duplication anywhere -- single nested loss, counted once (not + # once per ancestor level). + hogA = merged_hog([leaf_hog("p1||sp1", self._sp("sp1"))], self.A) + hogB = merged_hog([leaf_hog("p3||sp3", self._sp("sp3")), leaf_hog("p4||sp4", self._sp("sp4"))], self.B) + hogM = merged_hog([hogA, hogB], self.M) + self.assertEqual(_hog_implied_losses(hogM, self.M), 1) + + def test_duplication_loss_hidden_behind_sibling_paralog_is_still_counted(self): + # two paralogous copies under B (both tax_now=B): copy1 complete (sp3,sp4), copy2 only + # in sp3. Flattened species of B are still {sp3,sp4} (copy1 covers sp4), so the flat + # _count_implied_losses would report 0 -- but copy2 truly lost its sp4 copy. + copy1 = merged_hog([leaf_hog("p3a||sp3", self._sp("sp3")), leaf_hog("p4a||sp4", self._sp("sp4"))], self.B) + copy2 = merged_hog([leaf_hog("p3b||sp3", self._sp("sp3"))], self.B) + hogB = merged_hog([copy1, copy2], self.B) + + flat_species = _member_species(hogB.get_members()) + self.assertEqual(flat_species, {"sp3", "sp4"}) + self.assertEqual(_count_implied_losses(self.B, flat_species), 0) + + self.assertEqual(_hog_implied_losses(hogB, self.B), 1) + + +class SpeciesMrcaOnPrunedTreeTests(TestCase): + """A rootHOG never containing a given species causes prepare_species_tree() to prune that + species out of the working tree entirely. Scores must still be computed against the real, + unpruned species tree topology so such species count as losses.""" + + def setUp(self): + self.full_tree = Tree("((AQUAE,CHLTR)inter1,MYCGE)inter2;", format=1) + for n in self.full_tree.traverse(): + n.add_feature("size", len(n)) + self.full_tree_copy = self.full_tree.copy() + + # simulate prepare_species_tree(): rootHOG only ever contained CHLTR and MYCGE, so AQUAE + # (and the now-singleton inter1 node) is pruned away from the working tree. + self.pruned = self.full_tree.get_common_ancestor({"CHLTR", "MYCGE"}) + self.pruned.prune({"CHLTR", "MYCGE"}, preserve_branch_length=True) + + def test_pruned_tree_hides_the_missing_species(self): + self.assertEqual(sorted(c.name for c in self.pruned.children), ["CHLTR", "MYCGE"]) + self.assertEqual(_count_implied_losses(self.pruned, {"CHLTR", "MYCGE"}), 0) + + def test_mrca_on_the_full_tree_recovers_the_real_topology(self): + node = _species_mrca(self.full_tree_copy, {"CHLTR", "MYCGE"}) + self.assertEqual(node.name, self.pruned.name) # same node identity as the pruned mrca ... + self.assertEqual(sorted(c.name for c in node.children), ["MYCGE", "inter1"]) # ... but real topology + self.assertEqual(_count_implied_losses(node, {"CHLTR", "MYCGE"}), 1) + + def test_attach_scores_reports_the_loss_when_scored_against_the_full_tree(self): + elem = ET.Element("orthologGroup") + species_of_members = {"CHLTR", "MYCGE"} + node = _species_mrca(self.full_tree_copy, species_of_members) + hog = merged_hog( + [leaf_hog("p1||CHLTR", self.full_tree_copy.search_nodes(name="CHLTR")[0]), + leaf_hog("p2||MYCGE", self.full_tree_copy.search_nodes(name="MYCGE")[0])], + node, + ) + attach_scores(elem, hog, node, species_of_members) + + self.assertEqual(score_value(elem, "CompletenessScore"), "0.6667") + self.assertEqual(score_value(elem, "ImpliedLosses"), "1") + self.assertEqual(score_value(elem, "TCSScore"), "0.0") + + +class GeneTreeMathTCSTests(TestCase): + """Exercises attach_scores()/_tax_overlap() as pure scoring functions against + tests/data/tcs_tree*.nwk mapped directly onto species-tree MRCAs via + load_hog_from_genetree(), using the asymmetric, unbalanced species tree built from + tests/data/tcs_taxonomy.tsv (unlike the balanced synthetic tree in BalancedSpeciesTreeTests). + + NOTE: a real FastOMA-inferred HOG's sub-hog nesting always follows the species tree exactly + -- infer_hogs_for_rhog_levels_recursively() recurses over the real species tree, and + reconciliation only decides whether sibling sub-hogs at a given species-tree node get merged + (speciation) or kept apart (duplication), never reordering across species-tree clades. So + tcs_tree2/tcs_tree3 below -- whose topologies genuinely cross species-tree clades -- are not + structures the real pipeline could ever produce as a HOG's sub-hog tree; they're here purely + to check the score arithmetic in isolation. Only tcs_tree1 (fully congruent) reflects a + realistic pipeline output.""" + + def setUp(self): + self.sptree = load_species_tree_from_lineage_tsv(os.path.join(DATA_DIR, "tcs_taxonomy.tsv")) + + def _tcs_score_of(self, tree_filename): + hog = load_hog_from_genetree(os.path.join(DATA_DIR, tree_filename), self.sptree) + species_of_members = _member_species(hog.get_members()) + mrca = _species_mrca(self.sptree, species_of_members) + + elem = ET.Element("orthologGroup") + attach_scores(elem, hog, mrca, species_of_members) + return score_value(elem, "TCSScore") + + def test_gene_tree_congruent_with_species_tree_scores_highest(self): + # tcs_tree1: ((SP1,SP2),(SP3,SP4)) matches the real topology -- (SP1,SP2) under genus 4, + # (SP3,SP4) under order 7 -- exactly. The one case realistically reachable by the + # pipeline, and (see the other two tests below) the highest-scoring of the three. + self.assertEqual(self._tcs_score_of("tcs_tree1.nwk"), "2.0") + + def test_gene_tree_partially_incongruent_scores_lower(self): + # tcs_tree2: (((SP1,SP2),SP3),SP4) nests SP3 with the (SP1,SP2) clade instead of with + # SP4 as the real species tree has it, so it recovers less nested structure than + # tcs_tree1's perfect match (1.5 < 2.0) -- but more than tcs_tree3's total mismatch. + self.assertEqual(self._tcs_score_of("tcs_tree2.nwk"), "1.5") + + def test_gene_tree_fully_incongruent_scores_0(self): + # tcs_tree3: ((SP1,SP3),(SP2,SP4)) crosses both real clades ((SP1,SP2) and (SP3,SP4)), + # so no split in the gene tree agrees with the species tree at all -- every internal + # node's nset collapses to just the family's own top-level match, contributing nothing + # beyond what a single-copy family would already get for free (same cancellation as + # test_direct_sibling_species_score_zero, just one level deeper): 0.0. + self.assertEqual(self._tcs_score_of("tcs_tree3.nwk"), "0.0") + + +class DuplicationAsymmetricRetentionTCSTests(TestCase): + r"""A duplication scenario that *is* reachable by the real pipeline (unlike tcs_tree2/tcs_tree3 + above): under the order-7 clade, two paralogous copies exist -- copy1 retained cleanly in + both SP3 and SP4, copy2 retained only in SP3 (lost in SP4). + + Species tree: + 1 (class) + / \ + 2 (order) 7 (order) + | / \ + 3 (family) 8 (family) 11 (family) + | | | + 4 (genus) 9 (genus) 12 (subfamily) + / \ | | + SP1 SP2 10 (species) 13 (genus) + | | + SP3 14 (species) + | + 15 (subspecies) + | + SP4 + + HOG topology built in the test method below -- copy1 and copy2 both get _tax_now=7 (node + "7"), matching how the real pipeline's merge_subhogs() always bumps every surviving hog's + _tax_now to the level it was last processed at, whether or not it found a merge partner + there. That shared _tax_now is what marks them as a recognized duplication (paralogGroup): + + 1 (top HOG) + / \ + 4 (hogA: SP1,SP2) 7 (hog7: duplication, paralogGroup) + / \ + copy1 (tax_now=7) copy2 (tax_now=7) + / \ | + SP3 SP4 SP3 (2nd paralog, + (protA_SP3) (protA_SP4) protB_SP3 -- no SP4 + counterpart: lost) + + Every species is still present overall (some copy reaches both SP3 and SP4), so + CompletenessScore reads 1.0 -- it only sees species-level presence, not per-paralog + retention. But ImpliedLosses is computed hierarchically (_hog_implied_losses), treating + copy1 and copy2 as independent lineages so copy1's presence in SP4 can't mask copy2's loss + there: 1 implied loss (copy2 never reaching SP4). + + TCSScore (see attach_scores' docstring for the formula) works out, per _tax_overlap(), to: + hogA: nset={4,3,2,1} leaf_size=2 leaf_acc=2 tax_score=0+4*2=8 + copy1: nset={7,1} leaf_size=2 leaf_acc=2 tax_score=0+2*2=4 + copy2: nset={SP3,9,8,7,1} leaf_size=1 leaf_acc=0 tax_score=0 (bare pass-through) + hog7: fold(copy1,copy2) leaf_size=3 leaf_acc=5 tax_score=4+len({7,1})*3=10 + top: fold(hogA,hog7) leaf_size=5 leaf_acc=12 tax_score=18+len({1})*5=23 + TCSScore = (23 - 12*len({1})) / 5 = (23 - 12) / 5 = 2.2.""" + + def setUp(self): + self.sptree = load_species_tree_from_lineage_tsv(os.path.join(DATA_DIR, "tcs_taxonomy.tsv")) + + def _node(self, name): + return self.sptree.search_nodes(name=name)[0] + + def test_duplication_with_asymmetric_paralog_retention(self): + hogA = merged_hog( + [leaf_hog("p1||SP1", self._node("SP1")), leaf_hog("p2||SP2", self._node("SP2"))], + self._node("4"), + ) + copy1 = merged_hog( + [leaf_hog("p3a||SP3", self._node("SP3")), leaf_hog("p4a||SP4", self._node("SP4"))], + self._node("7"), + ) + # duplicate of copy1, lost in SP4; wrapped at tax_now=7 like copy1 (see docstring) so + # both are recognized as one duplication group rather than two unrelated lineages. + copy2 = merged_hog([leaf_hog("p3b||SP3", self._node("SP3"))], self._node("7")) + hog7 = merged_hog([copy1, copy2], self._node("7")) + top = merged_hog([hogA, hog7], self._node("1")) + + species_of_members = _member_species(top.get_members()) + self.assertEqual(species_of_members, {"SP1", "SP2", "SP3", "SP4"}) + mrca = _species_mrca(self.sptree, species_of_members) + + elem = ET.Element("orthologGroup") + attach_scores(elem, top, mrca, species_of_members) + + self.assertEqual(score_value(elem, "CompletenessScore"), "1.0") + self.assertEqual(score_value(elem, "ImpliedLosses"), "1") + self.assertEqual(score_value(elem, "TCSScore"), "2.2") From 31da9b01e85c3e8f33e55e1d9e8b1ff7dec526b1 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Sat, 29 Aug 2026 07:04:07 +0200 Subject: [PATCH 3/4] Add per-score CLI flags and only emit scoreDef entries that are present Adds --store-completeness-score/--store-implied-losses-score/--store-tcs-score options to fastoma-infer-subhogs (wired through FastOMA.nf as store_completness_score/store_implied_losses_score/store_tcs_score params) so score computation can be toggled per run. collect_subhogs now looks ahead at the first HOG in the pickle folder to only write scoreDef entries for scores actually present in the run's output, instead of always writing all three. compute time for score computation is tracked during the to_orthoxml() call. Serialization is generally quick, score computation uses roughly 40% of the time. --- FastOMA.nf | 10 +- FastOMA/_hog_class.py | 103 +++++++++++++++---- FastOMA/_infer_subhog.py | 12 ++- FastOMA/collect_subhogs.py | 31 ++++-- FastOMA/infer_subhogs.py | 6 ++ FastOMA/zoo/wrappers/treebuilders/parsers.py | 5 +- conf/test-fungi.config | 4 +- conf/test-mammalia.config | 3 +- nextflow.config | 4 + tests/test_hog_scores.py | 37 +++++++ 10 files changed, 177 insertions(+), 38 deletions(-) diff --git a/FastOMA.nf b/FastOMA.nf index 57bce58..fe101c9 100644 --- a/FastOMA.nf +++ b/FastOMA.nf @@ -222,7 +222,10 @@ process hog_big{ --gap-ratio-col ${params.filter_gap_ratio_col} \ --number-of-samples-per-hog ${params.nr_repr_per_hog} \ ${ params.write_msas ? "--msa-write" : ""} \ - ${ params.write_genetrees ? "--gene-trees-write" : ""} + ${ params.write_genetrees ? "--gene-trees-write" : ""} \ + ${ params.store_completness_score ? "" : "--no-store-completeness-score" } \ + ${ params.store_implied_losses_score ? "" : "--no-store-implied-losses-score" } \ + ${ params.store_tcs_score ? "" : "--no-store-tcs-score" } """ } @@ -255,7 +258,10 @@ process hog_rest{ --gap-ratio-col ${params.filter_gap_ratio_col} \ --number-of-samples-per-hog ${params.nr_repr_per_hog} \ ${ params.write_msas ? "--msa-write" : ""} \ - ${ params.write_genetrees ? "--gene-trees-write" : ""} + ${ params.write_genetrees ? "--gene-trees-write" : ""} \ + ${ params.store_completness_score ? "" : "--no-store-completeness-score" } \ + ${ params.store_implied_losses_score ? "" : "--no-store-implied-losses-score" } \ + ${ params.store_tcs_score ? "" : "--no-store-tcs-score" } """ } diff --git a/FastOMA/_hog_class.py b/FastOMA/_hog_class.py index 9e57ee0..9e80586 100644 --- a/FastOMA/_hog_class.py +++ b/FastOMA/_hog_class.py @@ -1,9 +1,10 @@ +import time import xml.etree.ElementTree as ET from Bio.Align import MultipleSeqAlignment from Bio.SeqRecord import SeqRecord from random import sample -from typing import Optional, List, Union +from typing import Optional, List, Union, NamedTuple from ete3 import Tree, TreeNode import random from ._utils_subhog import MSAFilter @@ -172,27 +173,69 @@ def _tax_overlap(hog: "HOG", species_lineage_index: dict): 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: +class ScoreFlags(NamedTuple): + """Which of the three orthoxml scores to compute/attach, threaded as a single immutable value + through to_orthoxml()/_to_orthoxml()/attach_scores() instead of one bool param per score.""" + store_completeness_score: bool = True + store_implied_losses_score: bool = True + store_tcs_score: bool = True + + +class _ScoreTimings: + """Dummy container threaded through one rootHOG's HOG.to_orthoxml() recursion. + + attach_scores() is called once per orthologGroup node produced during that recursion (not + just at the root), so timing a single call doesn't tell you the cost of a score across the + whole rootHOG. Passing the same instance down through every recursive to_orthoxml()/ + attach_scores() call lets each call add its own elapsed time, giving the total time spent + computing each score for the entire rootHOG once the top-level to_orthoxml() call returns.""" + + __slots__ = ("completeness_score", "implied_losses", "tcs_score") + + def __init__(self): + self.completeness_score = 0.0 + self.implied_losses = 0.0 + self.tcs_score = 0.0 + + +def attach_scores(hog_element: ET.Element, hog: "HOG", mrca: TreeNode, species_of_members: set, + timings: Optional[_ScoreTimings] = None, score_flags: ScoreFlags = ScoreFlags()) -> None: """Computes and attaches CompletenessScore, TCSScore and ImpliedLosses as - sub-elements of `hog_element`. + sub-elements of `hog_element`, each gated by its corresponding flag in `score_flags`. 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)}) + normalized by leaf_size -- the gene count of `hog` itself. + + If `timings` (a _ScoreTimings) is given, the wall time of each of the three score + computations is added to it -- see _ScoreTimings for why this needs to be an accumulator + rather than a local measurement.""" + if score_flags.store_completeness_score: + t0 = time.perf_counter() + completeness_score = round(len(species_of_members) / mrca.size, 4) + ET.SubElement(hog_element, "score", attrib={"id": "CompletenessScore", "value": str(completeness_score)}) + if timings is not None: + timings.completeness_score += time.perf_counter() - t0 + + if score_flags.store_implied_losses_score: + t0 = time.perf_counter() + 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)}) + if timings is not None: + timings.implied_losses += time.perf_counter() - t0 - 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))}) + if score_flags.store_tcs_score: + t0 = time.perf_counter() + 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))}) + if timings is not None: + timings.tcs_score += time.perf_counter() - t0 # from .infer_subhogs import conf_infer_subhhogs #fastoma_infer_subhogs # @@ -429,7 +472,27 @@ def merge_prots_msa(self, merged_fragment_name, merged_msa_new): # merged_frag # self._msa = MultipleSeqAlignment(msa_new) # return 1 - def to_orthoxml(self, full_species_tree: Optional[TreeNode] = None): + def to_orthoxml(self, full_species_tree: Optional[TreeNode] = None, score_flags: ScoreFlags = ScoreFlags()): + """Public entry point: builds the orthoxml element tree for this HOG and, once the whole + (recursive) build is done, logs how much of the total time went into computing scores + (see _ScoreTimings and attach_scores). `score_flags` is forwarded to every attach_scores() + call made during the recursion.""" + timings = _ScoreTimings() + start = time.perf_counter() + result = self._to_orthoxml(full_species_tree, timings, score_flags) + elapsed = time.perf_counter() - start + if elapsed > 0: + score_total = timings.completeness_score + timings.implied_losses + timings.tcs_score + logger.info( + "to_orthoxml for rootHOG %s took %.3fs; scoring took %.3fs (%.1f%%) of which " + "CompletenessScore=%.3fs ImpliedLosses=%.3fs TCSScore=%.3fs", + self._rhogid, elapsed, score_total, 100 * score_total / elapsed, + timings.completeness_score, timings.implied_losses, timings.tcs_score, + ) + return result + + def _to_orthoxml(self, full_species_tree: Optional[TreeNode], timings: _ScoreTimings, + score_flags: ScoreFlags = ScoreFlags()): if len(self._subhogs) == 0: list_member = list(self._members) if len(list_member) == 1: @@ -485,7 +548,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(full_species_tree) + element_p = sh._to_orthoxml(full_species_tree, timings, score_flags) if str(element_p): paralog_element.append(element_p) # ,**gene_id_name indent+2 else: @@ -495,7 +558,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(full_species_tree) + element = subhog._to_orthoxml(full_species_tree, timings, score_flags) if str(element): # element could be element_list.append(element) # indent+2 else: @@ -515,8 +578,8 @@ def _sorter_key(sh): logger.info(f"mrca ({mrca.name}) != self.taxlevel ({self.taxlevel.name})") logger.info(f"<{hog_elemnt.tag} {hog_elemnt.attrib}>") - attach_scores(hog_elemnt, self, mrca, species_of_members) - property_element = ET.SubElement(hog_elemnt, "property", attrib={"name": "TaxRange", "value": str(mrca.name)}) + attach_scores(hog_elemnt, self, mrca, species_of_members, timings=timings, score_flags=score_flags) + ET.SubElement(hog_elemnt, "property", attrib={"name": "TaxRange", "value": str(mrca.name)}) for element in element_list: hog_elemnt.append(element) diff --git a/FastOMA/_infer_subhog.py b/FastOMA/_infer_subhog.py index 93a4455..ad14bd6 100644 --- a/FastOMA/_infer_subhog.py +++ b/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, attach_scores, _member_species, _species_name_index +from ._hog_class import HOG, Representative, split_hog, attach_scores, _member_species, _species_name_index, ScoreFlags from ._utils_subhog import MSAFilter, MSAFilterElbow, MSAFilterTrimAL from .zoo.utils import unique @@ -96,13 +96,19 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi if not keep_subhog_each_pickle: shutil.rmtree(pickles_subhog_folder) + score_flags = ScoreFlags( + store_completeness_score=conf_infer_subhhogs.store_completeness_score, + store_implied_losses_score=conf_infer_subhhogs.store_implied_losses_score, + store_tcs_score=conf_infer_subhhogs.store_tcs_score, + ) + tot_genes, placed_genes = 0, 0 hogs_rhogs_xml = [] for hog_i in hogs_a_rhog: 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(full_species_tree) # + hogs_a_rhog_xml_raw = hog_i.to_orthoxml(full_species_tree, score_flags) # 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)}) @@ -110,7 +116,7 @@ def read_infer_xml_rhog(rhogid, inferhog_concurrent_on, pickles_rhog_folder, pi 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) + attach_scores(hog_elemnt, hog_i, scoring_node, species_of_members, score_flags=score_flags) 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 diff --git a/FastOMA/collect_subhogs.py b/FastOMA/collect_subhogs.py index e21a6a3..3a2d0ee 100644 --- a/FastOMA/collect_subhogs.py +++ b/FastOMA/collect_subhogs.py @@ -26,6 +26,24 @@ # This code collect subhogs and writes outputs. +SCORE_DEFS = { + "CompletenessScore": "Fraction of expected species with genes in the (Sub)HOG", + "TCSScore": "Taxonomic Congruence Score: how well the (Sub)HOG structure matches the species tree topology", + "ImpliedLosses": "Number of implied gene loss events (Dollo parsimony) within the (Sub)HOG's taxonomic range", +} + + +def peek_score_ids(pickle_folder: Path) -> set: + """Returns the set of score ids attached to the first HOG found in pickle_folder. + + Score ids depend only on the (uniform, run-wide) --store-*-score flags passed to + fastoma-infer-subhogs, not on individual HOGs, so a single HOG is representative of + the whole run and stops the scan after the first non-empty pickle file.""" + for hog in iter_hogs(pickle_folder): + return {child.get('id') for child in hog if child.tag == 'score'} + return set() + + def iter_hogs(pickle_folder: Path): cnt = 0 nr_hogs = 0 @@ -193,13 +211,12 @@ def write_hog_orthoxml(pickle_folder, output_xml_name, gene_id_pickle_file, id_t logger.debug("gene_xml is created.") orthoxml_file.append(taxonomy) - 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"}) + present_score_ids = peek_score_ids(Path(pickle_folder)) + if present_score_ids: + scores = ET.SubElement(orthoxml_file, "scores") + for score_id, desc in SCORE_DEFS.items(): + if score_id in present_score_ids: + ET.SubElement(scores, "scoreDef", {"id": score_id, "desc": desc}) # #### create the groups of orthoxml #### groups_xml = ET.SubElement(orthoxml_file, "groups") diff --git a/FastOMA/infer_subhogs.py b/FastOMA/infer_subhogs.py index 62d6379..0445330 100644 --- a/FastOMA/infer_subhogs.py +++ b/FastOMA/infer_subhogs.py @@ -49,6 +49,12 @@ def fastoma_infer_subhogs(): help="For trimming the MSA, the threshold of ratio of gaps for each column.") parser.add_argument("--min-col-trim", required=False, type=int, default=50, # todo min rows trim help="min no. columns in msa to consider for filtering") + parser.add_argument("--store-completeness-score", action=argparse.BooleanOptionalAction, default=True, + help="Store the CompletenessScore of each (Sub)HOG in the output orthoxml.") + parser.add_argument("--store-implied-losses-score", action=argparse.BooleanOptionalAction, default=True, + help="Store the ImpliedLosses score of each (Sub)HOG in the output orthoxml.") + parser.add_argument("--store-tcs-score", action=argparse.BooleanOptionalAction, default=True, + help="Store the TCSScore of each (Sub)HOG in the output orthoxml.") parser.add_argument('-v', action="count", default=0, help="Increase verbosity to info/debug") conf_infer_subhhogs = parser.parse_args() diff --git a/FastOMA/zoo/wrappers/treebuilders/parsers.py b/FastOMA/zoo/wrappers/treebuilders/parsers.py index 86e69e4..1cdff11 100644 --- a/FastOMA/zoo/wrappers/treebuilders/parsers.py +++ b/FastOMA/zoo/wrappers/treebuilders/parsers.py @@ -7,9 +7,8 @@ logger.addHandler(logging.StreamHandler()) - -FLOAT = Word(nums + '.-').setParseAction(lambda x: float(x[0])) -INT = Word(nums).setParseAction(lambda x: int(x[0])) +FLOAT = Word(nums + '.-').set_parse_action(lambda x: float(x[0])) +INT = Word(nums).set_parse_action(lambda x: int(x[0])) WORD = Word(alphanums + '_') SPACEDWORD = Word(alphanums + ' _') MODEL_CONCAT_WORD = Word(alphanums + "+_-") diff --git a/conf/test-fungi.config b/conf/test-fungi.config index 7025ab8..e257262 100644 --- a/conf/test-fungi.config +++ b/conf/test-fungi.config @@ -1,7 +1,7 @@ // Default configuration for Nextflow params { - test_data_url = "https://zenodo.org/records/17434495/files/fungi-30.tgz?download=1" + input = "https://zenodo.org/records/17434495/files/fungi-30.tgz?download=1" report = true - omamer_db = "${projectDir}/testdata/test.h5" + omamer_db = "https://zenodo.org/records/20814376/files/saccharomyceta.h5?download=1" } diff --git a/conf/test-mammalia.config b/conf/test-mammalia.config index b29ca3d..f6a4302 100644 --- a/conf/test-mammalia.config +++ b/conf/test-mammalia.config @@ -1,6 +1,7 @@ // Default configuration for Nextflow params { - test_data_url = "https://zenodo.org/records/17434495/files/mammalia-22.tgz?download=1" + input = "https://zenodo.org/records/17434495/files/mammalia-22.tgz?download=1" + omamer_db = "https://zenodo.org/records/20814376/files/Metazoa.h5?download=1" report = true } diff --git a/nextflow.config b/nextflow.config index e88ebb4..b06cc96 100644 --- a/nextflow.config +++ b/nextflow.config @@ -63,6 +63,10 @@ params { nr_repr_per_hog = 5 min_sequence_length = 20 + // scores written to FastOMA_HOGs.orthoxml + store_completness_score = true + store_implied_losses_score = true + store_tcs_score = true // other parameters debug_enabled = false // generates additional files for debugging diff --git a/tests/test_hog_scores.py b/tests/test_hog_scores.py index 9a69549..2d063a6 100644 --- a/tests/test_hog_scores.py +++ b/tests/test_hog_scores.py @@ -6,6 +6,7 @@ from FastOMA._hog_class import ( HOG, + ScoreFlags, attach_scores, _species_mrca, _species_name_index, @@ -247,6 +248,42 @@ def test_direct_sibling_species_score_zero(self): self.assertIsNotNone(score_value(elem, "CompletenessScore")) self.assertIsNotNone(score_value(elem, "ImpliedLosses")) + def test_score_flags_disable_completeness_score(self): + hogM = self._make_full_hog() + species_of_members = _member_species(hogM.get_members()) + + elem = ET.Element("orthologGroup") + attach_scores(elem, hogM, self.M, species_of_members, + score_flags=ScoreFlags(store_completeness_score=False)) + + self.assertIsNone(score_value(elem, "CompletenessScore")) + self.assertIsNotNone(score_value(elem, "ImpliedLosses")) + self.assertIsNotNone(score_value(elem, "TCSScore")) + + def test_score_flags_disable_implied_losses_score(self): + hogM = self._make_full_hog() + species_of_members = _member_species(hogM.get_members()) + + elem = ET.Element("orthologGroup") + attach_scores(elem, hogM, self.M, species_of_members, + score_flags=ScoreFlags(store_implied_losses_score=False)) + + self.assertIsNotNone(score_value(elem, "CompletenessScore")) + self.assertIsNone(score_value(elem, "ImpliedLosses")) + self.assertIsNotNone(score_value(elem, "TCSScore")) + + def test_score_flags_disable_tcs_score(self): + hogM = self._make_full_hog() + species_of_members = _member_species(hogM.get_members()) + + elem = ET.Element("orthologGroup") + attach_scores(elem, hogM, self.M, species_of_members, + score_flags=ScoreFlags(store_tcs_score=False)) + + self.assertIsNotNone(score_value(elem, "CompletenessScore")) + self.assertIsNotNone(score_value(elem, "ImpliedLosses")) + self.assertIsNone(score_value(elem, "TCSScore")) + class ImpliedLossesUnitTests(TestCase): def test_whole_missing_clade_counts_as_one_loss(self): From 7d304416e609d2311eab6227e52a26ec09685638 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Sat, 29 Aug 2026 08:14:51 +0200 Subject: [PATCH 4/4] Rename score store_* params to disable_* and add to nextflow_schema.json The previous --store_*_score params (default true) were missing from nextflow_schema.json, so nf-schema parsed CLI-supplied values as plain strings instead of booleans -- a non-empty string is truthy in Groovy, so --store_tcs_score false never actually disabled anything. Switch to --disable_*_score (default false): setting the flag means the user actively wants to disable the score, so bare presence of the flag conveys the intent without needing a boolean value at all. --- FastOMA.nf | 12 ++++++------ nextflow.config | 8 ++++---- nextflow_schema.json | 12 ++++++++++++ 3 files changed, 22 insertions(+), 10 deletions(-) diff --git a/FastOMA.nf b/FastOMA.nf index fe101c9..e3a27d2 100644 --- a/FastOMA.nf +++ b/FastOMA.nf @@ -223,9 +223,9 @@ process hog_big{ --number-of-samples-per-hog ${params.nr_repr_per_hog} \ ${ params.write_msas ? "--msa-write" : ""} \ ${ params.write_genetrees ? "--gene-trees-write" : ""} \ - ${ params.store_completness_score ? "" : "--no-store-completeness-score" } \ - ${ params.store_implied_losses_score ? "" : "--no-store-implied-losses-score" } \ - ${ params.store_tcs_score ? "" : "--no-store-tcs-score" } + ${ params.disable_completeness_score.toString().toBoolean() ? "--no-store-completeness-score" : "" } \ + ${ params.disable_implied_losses_score.toString().toBoolean() ? "--no-store-implied-losses-score" : "" } \ + ${ params.disable_tcs_score.toString().toBoolean() ? "--no-store-tcs-score" : "" } """ } @@ -259,9 +259,9 @@ process hog_rest{ --number-of-samples-per-hog ${params.nr_repr_per_hog} \ ${ params.write_msas ? "--msa-write" : ""} \ ${ params.write_genetrees ? "--gene-trees-write" : ""} \ - ${ params.store_completness_score ? "" : "--no-store-completeness-score" } \ - ${ params.store_implied_losses_score ? "" : "--no-store-implied-losses-score" } \ - ${ params.store_tcs_score ? "" : "--no-store-tcs-score" } + ${ params.disable_completeness_score.toString().toBoolean() ? "--no-store-completeness-score" : "" } \ + ${ params.disable_implied_losses_score.toString().toBoolean() ? "--no-store-implied-losses-score" : "" } \ + ${ params.disable_tcs_score.toString().toBoolean() ? "--no-store-tcs-score" : "" } """ } diff --git a/nextflow.config b/nextflow.config index b06cc96..0efae03 100644 --- a/nextflow.config +++ b/nextflow.config @@ -63,10 +63,10 @@ params { nr_repr_per_hog = 5 min_sequence_length = 20 - // scores written to FastOMA_HOGs.orthoxml - store_completness_score = true - store_implied_losses_score = true - store_tcs_score = true + // scores written to FastOMA_HOGs.orthoxml -- set to disable a given score + disable_completeness_score = false + disable_implied_losses_score = false + disable_tcs_score = false // other parameters debug_enabled = false // generates additional files for debugging diff --git a/nextflow_schema.json b/nextflow_schema.json index 7e4c827..7be0e30 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -171,6 +171,18 @@ "type": "boolean", "description": "Force generation of pairwise orthologs even for large datasets (may be slow)", "help": "By default, pairwise orthologs are not generated for datasets with more than 25 species to avoid long runtimes. Enable this option to override this behavior." + }, + "disable_completeness_score": { + "type": "boolean", + "description": "Disable computing/storing the CompletenessScore of each (Sub)HOG in the output orthoxml" + }, + "disable_implied_losses_score": { + "type": "boolean", + "description": "Disable computing/storing the ImpliedLosses score of each (Sub)HOG in the output orthoxml" + }, + "disable_tcs_score": { + "type": "boolean", + "description": "Disable computing/storing the TCSScore of each (Sub)HOG in the output orthoxml" } } },