From eb40beed7d8fbcf83e9a983396671524796bfd49 Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Tue, 18 Aug 2026 13:23:45 +0200 Subject: [PATCH 1/2] [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/2] 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")