analysis

Extraction of evolutionary information from simulated scenarios.

best_matches

Orthology graph, BMG, RBMG.

This module provides classes concerning colored graphs and colored (phylogenetic) trees, including a BMG(Tree)-function, and computation of anLRT from a (not necessarily valid) BMG or from a (pruned) gene tree.

Phylogenetic best matches of a gene x of species X are defined as those genes y of another species Y that share the lowest common ancestor with x in the gene tree among all genes in that species. In contrast, two genes are orthologs if their last common ancestor was a speciation event. Orthology and reciprocal best matches are closely related. The module asymmetree.analysis. BestMatches contains methods for extracting BMGs from trees, least resolved trees of BMG, and orthology graphs from trees.

TwoColoredLRT

Construct the least-resolved tree of a 2-colored best match graph.

Implementation of the efficient algorithm in [1].

References

1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Least resolved trees for two-colored best match graphs. In: Journal of Graph Algorithms and Applications, 25(1):397-416, 2021. doi: 10.7155/jgaa.00564.

Source code in asymmetree/analysis/best_matches.py

class TwoColoredLRT:
    """Construct the least-resolved tree of a 2-colored best match graph.

    Implementation of the efficient algorithm in [1].

    References:
        1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Least resolved trees for
           two-colored best match graphs. In: Journal of Graph Algorithms and Applications,
           25(1):397-416, 2021. doi: 10.7155/jgaa.00564.
    """

    def __init__(self, digraph: nx.DiGraph) -> None:
        """Construct the least-resolved tree of a 2-colored best match graph.

        Args:
            digraph: A digraph whose nodes have the 'color' attribute.
        """
        if not isinstance(digraph, nx.DiGraph):
            raise TypeError("input must be a digraph")

        self.digraph = digraph
        self.color_dict = sort_by_colors(digraph)

    def build_tree(self) -> Tree | None:
        """Construct the least-resolved tree.

        Returns:
            The least-resolved tree if the graph was a 2-colored BMG, and None otherwise.
        """
        if not is_properly_colored(self.digraph):
            raise RuntimeError("not a properly colored digraph")
        if len(self.color_dict) > 2:
            raise RuntimeError("more than 2 colors")

        # star tree if there is only one color
        if len(self.color_dict) == 1:
            root = TreeNode()
            for v in self.digraph.nodes():
                root.add_child(TreeNode(label=v))
        # 2 colors
        else:
            subtrees = []
            for wcc in nx.weakly_connected_components(self.digraph):
                if len(wcc) == 1:
                    return

                subroot = self._build_tree(self.digraph.subgraph(wcc).copy())

                if not subroot:
                    return
                else:
                    subtrees.append(subroot)

            if len(subtrees) == 1:
                root = subtrees[0]
            else:
                root = TreeNode()
                for subroot in subtrees:
                    root.add_child(subroot)

        return _finalize(root, self.digraph)

    def _build_tree(self, G: nx.DiGraph) -> TreeNode | None:
        """Recursive helper method for build_tree.

        Args:
            G: A digraph whose nodes have the 'color' attribute.

        Returns:
            A TreeNode representing the least-resolved tree of G if G is a 2-colored BMG, and None
                otherwise.
        """
        color_count = {c: 0 for c in self.color_dict.keys()}
        for v in G.nodes():
            color_count[G.nodes[v]["color"]] += 1
        other_color = {color: G.order() - count for color, count in color_count.items()}

        umbrella = {v for v in G.nodes() if other_color[G.nodes[v]["color"]] == G.out_degree(v)}
        S_1 = {v for v in umbrella if umbrella.issuperset(G.predecessors(v))}
        S_2 = {v for v in S_1 if S_1.issuperset(G.predecessors(v))}

        if not S_2 or len(S_1) != len(S_2):
            return

        node = TreeNode()
        for v in S_2:
            node.add_child(TreeNode(label=v))
            G.remove_node(v)

        for wcc in nx.weakly_connected_components(G):
            if len(wcc) == 1:
                return

            child = self._build_tree(G.subgraph(wcc).copy())

            if not child:
                return
            else:
                node.add_child(child)

        return node

__init__(digraph: nx.DiGraph) -> None

Construct the least-resolved tree of a 2-colored best match graph.

Parameters:

Name	Type	Description	Default
digraph	DiGraph	A digraph whose nodes have the 'color' attribute.	required

Source code in asymmetree/analysis/best_matches.py

def __init__(self, digraph: nx.DiGraph) -> None:
    """Construct the least-resolved tree of a 2-colored best match graph.

    Args:
        digraph: A digraph whose nodes have the 'color' attribute.
    """
    if not isinstance(digraph, nx.DiGraph):
        raise TypeError("input must be a digraph")

    self.digraph = digraph
    self.color_dict = sort_by_colors(digraph)

build_tree() -> Tree | None

Construct the least-resolved tree.

Returns:

Type	Description
Tree | None	The least-resolved tree if the graph was a 2-colored BMG, and None otherwise.

Source code in asymmetree/analysis/best_matches.py

def build_tree(self) -> Tree | None:
    """Construct the least-resolved tree.

    Returns:
        The least-resolved tree if the graph was a 2-colored BMG, and None otherwise.
    """
    if not is_properly_colored(self.digraph):
        raise RuntimeError("not a properly colored digraph")
    if len(self.color_dict) > 2:
        raise RuntimeError("more than 2 colors")

    # star tree if there is only one color
    if len(self.color_dict) == 1:
        root = TreeNode()
        for v in self.digraph.nodes():
            root.add_child(TreeNode(label=v))
    # 2 colors
    else:
        subtrees = []
        for wcc in nx.weakly_connected_components(self.digraph):
            if len(wcc) == 1:
                return

            subroot = self._build_tree(self.digraph.subgraph(wcc).copy())

            if not subroot:
                return
            else:
                subtrees.append(subroot)

        if len(subtrees) == 1:
            root = subtrees[0]
        else:
            root = TreeNode()
            for subroot in subtrees:
                root.add_child(subroot)

    return _finalize(root, self.digraph)

orthology_from_tree(tree: Tree) -> nx.Graph

Constructs the true orthology graph from an event-labeled tree.

Parameters:

Name	Type	Description	Default
tree	Tree	A tree with leaves that have the label, and inner vertices that have the event attribute set.	required

Returns:

Type	Description
Graph	The orthology graph explained by the tree, i.e., a graph on the set of leaf labels and edges between two of them if and only if their last common ancestor was a speciation event.

Source code in asymmetree/analysis/best_matches.py

def orthology_from_tree(tree: Tree) -> nx.Graph:
    """Constructs the true orthology graph from an event-labeled tree.

    Args:
        tree: A tree with leaves that have the `label`, and inner vertices that have the `event`
            attribute set.

    Returns:
        The orthology graph explained by the tree, i.e., a graph on the set of leaf labels and edges
            between two of them if and only if their last common ancestor was a speciation event.
    """
    if not isinstance(tree, Tree):
        raise TypeError("not of type 'Tree'")

    leaves = tree.leaf_dict()
    TOG = nx.Graph()

    for v in leaves[tree.root]:
        TOG.add_node(v.label, color=v.reconc)

    for node in tree.preorder():
        if node.event != "S":
            continue

        for c1, c2 in itertools.combinations(node.children, 2):
            for u in leaves[c1]:
                for v in leaves[c2]:
                    TOG.add_edge(u.label, v.label)

    return TOG

bmg_from_tree(tree: Tree, supply_rbmg: bool = False) -> nx.DiGraph | tuple[nx.DiGraph, nx.Graph]

Construct a BMG (and optionally RBMG) from a given tree.

Parameters:

Name	Type	Description	Default
tree	Tree	A tree with leaves that have the label and reconc attribute set.	required
supply_rbmg	bool	If True, also return the symmetric part of the constructed BMG.	False

Returns:

Type	Description
DiGraph | tuple[DiGraph, Graph]	The constructed BMG and optionally its symmetric part.

References

1. M. Geiß, E. Chávez, M. González Laffitte, A. López Sánchez, B. M. R. Stadler, D. I. Valdivia, M. Hellmuth, M. Hernández Rosales, and P. F.Stadler. Best match graphs. In: Journal of Mathematical Biology, 78:2015-2057, 2019. doi:10.1007/s00285-019-01332-9.

Source code in asymmetree/analysis/best_matches.py

def bmg_from_tree(
    tree: Tree,
    supply_rbmg: bool = False,
) -> nx.DiGraph | tuple[nx.DiGraph, nx.Graph]:
    """Construct a BMG (and optionally RBMG) from a given tree.

    Args:
        tree: A tree with leaves that have the `label` and `reconc` attribute set.
        supply_rbmg: If True, also return the symmetric part of the constructed BMG.

    Returns:
        The constructed BMG and optionally its symmetric part.

    References:
        1. M. Geiß, E. Chávez, M. González Laffitte, A. López Sánchez, B. M. R. Stadler, D. I.
           Valdivia, M. Hellmuth, M. Hernández Rosales, and P. F.Stadler. Best match graphs.
           In: Journal of Mathematical Biology, 78:2015-2057, 2019. doi:10.1007/s00285-019-01332-9.
    """
    if not isinstance(tree, Tree):
        raise TypeError("not of type 'Tree'")

    leaves = tree.leaf_dict()
    bmg = nx.DiGraph()
    colors = set()

    for v in leaves[tree.root]:
        colors.add(v.reconc)
        bmg.add_node(v.label, color=v.reconc)

    for u in leaves[tree.root]:
        remaining = colors - set([u.reconc])  # colors to which no best match has yet been found
        parent = u.parent  # start with direct parent of each node

        while remaining and parent:
            colors_here = set()
            for v in leaves[parent]:
                if v.reconc in remaining:  # best match found
                    colors_here.add(v.reconc)
                    bmg.add_edge(u.label, v.label)
            remaining -= colors_here
            parent = parent.parent

    if not supply_rbmg:
        return bmg
    else:
        return bmg, bmg.to_undirected(reciprocal=True)

extended_best_hits(leaves: list[TreeNode], D: np.ndarray, epsilon: float = 1e-08, supply_rbmg: bool = False) -> nx.DiGraph | tuple[nx.DiGraph, nx.Graph]

Compute BMG and RBMG from a distances matrix D.

Apporximation of best matches as the genes within a tolerance range from the gene with the smallest distance, see [1].

Parameters:

Name	Type	Description	Default
leaves	list[TreeNode]	List of TreeNode instances with label attribute.	required
D	ndarray	The distances of the leaves where the order of the rows/columns corresponds to that in the list.	required
epsilon	float	Relative threshold for the inclusion of suboptimal hits, i.e., (x,y) in BMG if D(x,y) <= (1 + epsilon) * min d(x,y').	1e-08
supply_rbmg	bool	If True, also return the symmetric part of the constructed BMG.	False

Returns:

Type	Description
DiGraph | tuple[DiGraph, Graph]	The constructed BMG and optionally its symmetric part.

References

1. P. F. Stadler, M. Geiß, D. Schaller, A. López Sánchez, M. Gonzalez Laffitte, D. I. Valdivia, M. Hellmuth, and M. Hernández Rosales. From pairs of most similar sequences to phylogenetic best matches. In: Algorithms for Molecular Biology, 15:5, 2020. doi: 10.1186/s13015-020-00165-2.

Source code in asymmetree/analysis/best_matches.py

def extended_best_hits(
    leaves: list[TreeNode],
    D: np.ndarray,
    epsilon: float = 1e-8,
    supply_rbmg: bool = False,
) -> nx.DiGraph | tuple[nx.DiGraph, nx.Graph]:
    """Compute BMG and RBMG from a distances matrix D.

    Apporximation of best matches as the genes within a tolerance range from the gene with the
    smallest distance, see [1].

    Args:
        leaves: List of TreeNode instances with `label` attribute.
        D: The distances of the leaves where the order of the rows/columns corresponds to that in
            the list.
        epsilon: Relative threshold for the inclusion of suboptimal hits, i.e., (x,y) in BMG if
            D(x,y) <= (1 + epsilon) * min d(x,y').
        supply_rbmg: If True, also return the symmetric part of the constructed BMG.

    Returns:
        The constructed BMG and optionally its symmetric part.

    References:
        1. P. F. Stadler, M. Geiß, D. Schaller, A. López Sánchez, M. Gonzalez Laffitte, D. I.
           Valdivia, M. Hellmuth, and M. Hernández Rosales. From pairs of most similar sequences to
           phylogenetic best matches. In: Algorithms for Molecular Biology, 15:5, 2020.
           doi: 10.1186/s13015-020-00165-2.
    """
    bmg = nx.DiGraph()
    colors = set()
    relative_threshold = 1 + epsilon

    for v in leaves:
        bmg.add_node(v.label, color=v.reconc)
        colors.add(v.reconc)

    # build BMG
    for i, u in enumerate(leaves):
        minima = {c: float("inf") for c in colors}
        for j, v in enumerate(leaves):
            if D[i, j] < minima[v.reconc]:
                minima[v.reconc] = D[i, j]

        for j, v in enumerate(leaves):
            if u.reconc != v.reconc and D[i, j] <= relative_threshold * minima[v.reconc]:
                bmg.add_edge(u.label, v.label, distance=D[i, j])

    if not supply_rbmg:
        return bmg
    else:
        return bmg, bmg.to_undirected(reciprocal=True)

informative_triples(graph: nx.DiGraph, color_dict: dict | None = None) -> list[tuple]

Compute the informative triples of a colored digraph.

Parameters:

Name	Type	Description	Default
graph	DiGraph	A digraph whose nodes have the 'color' attribute.	required
color_dict	dict | None	A dict containing the colors that appear in the graph as keys and the list of nodes with the respective color as values. If not provided, this dict is computed internally.	None

Returns:

Type	Description
list[tuple]	A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is an informative triple in the graph, see [1].

References

1. M. Geiß, E. Chávez, M. González Laffitte, A. López Sánchez, B. M. R. Stadler, D. I. Valdivia, M. Hellmuth, M. Hernández Rosales, and P. F. Stadler. Best match graphs. In: Journal of Mathematical Biology, 78:2015-2057, 2019. doi:10.1007/s00285-019-01332-9.

Source code in asymmetree/analysis/best_matches.py

def informative_triples(
    graph: nx.DiGraph,
    color_dict: dict | None = None,
) -> list[tuple]:
    """Compute the informative triples of a colored digraph.

    Args:
        graph: A digraph whose nodes have the 'color' attribute.
        color_dict: A dict containing the colors that appear in the graph as keys and the list of
            nodes with the respective color as values. If not provided, this dict is computed
            internally.

    Returns:
        A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is an
            informative triple in the graph, see [1].

    References:
        1. M. Geiß, E. Chávez, M. González Laffitte, A. López Sánchez, B. M. R. Stadler, D. I.
           Valdivia, M. Hellmuth, M. Hernández Rosales, and P. F. Stadler. Best match graphs.
           In: Journal of Mathematical Biology, 78:2015-2057, 2019. doi:10.1007/s00285-019-01332-9.
    """
    if not isinstance(graph, nx.DiGraph):
        raise TypeError("must be a NetworkX 'Digraph'")

    if not color_dict:
        color_dict = sort_by_colors(graph)

    R = []

    for c1, c2 in itertools.permutations(color_dict.keys(), 2):
        for a in color_dict[c1]:
            for b1, b2 in itertools.permutations(color_dict[c2], 2):
                if graph.has_edge(a, b1) and (not graph.has_edge(a, b2)):
                    R.append((a, b1, b2))

    return R

forbidden_triples(graph: nx.DiGraph, color_dict: dict | None = None) -> list[tuple]

Compute the forbidden triples of a colored digraph.

Parameters:

Name	Type	Description	Default
graph	DiGraph	A digraph whose nodes have the 'color' attribute.	required
color_dict	dict | None	A dict containing the colors that appear in the graph as keys and the list of nodes with the respective color as values. If not provided, this dict is computed internally.	None

Returns:

Type	Description
list[tuple]	A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is a forbidden triple in the graph, see [1].

References

1. D. Schaller, P. F. Stadler, and M. Hellmuth. Complexity of modification problems for best match graphs. In: Theoretical Computer Science, 865:63-84, 2021. doi:10.1016/j.tcs.2021.02.037.

Source code in asymmetree/analysis/best_matches.py

def forbidden_triples(graph: nx.DiGraph, color_dict: dict | None = None) -> list[tuple]:
    """Compute the forbidden triples of a colored digraph.

    Args:
        graph: A digraph whose nodes have the 'color' attribute.
        color_dict: A dict containing the colors that appear in the graph as keys and the list of
            nodes with the respective color as values. If not provided, this dict is computed
            internally.

    Returns:
        A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is a
            forbidden triple in the graph, see [1].

    References:
        1. D. Schaller, P. F. Stadler, and M. Hellmuth. Complexity of modification problems for
           best match graphs. In: Theoretical Computer Science, 865:63-84, 2021.
           doi:10.1016/j.tcs.2021.02.037.
    """
    if not isinstance(graph, nx.DiGraph):
        raise TypeError("must be a NetworkX 'Digraph'")

    if not color_dict:
        color_dict = sort_by_colors(graph)

    F = []

    for c1, c2 in itertools.permutations(color_dict.keys(), 2):
        for a in color_dict[c1]:
            for b1, b2 in itertools.combinations(color_dict[c2], 2):
                if graph.has_edge(a, b1) and graph.has_edge(a, b2):
                    F.append((a, b1, b2))
                    F.append((a, b2, b1))

    return F

informative_forbidden_triples(graph: nx.DiGraph, color_dict: dict | None = None) -> tuple[list[tuple], list[tuple]]

Compute the informative and forbidden triples of a colored digraph.

Parameters:

Name	Type	Description	Default
graph	DiGraph	A digraph whose nodes have the 'color' attribute.	required
color_dict	dict | None	A dict containing the colors that appear in the graph as keys and the list of nodes with the respective color as values. If not provided, this dict is computed internally.	None

Returns:

Type	Description
list[tuple]	A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is an informative triple in the graph, see [1].
list[tuple]	A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is a forbidden triple in the graph, see [1].

References

1. D. Schaller, P. F. Stadler, and M. Hellmuth. Complexity of modification problems for best match graphs. In: Theoretical Computer Science, 865:63-84, 2021. doi:10.1016/j.tcs.2021.02.037.

Source code in asymmetree/analysis/best_matches.py

def informative_forbidden_triples(
    graph: nx.DiGraph, color_dict: dict | None = None
) -> tuple[list[tuple], list[tuple]]:
    """Compute the informative and forbidden triples of a colored digraph.

    Args:
        graph: A digraph whose nodes have the 'color' attribute.
        color_dict: A dict containing the colors that appear in the graph as keys and the list of
            nodes with the respective color as values. If not provided, this dict is computed
            internally.

    Returns:
        A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is an
            informative triple in the graph, see [1].
        A list where each item is a tuple (a, b1, b2) of three nodes such that a b1 | b2 is a
            forbidden triple in the graph, see [1].

    References:
        1. D. Schaller, P. F. Stadler, and M. Hellmuth. Complexity of modification problems for
           best match graphs. In: Theoretical Computer Science, 865:63-84, 2021.
           doi:10.1016/j.tcs.2021.02.037.
    """
    if not isinstance(graph, nx.DiGraph):
        raise TypeError("must be a NetworkX 'Digraph'")

    if not color_dict:
        color_dict = sort_by_colors(graph)

    R, F = [], []

    for c1, c2 in itertools.permutations(color_dict.keys(), 2):
        for a in color_dict[c1]:
            for b1, b2 in itertools.combinations(color_dict[c2], 2):
                if graph.has_edge(a, b1) and graph.has_edge(a, b2):
                    F.append((a, b1, b2))
                    F.append((a, b2, b1))
                elif graph.has_edge(a, b1):
                    R.append((a, b1, b2))
                elif graph.has_edge(a, b2):
                    R.append((a, b2, b1))

    return R, F

binary_explainable_triples(graph: nx.DiGraph, color_dict: dict | None = None) -> list[tuple]

Extended informative triple set for binary-explainable graphs.

Parameters:

Name	Type	Description	Default
graph	DiGraph	A digraph whose nodes have the 'color' attribute.	required
color_dict	dict | None	A dict containing the colors that appear in the graph as keys and the list of nodes with the respective color as values. If not provided, this dict is computed internally.	None

Returns:

Type	Description
list[tuple]	List of tuples where each item is a tuple (a, b, c) of three nodes such that ab|c is a triple displayed by every binary tree that explains the graph (if existent), see [1].

References

1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Best match graphs with binary trees. In: C. Martín-Vide, M. A. Vega-Rodríguez, and T. Wheeler, editors, Algorithms for Computational Biology, 8th AlCoB, volume 12715 of Lecture Notes in Computer Science, pages 82-93, 2021. doi: 10.1007/978-3-030-74432-8_6.

Source code in asymmetree/analysis/best_matches.py

def binary_explainable_triples(
    graph: nx.DiGraph,
    color_dict: dict | None = None,
) -> list[tuple]:
    """Extended informative triple set for binary-explainable graphs.

    Args:
        graph: A digraph whose nodes have the 'color' attribute.
        color_dict: A dict containing the colors that appear in the graph as keys and the list of
            nodes with the respective color as values. If not provided, this dict is computed
            internally.

    Returns:
        List of tuples where each item is a tuple (a, b, c) of three nodes such that ab|c is a
            triple displayed by every binary tree that explains the graph (if existent), see [1].

    References:
        1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Best match graphs with binary
           trees. In: C. Martín-Vide, M. A. Vega-Rodríguez, and T. Wheeler, editors, Algorithms for
           Computational Biology, 8th AlCoB, volume 12715 of Lecture Notes in Computer Science,
           pages 82-93, 2021. doi: 10.1007/978-3-030-74432-8_6.
    """
    if not isinstance(graph, nx.DiGraph):
        raise TypeError("must be a NetworkX 'Digraph'")

    if not color_dict:
        color_dict = sort_by_colors(graph)

    R_binary = []

    for c1, c2 in itertools.permutations(color_dict.keys(), 2):
        for a in color_dict[c1]:
            for b1, b2 in itertools.combinations(color_dict[c2], 2):
                if graph.has_edge(a, b1) and graph.has_edge(a, b2):
                    try:
                        if b2 < b1:
                            b1, b2 = b2, b1
                    except NotImplemented:
                        pass
                    R_binary.append((b1, b2, a))
                elif graph.has_edge(a, b1):
                    R_binary.append((a, b1, b2))
                elif graph.has_edge(a, b2):
                    R_binary.append((a, b2, b1))

    return R_binary

lrt_from_tree(T: Tree) -> Tree

Computes the Least Resolved Tree from a tree.

The unique Least Resolved Tree from a (pruned) gene tree is computed by contraction of all redundant edges.

Parameters:

Name	Type	Description	Default
T	Tree	A tree whose leaf nodes have the 'reconc' attribute.	required

Returns:

Type	Description
Tree	The unique least resolved tree (w.r.t. best matches), see [1].

References

1. D. Schaller, M. Geiß, P. F. Stadler, and M. Hellmuth. Complete characterization of incorrect orthology assignments in best match graphs. In: Journal of Mathematical Biology, 82:20, 2021. doi: 10.1007/s00285-021-01564-8.

Source code in asymmetree/analysis/best_matches.py

def lrt_from_tree(T: Tree) -> Tree:
    """Computes the Least Resolved Tree from a tree.

    The unique Least Resolved Tree from a (pruned) gene tree is computed by contraction of all
    redundant edges.

    Args:
        T: A tree whose leaf nodes have the 'reconc' attribute.

    Returns:
        The unique least resolved tree (w.r.t. best matches), see [1].

    References:
        1. D. Schaller, M. Geiß, P. F. Stadler, and M. Hellmuth. Complete characterization of
           incorrect orthology assignments in best match graphs. In: Journal of Mathematical
           Biology, 82:20, 2021. doi: 10.1007/s00285-021-01564-8.
    """
    lrt = T.copy()
    if not lrt.root:
        return lrt

    # remove planted root if existent
    if len(lrt.root.children) == 1:
        new_root = lrt.root.children[0]
        new_root.detach()
        lrt.root = new_root
        new_root.dist = 0.0

    # assign list of leaves to each node
    leaves = lrt.leaf_dict()

    subtree_reconcs = {}
    for v in lrt.preorder():
        subtree_reconcs[v] = {leaf.reconc for leaf in leaves[v]}

    arc_colors = _arc_colors(lrt, leaves, subtree_reconcs)
    redundant_edge_list = redundant_edges(lrt, subtree_reconcs, arc_colors)
    lrt.contract(redundant_edge_list)
    lrt = topology_only(lrt)

    return lrt

redundant_edges(T: Tree, subtree_reconcs: dict[TreeNode, set], arc_colors: dict[TreeNode, set]) -> list[tuple[TreeNode, TreeNode]]

Compute the redundant edges of a tree.

An inner edge (u,v) is redundant if there is no arc (x,y) in the BMG with lca(x,y) = u and y has a color that is not in the subtree of v, see [1].

Parameters:

Name	Type	Description	Default
T	Tree	A tree whose leaf nodes have the 'reconc' attribute.	required
subtree_reconcs	dict[TreeNode, set]	A dict mapping each node to the set of reconciliations of its leaf descendants.	required
arc_colors	dict[TreeNode, set]	A dict mapping each inner vertex to the set of colors s such that there is an arc (x,y) in the BMG with lca(x,y) = v and y has color s.	required

Returns:

Type	Description
list[tuple[TreeNode, TreeNode]]	A list of tuples representing the redundant edges of the tree.

Source code in asymmetree/analysis/best_matches.py

def redundant_edges(
    T: Tree,
    subtree_reconcs: dict[TreeNode, set],
    arc_colors: dict[TreeNode, set],
) -> list[tuple[TreeNode, TreeNode]]:
    """Compute the redundant edges of a tree.

    An inner edge (u,v) is redundant if there is no arc (x,y) in the BMG with lca(x,y) = u and y
    has a color that is not in the subtree of v, see [1].

    Args:
        T: A tree whose leaf nodes have the 'reconc' attribute.
        subtree_reconcs: A dict mapping each node to the set of reconciliations of its leaf
            descendants.
        arc_colors: A dict mapping each inner vertex to the set of colors s such that there is an
            arc (x,y) in the BMG with lca(x,y) = v and y has color s.

    Returns:
        A list of tuples representing the redundant edges of the tree.
    """
    redundant_edge_list = []

    for u, v in T.inner_edges():
        # colors s in sigma( L(T(u) \ T(v)) )
        aux_set = set()

        for v2 in u.children:
            if v2 is not v:
                aux_set.update(subtree_reconcs[v2])

        if not arc_colors[v].intersection(aux_set):
            redundant_edge_list.append((u, v))

    return redundant_edge_list

is_bmg(G: nx.DiGraph) -> Tree | None

Determine whether a colored digraph is a BMG.

If the graph is a BMG, then its LRT is returned, see [1].

Parameters:

Name	Type	Description	Default
G	DiGraph	A digraph whose nodes have the 'color' attribute.	required

Returns:

Type	Description
Tree | None	The least-resolved tree if the graph was a BMG, and None otherwise.

References

1. D. Schaller, M. Geiß, E. Chávez, M. González Laffitte, A. López Sánchez, B. M. R. Stadler, D. I. Valdivia, M. Hellmuth, M. Hernández Rosales, and P. F. Stadler. Corrigendum to 'Best match graphs'. In: Journal of Mathematical Biology, 82:47, 2021. doi: 10.1007/s00285-021-01601-6.

Source code in asymmetree/analysis/best_matches.py

def is_bmg(G: nx.DiGraph) -> Tree | None:
    """Determine whether a colored digraph is a BMG.

    If the graph is a BMG, then its LRT is returned, see [1].

    Args:
        G: A digraph whose nodes have the 'color' attribute.

    Returns:
        The least-resolved tree if the graph was a BMG, and None otherwise.

    References:
        1. D. Schaller, M. Geiß, E. Chávez, M. González Laffitte, A. López Sánchez, B. M. R.
           Stadler, D. I. Valdivia, M. Hellmuth, M. Hernández Rosales, and P. F. Stadler.
           Corrigendum to 'Best match graphs'. In: Journal of Mathematical Biology, 82:47, 2021.
           doi: 10.1007/s00285-021-01601-6.
    """
    if not isinstance(G, nx.DiGraph):
        raise TypeError("not a digraph")
    if not is_properly_colored(G):
        raise RuntimeError("not a properly colored digraph")

    subtrees = []
    colors = None

    for wcc in nx.weakly_connected_components(G):
        sg = G.subgraph(wcc).copy()
        color_dict = sort_by_colors(sg)

        # the color sets of all components must be equal
        if colors is None:
            colors = set(color_dict)
        elif colors != set(color_dict):
            return

        L = {v for v in sg.nodes()}
        R = informative_triples(sg, color_dict=color_dict)
        build = Build(R, L, mincut=False)
        subtree = build.build_tree()

        if not subtree:
            return
        else:
            # a digraph is a BMG iff its equal to the BMG of BUILD(R)
            reconstruct_reconc_from_graph(subtree, sg)
            if not graphs_equal(sg, bmg_from_tree(subtree)):
                return
            subtrees.append(subtree)

    if len(subtrees) == 1:
        root = subtrees[0].root
    else:
        root = TreeNode()
        for subtree in subtrees:
            root.add_child(subtree.root)

    return _finalize(root, G)

lrt_from_colored_graph(G: nx.DiGraph, mincut: bool = False, weighted_mincut: bool = False) -> Tree | None

Compute the least resolved tree from a colored digraph.

If the graph is a BMG, then the unique least resolved tree is returned. If not, then a heuristic least resolved tree is returned if mincut is True, and None otherwise.

Parameters:

Name	Type	Description	Default
G	DiGraph	A digraph whose nodes have the 'color' attribute.	required
mincut	bool	Handle inconsistencies that occur if the supplied graph is not a BMG.	False
weighted_mincut	bool	If True, apply a weighted mincut heuristic.	False

Returns:

Type	Description
Tree | None	The least resolved tree if the graph was a BMG, a heuristic least resolved tree if mincut is True, and None otherwise.

Source code in asymmetree/analysis/best_matches.py

def lrt_from_colored_graph(
    G: nx.DiGraph,
    mincut: bool = False,
    weighted_mincut: bool = False,
) -> Tree | None:
    """Compute the least resolved tree from a colored digraph.

    If the graph is a BMG, then the unique least resolved tree is returned. If not, then a heuristic
    least resolved tree is returned if `mincut` is True, and None otherwise.

    Args:
        G: A digraph whose nodes have the 'color' attribute.
        mincut: Handle inconsistencies that occur if the supplied graph is not a BMG.
        weighted_mincut: If True, apply a weighted mincut heuristic.

    Returns:
        The least resolved tree if the graph was a BMG, a heuristic least resolved tree if `mincut`
            is True, and None otherwise.
    """
    L = {v for v in G.nodes()}
    R = informative_triples(G)

    build = Build(R, L, mincut=mincut, weighted_mincut=weighted_mincut)
    tree = build.build_tree()

    return _finalize(tree, G)

correct_bmg(G: nx.DiGraph) -> nx.DiGraph

Build the LRT (using min cut in BUILD algorithm) and return its BMG.

Parameters:

Name	Type	Description	Default
G	DiGraph	A digraph whose nodes have the 'color' attribute.	required

Returns:

Type	Description
DiGraph	The graph edited to be a BMG.

Raises:

Type	Description
RuntimeError	If the mincut heuristic fails to construct a tree.

Source code in asymmetree/analysis/best_matches.py

def correct_bmg(G: nx.DiGraph) -> nx.DiGraph:
    """Build the LRT (using min cut in BUILD algorithm) and return its BMG.

    Args:
        G: A digraph whose nodes have the 'color' attribute.

    Returns:
        The graph edited to be a BMG.

    Raises:
        RuntimeError: If the mincut heuristic fails to construct a tree.
    """
    if G.order() == 0:
        return nx.DiGraph()

    subtrees = []
    for sg in (G.subgraph(c) for c in nx.weakly_connected_components(G)):
        tree = lrt_from_colored_graph(sg, mincut=True)
        if not tree:
            raise RuntimeError("mincut heuristic failed to construct a tree")

        subtrees.append(tree)

    if len(subtrees) == 1:
        tree = subtrees[0]
    else:
        tree = Tree(TreeNode())
        for subtree in subtrees:
            tree.root.add_child(subtree.root)

    return bmg_from_tree(tree)

lrt_from_2bmg(G: nx.DiGraph) -> Tree | None

Effieciently constructs the LRT for a 2-BMG via the support vertices.

Implementation of the efficient algorithm in [1].

Parameters:

Name	Type	Description	Default
G	DiGraph	A digraph whose nodes have the 'color' attribute.	required

Returns:

Type	Description
Tree | None	The least-resolved tree if the graph was a 2-colored BMG, and None otherwise.

References

1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Least resolved trees for two-colored best match graphs. In: Journal of Graph Algorithms and Applications, 25(1):397-416, 2021. doi: 10.7155/jgaa.00564.

Source code in asymmetree/analysis/best_matches.py

def lrt_from_2bmg(G: nx.DiGraph) -> Tree | None:
    """Effieciently constructs the LRT for a 2-BMG via the support vertices.

    Implementation of the efficient algorithm in [1].

    Args:
        G: A digraph whose nodes have the 'color' attribute.

    Returns:
        The least-resolved tree if the graph was a 2-colored BMG, and None otherwise.

    References:
        1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Least resolved trees for
           two-colored best match graphs. In: Journal of Graph Algorithms and Applications,
           25(1):397-416, 2021. doi: 10.7155/jgaa.00564.
    """
    builder = TwoColoredLRT(G)

    return builder.build_tree()

binary_refinable_tree(G: nx.DiGraph, mincut: bool = False, weighted_mincut: bool = False) -> Tree | None

Compute the binary-refinable tree (BRT) for a graph.

The BRT of a BMG can, if it exists, be refined arbitrarily and such that it still explains the BMG. In particular, all binary explanations are refinements of the BRT.

Parameters:

Name	Type	Description	Default
G	DiGraph	A digraph whose nodes have the color attribute.	required
mincut	bool	If True, handle inconsistencies that occur if the supplied graph is not a binary-explainable BMG.	False
weighted_mincut	bool	If True, apply a weighted mincut heuristic.	False

Returns:

Type	Description
Tree | None	The unique binary-resolvable tree if the graph is a binary-explainable BMG, see [1], or a heuristic tree if it is not but mincut is True, or None if neither is the case.

References

1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Best match graphs with binary trees. In: C. Martín-Vide, M. A. Vega-Rodríguez, and T. Wheeler, editors, Algorithms for Computational Biology, 8th AlCoB, volume 12715 of Lecture Notes in Computer Science, pages 82-93, 2021. doi: 10.1007/978-3-030-74432-8_6.

Source code in asymmetree/analysis/best_matches.py

def binary_refinable_tree(
    G: nx.DiGraph,
    mincut: bool = False,
    weighted_mincut: bool = False,
) -> Tree | None:
    """Compute the binary-refinable tree (BRT) for a graph.

    The BRT of a BMG can, if it exists, be refined arbitrarily and such that it still explains the
    BMG. In particular, all binary explanations are refinements of the BRT.

    Args:
        G: A digraph whose nodes have the `color` attribute.
        mincut: If True, handle inconsistencies that occur if the supplied graph is not a
            binary-explainable BMG.
        weighted_mincut: If True, apply a weighted mincut heuristic.

    Returns:
        The unique binary-resolvable tree if the graph is a binary-explainable BMG, see [1], or a
            heuristic tree if it is not but `mincut` is True, or None if neither is the case.

    References:
        1. D. Schaller, M. Geiß, M. Hellmuth, and P. F. Stadler. Best match graphs with binary
           trees. In: C. Martín-Vide, M. A. Vega-Rodríguez, and T. Wheeler, editors, Algorithms for
           Computational Biology, 8th AlCoB, volume 12715 of Lecture Notes in Computer Science,
           pages 82-93, 2021. doi: 10.1007/978-3-030-74432-8_6.
    """
    L = {v for v in G.nodes()}
    R = binary_explainable_triples(G)

    build = Build(R, L, mincut=mincut, weighted_mincut=weighted_mincut)
    brt = build.build_tree()

    return _finalize(brt, G)

augment_and_label(tree: Tree, inplace: bool = False) -> Tree

Augment tree and add event labeling based on color intersections.

Parameters:

Name	Type	Description	Default
tree	Tree	A tree whose nodes have the 'reconc' attribute.	required
inplace	bool	If True, the supplied tree instance is directly manipulated, otherwise the tree is copied first.	False

Returns:

Type	Description
Tree	The unique augmented and duplication/loss-labeled tree as defined in [1].

References

1. D. Schaller, M. Geiß, P. F. Stadler, and M. Hellmuth. Complete characterization of incorrect orthology assignments in best match graphs. In: Journal of Mathematical Biology, 82:20, 2021. doi: 10.1007/s00285-021-01564-8.

Source code in asymmetree/analysis/best_matches.py

def augment_and_label(tree: Tree, inplace: bool = False) -> Tree:
    """Augment tree and add event labeling based on color intersections.

    Args:
        tree: A tree whose nodes have the 'reconc' attribute.
        inplace: If True, the supplied tree instance is directly manipulated, otherwise the tree
            is copied first.

    Returns:
        The unique augmented and duplication/loss-labeled tree as defined in [1].

    References:
        1. D. Schaller, M. Geiß, P. F. Stadler, and M. Hellmuth. Complete characterization of
           incorrect orthology assignments in best match graphs. In: Journal of Mathematical
           Biology, 82:20, 2021. doi: 10.1007/s00285-021-01564-8.
    """
    if not inplace:
        tree = tree.copy()

    # precompute since tree will be modified
    inner_nodes = [u for u in tree.inner_nodes()]
    leaves = tree.leaf_dict()

    for u in inner_nodes:
        # should not appear except if tree is planted, otherwise ignore
        if len(u.children) == 1:
            continue

        C = _color_intersection_components(u, leaves)

        if len(C) == 1:
            u.event = "D"

        else:
            u.event = "S"

            for cc in C:
                if len(cc) > 1:
                    w = TreeNode(event="D")
                    u.add_child(w)
                    for v in cc:
                        u.remove_child(v)
                        w.add_child(v)

    return tree

classify_good_ugly(bmg: nx.DiGraph, rbmg: nx.Graph, fp: nx.Graph) -> nx.Graph

Classify edges in the false-positives graph as middle_in_good, first_in_ugly, first_in_bad.

The false positives in an RBMG (w.r.t. to orthology) can be classified into several categories, including middle edges of good quartets, first edges of ugly quartets, and first edges of bad quartets, see [1].

Parameters:

Name	Type	Description	Default
bmg	DiGraph	The best match graph.	required
rbmg	Graph	The reciprocal best match graph.	required
fp	Graph	The false-positives graph.	required

Returns:

Type	Description
Graph	The updated false-positives graph with edges classified as middle_in_good, first_in_ugly, and first_in_bad.

References

1. D. Schaller, M. Geiß, P. F. Stadler, and M. Hellmuth. Complete characterization of incorrect orthology assignments in best match graphs. In: Journal of Mathematical Biology, 82:20, 2021. doi: 10.1007/s00285-021-01564-8.

Source code in asymmetree/analysis/best_matches.py

def classify_good_ugly(bmg: nx.DiGraph, rbmg: nx.Graph, fp: nx.Graph) -> nx.Graph:
    """Classify edges in the false-positives graph as middle_in_good, first_in_ugly, first_in_bad.

    The false positives in an RBMG (w.r.t. to orthology) can be classified into several categories,
    including middle edges of good quartets, first edges of ugly quartets, and first edges of bad
    quartets, see [1].

    Args:
        bmg: The best match graph.
        rbmg: The reciprocal best match graph.
        fp: The false-positives graph.

    Returns:
        The updated false-positives graph with edges classified as `middle_in_good`,
            `first_in_ugly`, and `first_in_bad`.

    References:
        1. D. Schaller, M. Geiß, P. F. Stadler, and M. Hellmuth. Complete characterization of
           incorrect orthology assignments in best match graphs. In: Journal of Mathematical
           Biology, 82:20, 2021. doi: 10.1007/s00285-021-01564-8.
    """
    for x, y in fp.edges():
        fp.edges[x, y]["middle_in_good"] = False
        fp.edges[x, y]["first_in_ugly"] = False
        fp.edges[x, y]["first_in_bad"] = False

    color_dict = sort_by_colors(bmg)

    for x, y in fp.edges():
        # middle edges of good quartets
        for c in color_dict:
            if c == bmg.nodes[x]["color"] or c == bmg.nodes[y]["color"]:
                continue

            for z1, z2 in itertools.permutations(color_dict[c], 2):
                if (
                    rbmg.has_edge(z1, x)
                    and rbmg.has_edge(z2, y)
                    and bmg.has_edge(z1, y)
                    and not bmg.has_edge(y, z1)
                    and bmg.has_edge(z2, x)
                    and not bmg.has_edge(x, z2)
                ):
                    fp.edges[x, y]["middle_in_good"] = True

        # first edges of ugly/bad quartets
        for _ in range(2):
            for x2 in color_dict[bmg.nodes[x]["color"]]:
                if x == x2:
                    continue

                for z in rbmg.neighbors(x2):
                    if (
                        bmg.nodes[z]["color"] == bmg.nodes[x]["color"]
                        or bmg.nodes[z]["color"] == bmg.nodes[y]["color"]
                    ):
                        continue
                    # first edges of ugly quartets
                    if rbmg.has_edge(y, x2) and not rbmg.has_edge(z, y) and not rbmg.has_edge(z, x):
                        fp.edges[x, y]["first_in_ugly"] = True

                    # first edges of bad quartets
                    elif (
                        not rbmg.has_edge(y, x2) and rbmg.has_edge(z, y) and not bmg.has_edge(x, z)
                    ):
                        fp.edges[x, y]["first_in_bad"] = True

            # swap for second iteration
            x, y = y, x

    return fp

count_good_ugly(bmg: nx.DiGraph, rbmg: nx.Graph, fp: nx.Graph) -> tuple[int, int, int]

Count number of edge types in the false-positives graph.

Parameters:

Name	Type	Description	Default
bmg	DiGraph	The best match graph.	required
rbmg	Graph	The reciprocal best match graph.	required
fp	Graph	The false-positives graph.	required

Returns:

Type	Description
tuple[int, int, int]	A tuple containing the number of edges in the false-positives graph that are middle edges of good quartets, first edges of ugly quartets, and both.

Source code in asymmetree/analysis/best_matches.py

def count_good_ugly(bmg: nx.DiGraph, rbmg: nx.Graph, fp: nx.Graph) -> tuple[int, int, int]:
    """Count number of edge types in the false-positives graph.

    Args:
        bmg: The best match graph.
        rbmg: The reciprocal best match graph.
        fp: The false-positives graph.

    Returns:
        A tuple containing the number of edges in the false-positives graph that are middle edges of
            good quartets, first edges of ugly quartets, and both.
    """
    classify_good_ugly(bmg, rbmg, fp)

    good, ugly, both = 0, 0, 0

    for x, y in fp.edges():
        if fp.edges[x, y]["middle_in_good"] and fp.edges[x, y]["first_in_ugly"]:
            good += 1
            ugly += 1
            both += 1
        elif fp.edges[x, y]["middle_in_good"]:
            good += 1
        elif fp.edges[x, y]["first_in_ugly"]:
            ugly += 1

    return good, ugly, both

count_good_ugly_bad(bmg: nx.DiGraph, rbmg: nx.Graph, fp: nx.Graph) -> tuple[int, int, int, int, int, int, int]

Count number of edge types in the false-positives graph.

Parameters:

Name	Type	Description	Default
bmg	DiGraph	The best match graph.	required
rbmg	Graph	The reciprocal best match graph.	required
fp	Graph	The false-positives graph.	required

Returns:

Type	Description
tuple[int, int, int, int, int, int, int]	A tuple containing the number of edges in the false-positives graph that are middle edges of good quartets, first edges of ugly quartets, first edges of bad quartets, and all combinations thereof in in order: (gub, gu, gb, ub, g, u, b).

Source code in asymmetree/analysis/best_matches.py

def count_good_ugly_bad(
    bmg: nx.DiGraph,
    rbmg: nx.Graph,
    fp: nx.Graph,
) -> tuple[int, int, int, int, int, int, int]:
    """Count number of edge types in the false-positives graph.

    Args:
        bmg: The best match graph.
        rbmg: The reciprocal best match graph.
        fp: The false-positives graph.

    Returns:
        A tuple containing the number of edges in the false-positives graph that are middle edges of
            good quartets, first edges of ugly quartets, first edges of bad quartets, and all
            combinations thereof in in order: (gub, gu, gb, ub, g, u, b).
    """
    classify_good_ugly(bmg, rbmg, fp)

    gub, gu, gb, ub, g, u, b = 0, 0, 0, 0, 0, 0, 0

    for x, y in fp.edges():
        if (
            fp.edges[x, y]["middle_in_good"]
            and fp.edges[x, y]["first_in_ugly"]
            and fp.edges[x, y]["first_in_bad"]
        ):
            gub += 1
            gu += 1
            gb += 1
            ub += 1
            g += 1
            u += 1
            b += 1

        elif fp.edges[x, y]["middle_in_good"] and fp.edges[x, y]["first_in_ugly"]:
            gu += 1
            g += 1
            u += 1
        elif fp.edges[x, y]["middle_in_good"] and fp.edges[x, y]["first_in_bad"]:
            gb += 1
            g += 1
            b += 1
        elif fp.edges[x, y]["first_in_ugly"] and fp.edges[x, y]["first_in_bad"]:
            ub += 1
            u += 1
            b += 1
        elif fp.edges[x, y]["middle_in_good"]:
            g += 1
        elif fp.edges[x, y]["first_in_ugly"]:
            u += 1
        elif fp.edges[x, y]["first_in_bad"]:
            b += 1

    return gub, gu, gb, ub, g, u, b

list_path(t: int, k: int, path: list, G: nx.Graph, P4_list: list) -> None

List all paths of length k in G.

The paths are stored in P4_list as tuples of the form (v1, v2, v3, v4).

Parameters:

Name	Type	Description	Default
t	int	The current length of the path.	required
k	int	The desired length of the path.	required
path	list	The current path being explored.	required
G	Graph	The graph in which to find paths.	required
P4_list	list	A list to store the found paths.	required

Source code in asymmetree/analysis/best_matches.py

def list_path(t: int, k: int, path: list, G: nx.Graph, P4_list: list) -> None:
    """List all paths of length k in G.

    The paths are stored in P4_list as tuples of the form (v1, v2, v3, v4).

    Args:
        t: The current length of the path.
        k: The desired length of the path.
        path: The current path being explored.
        G: The graph in which to find paths.
        P4_list: A list to store the found paths.
    """
    for u in G[path[-1]]:
        G.nodes[u]["vis"] += 1

    for u in G[path[-1]]:
        if G.nodes[u]["vis"] == 1:
            if t < k:
                list_path(t + 1, k, path + [u], G, P4_list)
            elif path[0] < u:
                P4_list.append(tuple(path) + (u,))

    for u in G[path[-1]]:
        G.nodes[u]["vis"] -= 1

find_all_P4(G: nx.Graph) -> list[tuple]

Find all induced P4 in G.

Parameters:

Name	Type	Description	Default
G	Graph	A graph in which to find induced P4.	required

Returns:

Type	Description
list[tuple]	A list of tuples of the form (v1, v2, v3, v4) representing all induced P4 in G.

Source code in asymmetree/analysis/best_matches.py

def find_all_P4(G: nx.Graph) -> list[tuple]:
    """Find all induced P4 in G.

    Args:
        G: A graph in which to find induced P4.

    Returns:
        A list of tuples of the form (v1, v2, v3, v4) representing all induced P4 in G.
    """
    P4_list = []

    for v in G:
        G.nodes[v]["vis"] = 0

    for v in G:
        G.nodes[v]["vis"] = 1
        list_path(2, 4, [v], G, P4_list)
        G.nodes[v]["vis"] = 0

    return P4_list

is_good_quartet(path: tuple, bmg: nx.Graph) -> bool

Determine whether an induced P4 is a good quartet in the bmg.

Parameters:

Name	Type	Description	Default
path	tuple	A tuple representing the induced P4.	required
bmg	Graph	The best match graph.	required

Returns:

Type	Description
bool	True if the induced P4 is a good quartet, False otherwise.

Source code in asymmetree/analysis/best_matches.py

def is_good_quartet(path: tuple, bmg: nx.Graph) -> bool:
    """Determine whether an induced P4 is a good quartet in the bmg.

    Args:
        path: A tuple representing the induced P4.
        bmg: The best match graph.

    Returns:
        True if the induced P4 is a good quartet, False otherwise.
    """
    if (
        bmg.has_edge(path[0], path[2])
        and bmg.has_edge(path[3], path[1])
        and bmg.nodes[path[0]]["color"] == bmg.nodes[path[3]]["color"]
    ):
        return True
    else:
        return False

remove_all_P4(rbmg: nx.Graph, bmg: nx.Graph, P4_list: list[tuple] = None) -> nx.Graph

Find all good quartets and removes the middle edges.

Parameters:

Name	Type	Description	Default
rbmg	Graph	The reciprocal best match graph.	required
bmg	Graph	The best match graph.	required
P4_list	list[tuple]	A list of tuples representing all induced P4 in the rbmg. If not provided, it is computed internally.	None

Returns:

Type	Description
Graph	The graph obtained from the RBMG by removing the middle edges of all good quartets.

Source code in asymmetree/analysis/best_matches.py

def remove_all_P4(rbmg: nx.Graph, bmg: nx.Graph, P4_list: list[tuple] = None) -> nx.Graph:
    """Find all good quartets and removes the middle edges.

    Args:
        rbmg: The reciprocal best match graph.
        bmg: The best match graph.
        P4_list: A list of tuples representing all induced P4 in the rbmg. If not provided, it is
            computed internally.

    Returns:
        The graph obtained from the RBMG by removing the middle edges of all good quartets.
    """
    graph = rbmg.copy()

    if P4_list is None:
        P4_list = find_all_P4(graph)

    for path in P4_list:
        if is_good_quartet(path, bmg) and graph.has_edge(path[1], path[2]):
            graph.remove_edge(path[1], path[2])

    return graph

find_all_C6(G: nx.Graph, P4_list: list[tuple] | None = None) -> list[tuple]

Find induces C6 of the form .

Parameters:

Name	Type	Description	Default
G	Graph	A graph in which to find induced C6.	required
P4_list	list[tuple] | None	A list of tuples representing all induced P4 in G. If not provided, it is computed internally.	None

Returns:

Type	Description
list[tuple]	A list of tuples of the form (x1, y1, z1, x2, y2, z2) representing all induced C6 in G.

Raises:

Type	Description
RuntimeError	If a detected C6 does not satisfy the properties of an induced C6 in G.

Source code in asymmetree/analysis/best_matches.py

def find_all_C6(G: nx.Graph, P4_list: list[tuple] | None = None) -> list[tuple]:
    """Find induces C6 of the form <x1 y1 z1 x2 y2 z2>.

    Args:
        G: A graph in which to find induced C6.
        P4_list: A list of tuples representing all induced P4 in G. If not provided, it is computed
            internally.

    Returns:
        A list of tuples of the form (x1, y1, z1, x2, y2, z2) representing all induced C6 in G.

    Raises:
        RuntimeError: If a detected C6 does not satisfy the properties of an induced C6 in G.
    """
    if P4_list is None:
        P4_list = find_all_P4(G)

    endpoints: dict[tuple, dict[tuple, list[tuple]]] = {}
    for path in P4_list:
        a, b, c, d = path if path[0] < path[3] else path[::-1]
        a_col, b_col, c_col, d_col = (
            G.nodes[a]["color"],
            G.nodes[b]["color"],
            G.nodes[c]["color"],
            G.nodes[d]["color"],
        )

        # avoid multiple recording for the 3 different colors
        if a_col == d_col and a_col < b_col and a_col < c_col:
            if (a, d) not in endpoints:
                endpoints[(a, d)] = {}
            if (b_col, c_col) not in endpoints[(a, d)]:
                endpoints[(a, d)][(b_col, c_col)] = []
            endpoints[(a, d)][(b_col, c_col)].append(path)

    C6_list = []
    for colors in endpoints.values():
        for color_pair, paths in colors.items():
            y_color, z_color = color_pair

            # avoid multiple recording for the 2 directions:
            if y_color >= z_color or (z_color, y_color) not in colors:
                continue

            for path1 in paths:
                for path2 in colors[(z_color, y_color)]:
                    if (not G.has_edge(path1[1], path2[1])) and (
                        not G.has_edge(path1[2], path2[2])
                    ):
                        C6_list.append(path1 + (path2[2], path2[1]))
                        if not _validate_C6(G, C6_list[-1]):
                            raise RuntimeError("problem in C6 detection")

    return C6_list

graph_type(G: nx.Graph) -> tuple[str, list[tuple], list[tuple]]

Classify the graph as type A, B, or C based on the presence of induced P4 and C6.

A graph is of type A if it contains no induced P4, of type B if it contains an induced P4 but no induced C6, and of type C if it contains an induced C6.

Parameters:

Name	Type	Description	Default
G	Graph	A graph in which to classify the structure.	required

Returns:

Type	Description
str	The type of the graph ("A", "B", or "C").
list[tuple]	A list of all induced P4 in the graph.
list[tuple]	A list of all induced C6 in the graph.

Source code in asymmetree/analysis/best_matches.py

def graph_type(G: nx.Graph) -> tuple[str, list[tuple], list[tuple]]:
    """Classify the graph as type A, B, or C based on the presence of induced P4 and C6.

    A graph is of type A if it contains no induced P4, of type B if it contains an induced P4 but no
    induced C6, and of type C if it contains an induced C6.

    Args:
        G: A graph in which to classify the structure.

    Returns:
        The type of the graph ("A", "B", or "C").
        A list of all induced P4 in the graph.
        A list of all induced C6 in the graph.
    """
    P4_list = find_all_P4(G)
    C6_list = find_all_C6(G, P4_list=P4_list)

    if C6_list:
        return "C", P4_list, C6_list
    elif P4_list:
        return "B", P4_list, C6_list
    else:
        return "A", P4_list, C6_list

horizontal_gene_transfer

Horizontal gene transfer (HGT), xenology and Fitch graphs.

This module contains several functions for the analysis of horizontal gene transfer events in the simulated scenarios. In particular, the directed and undirected Fitch graph can be extracted, as well as the pairs of genes that diverged later than the respective species in which they reside, i.e., the later-divergence-time (LDT) graph, see [1]. The latter situation is indicative of the presence of HGT events in the scenario.

References

1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021, 83(1):10. doi: 10.1007/s00285-021-01631-0.

RsScenarioConstructor

Construct an rs-scenario for a given LDT graph.

Implementation of the algorithm presented in [1].

References

1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021, 83(1):10. doi: 10.1007/s00285-021-01631-0.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

class RsScenarioConstructor:
    """Construct an rs-scenario for a given LDT graph.

    Implementation of the algorithm presented in [1].

    References:
        1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect
           identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021,
           83(1):10. doi: 10.1007/s00285-021-01631-0.
    """

    def __init__(self, colored_cograph: nx.Graph, color_set: set | None = None) -> None:
        """Construct an instance of the RsScenarioConstructor class.

        Args:
            colored_cograph: A cograph whose nodes have the 'color' attribute.
            color_set: The set of colors appearing in the graph. The default is None, in which case
                such a set is created internally.
        """
        self.G = colored_cograph

        self.L = [x for x in self.G.nodes()]

        if color_set is None:
            color_set = set()
            for v in self.L:
                color_set.add(self.G.nodes[v]["color"])
        elif not isinstance(color_set, (set, dict)):
            color_set = set(color_set)
        self.color_set = color_set

    def run(self) -> tuple[Tree, Tree] | None:
        """Construct an rs-scenario consisting of a species and gene tree.

        Returns:
            The species and gene tree, or None if the supplied graph was not an LDT graph.
        """
        self.S = self._species_tree()
        if not self.S:
            return

        self.epsilon = float("inf")
        for u, v in self.S.edges():
            self.epsilon = min(self.epsilon, abs(u.tstamp - v.tstamp))
        self.epsilon /= 3.0

        self.S_leaves = self.S.leaf_dict()

        # top-level call on full leaf set and rho_S
        self.T = Tree(self._build_gene_tree(self.L, self.S.root.children[0]))

        # --- post-processing of T ---

        # planted root
        add_planted_root(self.T)
        self.T.root.tstamp = self.S.root.tstamp
        self.T.root.reconc = self.S.root.label

        # suppress all vertices with a single child except the planted root
        to_suppress = []
        for u in self.T.preorder():
            if u.parent and len(u.children) == 1:
                to_suppress.append((u.parent, u))
        self.T.contract(to_suppress)

        distance_from_timing(self.T)

        return self.S, self.T

    def _species_tree(self) -> Tree | None:
        """Construct the species tree from the input graph.

        Returns:
            The species tree, or None if the supplied graph was not an LDT graph.
        """
        cotree = to_cotree(self.G)

        if not cotree:
            return

        triples = set()

        for a, b, c in paths_of_length_2(cotree):
            A = self.G.nodes[a.label]["color"]
            B = self.G.nodes[b.label]["color"]
            C = self.G.nodes[c.label]["color"]

            if A != B and A != C and B != C:
                # sort to avoid redundancy in triple set
                if A > C:
                    A, C = C, A

                triples.add((A, C, B))

        build = Build(triples, self.color_set, mincut=False)
        S = build.build_tree()

        if not S:
            return

        phylo_tree_attributes(S, inplace=True)
        add_planted_root(S)
        assign_missing_labels(S)
        _simulate_timing(S)
        distance_from_timing(S)

        return S

    def _build_gene_tree(self, L: set, u_S: TreeNode) -> TreeNode:
        """Build the gene tree recursively.

        Args:
            L: The set of leaf labels that shall be included in the gene tree.
            u_S: The current vertex in the species tree that is being processed.

        Returns:
            The root of the gene tree for the current vertex u_S and leaf set L.
        """
        u_T = TreeNode(
            label="",
            event="D",
            reconc=(u_S.parent.label, u_S.label),
            tstamp=u_S.tstamp + self.epsilon,
        )

        # u_S is a leaf
        if not u_S.children:
            for x in L:
                leaf = TreeNode(label=x, reconc=self.G.nodes[x]["color"], tstamp=0.0)
                u_T.add_child(leaf)

            return u_T

        # otherwise u_S is an inner vertex

        # maps color to the respective child of u_S
        color_to_v_S = {}
        for v_S in u_S.children:
            for x in self.S_leaves[v_S]:
                color_to_v_S[x.label] = v_S

        # connected components C of G[L']
        for C in nx.connected_components(self.G.subgraph(L)):
            # choose v_S_star s.t. L(S(v_S_star)) \cap sigma(C) is non-empty
            v_S_star = color_to_v_S[self.G.nodes[next(iter(C))]["color"]]
            v_T = TreeNode(
                label="",
                event="H",
                reconc=(u_S.label, v_S_star.label),
                tstamp=u_S.tstamp - self.epsilon,
            )
            u_T.add_child(v_T)

            # aux. graph for equivalence classes
            aux_graph = nx.Graph()
            aux_graph.add_nodes_from(C)

            for x, y in itertools.combinations(C, 2):
                if color_to_v_S[self.G.nodes[x]["color"]] is color_to_v_S[self.G.nodes[y]["color"]]:
                    aux_graph.add_edge(x, y)

            # for each equivalence class K that is a subset of C
            for K in nx.connected_components(aux_graph):
                # choose v_S s.t. \sigma(K) \subseteq L(S(v_S))
                v_S = color_to_v_S[self.G.nodes[next(iter(K))]["color"]]

                w_K = self._build_gene_tree(K, v_S)
                v_T.add_child(w_K)

                if v_S is not v_S_star:
                    w_K.transferred = 1

        return u_T

__init__(colored_cograph: nx.Graph, color_set: set | None = None) -> None

Construct an instance of the RsScenarioConstructor class.

Parameters:

Name	Type	Description	Default
colored_cograph	Graph	A cograph whose nodes have the 'color' attribute.	required
color_set	set | None	The set of colors appearing in the graph. The default is None, in which case such a set is created internally.	None

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def __init__(self, colored_cograph: nx.Graph, color_set: set | None = None) -> None:
    """Construct an instance of the RsScenarioConstructor class.

    Args:
        colored_cograph: A cograph whose nodes have the 'color' attribute.
        color_set: The set of colors appearing in the graph. The default is None, in which case
            such a set is created internally.
    """
    self.G = colored_cograph

    self.L = [x for x in self.G.nodes()]

    if color_set is None:
        color_set = set()
        for v in self.L:
            color_set.add(self.G.nodes[v]["color"])
    elif not isinstance(color_set, (set, dict)):
        color_set = set(color_set)
    self.color_set = color_set

run() -> tuple[Tree, Tree] | None

Construct an rs-scenario consisting of a species and gene tree.

Returns:

Type	Description
tuple[Tree, Tree] | None	The species and gene tree, or None if the supplied graph was not an LDT graph.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def run(self) -> tuple[Tree, Tree] | None:
    """Construct an rs-scenario consisting of a species and gene tree.

    Returns:
        The species and gene tree, or None if the supplied graph was not an LDT graph.
    """
    self.S = self._species_tree()
    if not self.S:
        return

    self.epsilon = float("inf")
    for u, v in self.S.edges():
        self.epsilon = min(self.epsilon, abs(u.tstamp - v.tstamp))
    self.epsilon /= 3.0

    self.S_leaves = self.S.leaf_dict()

    # top-level call on full leaf set and rho_S
    self.T = Tree(self._build_gene_tree(self.L, self.S.root.children[0]))

    # --- post-processing of T ---

    # planted root
    add_planted_root(self.T)
    self.T.root.tstamp = self.S.root.tstamp
    self.T.root.reconc = self.S.root.label

    # suppress all vertices with a single child except the planted root
    to_suppress = []
    for u in self.T.preorder():
        if u.parent and len(u.children) == 1:
            to_suppress.append((u.parent, u))
    self.T.contract(to_suppress)

    distance_from_timing(self.T)

    return self.S, self.T

true_transfer_edges(T: Tree) -> set[TreeNode]

Returns a set containing v if (u, v) is labeled as a transfer edge.

Parameters:

Name	Type	Description	Default
T	Tree	A tree whose nodes have the transferred attribute.	required

Returns:

Type	Description
set[TreeNode]	The subset of TreeNode instances in the tree for which transferred is 1.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def true_transfer_edges(T: Tree) -> set[TreeNode]:
    """Returns a set containing v if (u, v) is labeled as a transfer edge.

    Args:
        T: A tree whose nodes have the `transferred` attribute.

    Returns:
        The subset of `TreeNode` instances in the tree for which `transferred` is 1.
    """
    return {v for _, v in T.edges() if v.transferred}

rs_transfer_edges(T: Tree, S: Tree, lca_S: LCA = None) -> set[TreeNode]

Transfer edges in T according to the relaxed scenario definition.

An edge (u,v) in T is an (rs-)transfer edge if u and v are mapped to incomparable nodes/edges in the species tree S.

Parameters:

Name	Type	Description	Default
T	Tree	The gene tree.	required
S	Tree	The species tree corresponding to the gene tree.	required
lca_S	LCA	Instance of LCA corresponding to the species tree, default is None in which case a new instance is created and used.	None

Returns:

Type	Description
set[TreeNode]	The subset of TreeNode instances in the tree for which the reconc attributes of its incident edges refer to incomparable branches in the species tree.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def rs_transfer_edges(T: Tree, S: Tree, lca_S: LCA = None) -> set[TreeNode]:
    """Transfer edges in T according to the relaxed scenario definition.

    An edge (u,v) in T is an (rs-)transfer edge if u and v are mapped to incomparable nodes/edges in
    the species tree S.

    Args:
        T: The gene tree.
        S: The species tree corresponding to the gene tree.
        lca_S: Instance of LCA corresponding to the species tree, default is None
            in which case a new instance is created and used.

    Returns:
        The subset of `TreeNode` instances in the tree for which the `reconc` attributes of its
            incident edges refer to incomparable branches in the species tree.
    """
    if not isinstance(lca_S, LCA):
        lca_S = LCA(S)

    transfer_edges = set()

    for u, v in T.edges():
        if not lca_S.are_comparable(u.reconc, v.reconc):
            transfer_edges.add(v)

    return transfer_edges

fitch(tree: Tree, transfer_edges: Iterable[TreeNode], supply_undirected: bool = False, lca_T: LCA | None = None) -> nx.DiGraph | tuple[nx.DiGraph, nx.Graph]

Returns the (directed) Fitch graph.

Parameters:

Name	Type	Description	Default
tree	Tree	The gene tree.	required
transfer_edges	Iterable[TreeNode]	The subset of TreeNode instances in the tree that are transfer edges.	required
supply_undirected	bool	If True, additionally return the undirected Fitch graph.	False
lca_T	LCA | None	Instance of LCA corresponding to the tree, default is False in which case a new instance is created and used.	None

Returns:

Type	Description
DiGraph | tuple[DiGraph, Graph]	The directed Fitch graph or a tuple containing the directed Fitch graph (networkx.DiGraph) and the undirected Fitch graph (networkx.Graph).

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def fitch(
    tree: Tree,
    transfer_edges: Iterable[TreeNode],
    supply_undirected: bool = False,
    lca_T: LCA | None = None,
) -> nx.DiGraph | tuple[nx.DiGraph, nx.Graph]:
    """Returns the (directed) Fitch graph.

    Args:
        tree: The gene tree.
        transfer_edges: The subset of TreeNode instances in the tree that are transfer edges.
        supply_undirected: If True, additionally return the undirected Fitch graph.
        lca_T: Instance of LCA corresponding to the tree, default is False in which case a new
            instance is created and used.

    Returns:
        The directed Fitch graph or a tuple containing the directed Fitch graph (`networkx.DiGraph`)
            and the undirected Fitch graph (`networkx.Graph`).
    """
    if not isinstance(lca_T, LCA):
        lca_T = LCA(tree)

    if not isinstance(transfer_edges, (set, dict)):
        transfer_edges = set(transfer_edges)

    leaves = [leaf for leaf in tree.leaves()]
    fitch = nx.DiGraph()

    # store for each leaf the first transfer edge on the way to the root
    first_transfer = {}

    for x in leaves:
        fitch.add_node(x.label, color=x.reconc)

        current = x
        while current:
            if current in transfer_edges:
                first_transfer[x] = current
                break
            current = current.parent

    for x, y in itertools.permutations(leaves, 2):
        if y in first_transfer and lca_T.ancestor_not_equal(lca_T(x, y), first_transfer[y]):
            fitch.add_edge(x.label, y.label)

    if not supply_undirected:
        return fitch
    else:
        return fitch, fitch.to_undirected()

undirected_fitch(tree: Tree, transfer_edges: Iterable[TreeNode], lca_T: LCA | None = None) -> nx.Graph

Returns the undirected Fitch graph.

Parameters:

Name	Type	Description	Default
tree	Tree	The gene tree.	required
transfer_edges	Iterable[TreeNode]	The subset of TreeNode instances in the tree that are transfer edges.	required
lca_T	LCA | None	Instance of LCA corresponding to the tree, default is None in which case a new instance is created and used.	None

Returns:

Type	Description
Graph	The undirected Fitch graph (networkx.Graph).

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def undirected_fitch(
    tree: Tree,
    transfer_edges: Iterable[TreeNode],
    lca_T: LCA | None = None,
) -> nx.Graph:
    """Returns the undirected Fitch graph.

    Args:
        tree: The gene tree.
        transfer_edges: The subset of TreeNode instances in the tree that are transfer edges.
        lca_T: Instance of LCA corresponding to the tree, default is None in which case a new
            instance is created and used.

    Returns:
        The undirected Fitch graph (`networkx.Graph`).
    """
    return fitch(tree, transfer_edges, supply_undirected=True, lca_T=lca_T)[1]

is_rs_fitch(G: nx.Graph, color_set: set | None = None) -> bool

Checks whether a given graph is an rs-Fitch graph.

This function checks whether the graph is the Fitch graph of some relaxed scenario (rs), see [1].

Parameters:

Name	Type	Description	Default
G	Graph	A graph whose nodes have the color attribute.	required
color_set	set | None	The set of colors appearing in the graph. The default is None, in which case such a set is created internally.	None

Returns:

Type	Description
bool	True if the graph is an rs-Fitch graph as defined in [1].

References

1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021, 83(1):10. doi: 10.1007/s00285-021-01631-0.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def is_rs_fitch(G: nx.Graph, color_set: set | None = None) -> bool:
    """Checks whether a given graph is an rs-Fitch graph.

    This function checks whether the graph is the Fitch graph of some relaxed scenario (rs), see
    [1].

    Args:
        G: A graph whose nodes have the `color` attribute.
        color_set: The set of colors appearing in the graph. The default is None, in which case such
            a set is created internally.

    Returns:
        True if the graph is an rs-Fitch graph as defined in [1].

    References:
        1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect
           identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021,
           83(1):10. doi: 10.1007/s00285-021-01631-0.
    """
    if color_set is None:
        color_set = set()
        for v in G.nodes():
            color_set.add(G.nodes[v]["color"])

    ind_sets = independent_sets(G)

    # not a complete multipartite graph
    if ind_sets is False:
        return False

    # trivial number of independent sets
    k = len(ind_sets)
    if k <= 1:
        return True

    # more than one independent set
    for ignore_set in range(-1, k):
        # ignore_set = -1 means that no set is ignored
        aux_graph = _rs_fitch_aux_graph(G, color_set, ind_sets, ignore_set)
        if not nx.is_connected(aux_graph):
            return True

    # if nothing was True
    return False

below_equal_above(T: Tree, S: Tree, lca_T: LCA | None = None, lca_S: LCA | None = None) -> tuple[nx.Graph, nx.Graph, nx.Graph]

Determine pairs of genes that diverged later, at the same time, or earlier as their species.

Returns three graphs with the leaves of the gene tree T as vertex set, and edges ab if and only if a and b diverged later, at the same time, or earlier, resp., than the corresponding species A and B in the species tree S.

Parameters:

Name	Type	Description	Default
T	Tree	The gene tree.	required
S	Tree	The corresponding species tree.	required
lca_T	LCA | None	Instance of LCA corresponding to the gene tree, default is None in which case a new instance is created and used.	None
lca_S	LCA | None	Instance of LCA corresponding to the species tree, default is None in which case a new instance is created and used.	None

Returns:

Type	Description
tuple[Graph, Graph, Graph]	The below, equal, and above relation as described above.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def below_equal_above(
    T: Tree,
    S: Tree,
    lca_T: LCA | None = None,
    lca_S: LCA | None = None,
) -> tuple[nx.Graph, nx.Graph, nx.Graph]:
    """Determine pairs of genes that diverged later, at the same time, or earlier as their species.

    Returns three graphs with the leaves of the gene tree T as vertex set, and edges ab if and only
    if a and b diverged later, at the same time, or earlier, resp., than the corresponding species
    A and B in the species tree S.

    Args:
        T: The gene tree.
        S: The corresponding species tree.
        lca_T: Instance of LCA corresponding to the gene tree, default is None in which case a new
            instance is created and used.
        lca_S: Instance of LCA corresponding to the species tree, default is None in which case a
            new instance is created and used.

    Returns:
        The below, equal, and above relation as described above.
    """
    L_T = [leaf for leaf in T.leaves()]
    L_S = {leaf.label: leaf for leaf in S.leaves()}

    below = nx.Graph()
    equal = nx.Graph()
    above = nx.Graph()
    for u in L_T:
        below.add_node(u.label, color=u.reconc)
        equal.add_node(u.label, color=u.reconc)
        above.add_node(u.label, color=u.reconc)

    if not isinstance(lca_T, LCA):
        lca_T = LCA(T)
    if not isinstance(lca_S, LCA):
        lca_S = LCA(S)

    for a, b in itertools.combinations(L_T, 2):
        t_ab = lca_T.get(a, b).tstamp
        t_AB = lca_S.get(L_S[a.reconc], L_S[b.reconc]).tstamp

        if t_ab < t_AB:
            below.add_edge(a.label, b.label)
        elif t_ab == t_AB:
            equal.add_edge(a.label, b.label)
        else:
            above.add_edge(a.label, b.label)

    return below, above, equal

ldt_graph(T: Tree, S: Tree, lca_T: LCA | None = None, lca_S: LCA | None = None) -> nx.Graph

Later-divergence-time graph, see [1].

Returns a graph with the leaves of the gene tree T as vertex set, and edges ab if and only if a and b diverged later than the corresponding species A and B in the species tree S.

Parameters:

Name	Type	Description	Default
T	Tree	The gene tree.	required
S	Tree	The corresponding species tree.	required
lca_T	LCA | None	Instance of LCA corresponding to the gene tree, default is None in which case a new instance is created and used.	None
lca_S	LCA | None	Instance of LCA corresponding to the species tree, default is None in which case a new instance is created and used.	None

References

1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021, 83(1):10. doi: 10.1007/s00285-021-01631-0.

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def ldt_graph(T: Tree, S: Tree, lca_T: LCA | None = None, lca_S: LCA | None = None) -> nx.Graph:
    """Later-divergence-time graph, see [1].

    Returns a graph with the leaves of the gene tree T as vertex set, and edges ab if and only if a
    and b diverged later than the corresponding species A and B in the species tree S.

    Args:
        T: The gene tree.
        S: The corresponding species tree.
        lca_T: Instance of LCA corresponding to the gene tree, default is None in which case a new
            instance is created and used.
        lca_S: Instance of LCA corresponding to the species tree, default is None in which case a
            new instance is created and used.

    References:
        1. D. Schaller, M. Lafond, P.F. Stadler, N. Wieseke, M. Hellmuth. Indirect
           identification of horizontal gene transfer. In: Journal of Mathematical Biology, 2021,
           83(1):10. doi: 10.1007/s00285-021-01631-0.
    """
    return below_equal_above(T, S, lca_T=lca_T, lca_S=lca_S)[0]

is_compatible(T: Tree, partition: list | tuple, lca: LCA | None = None) -> tuple[dict, list] | None

Checks whether a tree and a partition are compatible.

A tree with leaf set L and a partition of L are compatible if there is a subset of tree edges such that the forest obtained by their removal induces the given partition (in terms of membership to the same connected component). In this case, the vertices can be uniquely colored by set A of the partition if they lie on a path connecting two elements from A.

Parameters:

Name	Type	Description	Default
T	Tree	A tree with unique leaf labels.	required
partition	list | tuple	A partition of the labels of the trees' leaves.	required
lca	LCA | None	Precomputed LCA datastructure for the tree. The default is None in which case it is computed internally.	None

Returns:

Type	Description
tuple[dict, list] | None	Returns None if the tree and the partition are incompatible, and otherwise a tuple of a dict and a list. The dict contains vertices of the tree that lie on a path connecting two elements from the same set of the partition with the index of this set as value. The list contains the last common ancestor in the tree for each set in the partition.

References

1. M. Hellmuth, D. Schaller, P.F. Stadler. Compatibility of partitions with trees, hierarchies, and split systems. In: Discrete Applied Mathematics 314 (2022) 265-283. doi: https://doi.org/10.1016/j.dam.2022.03.014

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def is_compatible(
    T: Tree,
    partition: list | tuple,
    lca: LCA | None = None,
) -> tuple[dict, list] | None:
    """Checks whether a tree and a partition are compatible.

    A tree with leaf set L and a partition of L are compatible if there is a subset of tree edges
    such that the forest obtained by their removal induces the given partition (in terms of
    membership to the same connected component). In this case, the vertices can be uniquely colored
    by set A of the partition if they lie on a path connecting two elements from A.

    Args:
        T: A tree with unique leaf labels.
        partition: A partition of the labels of the trees' leaves.
        lca: Precomputed LCA datastructure for the tree. The default is None in which case it is
            computed internally.

    Returns:
        Returns None if the tree and the partition are incompatible, and otherwise a tuple of a dict
            and a list. The dict contains vertices of the tree that lie on a path connecting two
            elements from the same set of the partition with the index of this set as value. The
            list contains the last common ancestor in the tree for each set in the partition.

    References:
        1. M. Hellmuth, D. Schaller, P.F. Stadler. Compatibility of partitions with trees,
           hierarchies, and split systems. In: Discrete Applied Mathematics 314 (2022) 265-283.
           doi: https://doi.org/10.1016/j.dam.2022.03.014
    """
    if not isinstance(T, Tree):
        raise TypeError("T must be of type 'Tree'")

    if not isinstance(partition, (list, tuple)):
        raise TypeError("partition must be of type 'list' or 'tuple'")

    if not lca:
        lca = LCA(T)

    # map labels to leaves of the tree
    label_to_leaf = {v.label: v for v in T.leaves()}

    # color each vertex if lies on the path connecting two vertices from the same set
    vertex_coloring = {}

    # last common ancestors of the sets in the partition
    lcas = []

    for i, A in enumerate(partition):
        A_iterator = iter(A)
        x_1 = label_to_leaf[next(A_iterator)]
        current_lca = x_1
        visited = {x_1}
        vertex_coloring[x_1] = i

        for x_label in A_iterator:
            x = label_to_leaf[x_label]
            new_lca = lca(x, current_lca)

            starts = [x]
            if new_lca is not current_lca:
                starts.append(current_lca.parent)
                current_lca = new_lca

            for current in starts:
                while current:
                    if current in visited:
                        break

                    if current in vertex_coloring and vertex_coloring[current] != i:
                        return

                    vertex_coloring[current] = i
                    visited.add(current)

                    if current is new_lca:
                        break

                    current = current.parent

        lcas.append(current_lca)

    return vertex_coloring, lcas

is_refinement_compatible(T: Tree, partition: list | tuple, lca: LCA | None = None) -> tuple[dict, list] | None

Checks whether a tree and a partition are refinement-compatible.

A tree T with leaf set L and a partition of L are refinement-compatible if there is a refinement T' of T for which a subset of edges of T' exists such that the forest obtained by their removal induces the given partition (in terms of membership to the same connected component). In this case, the edges can be uniquely colored by set A of the partition if they lie on a path connecting two elements from A.

Parameters:

Name	Type	Description	Default
T	Tree	A tree with unique leaf labels.	required
partition	list | tuple	A partition of the labels of the trees' leaves.	required
lca	LCA | None	Precomputed LCA datastructure for the tree. The default is None in which case it is computed internally.	None

Returns:

Type	Description
tuple[dict, list] | None	Returns None if the tree and the partition are not refinement-compatible, and otherwise a tuple of a dict and a list. The dict contains nodes v of the tree such that the edge (parent(v), v) lies on a path connecting two elements from the same set of the partition with the index of this set as value. The list contains the last common ancestor in the tree for each set in the partition.

References

1. M. Hellmuth, D. Schaller, P.F. Stadler. Compatibility of partitions with trees, hierarchies, and split systems. In: Discrete Applied Mathematics 314 (2022) 265-283. doi: https://doi.org/10.1016/j.dam.2022.03.014

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def is_refinement_compatible(
    T: Tree,
    partition: list | tuple,
    lca: LCA | None = None,
) -> tuple[dict, list] | None:
    """Checks whether a tree and a partition are refinement-compatible.

    A tree T with leaf set L and a partition of L are refinement-compatible if there is a refinement
    T' of T for which a subset of edges of T' exists such that the forest obtained by their removal
    induces the given partition (in terms of membership to the same connected component). In this
    case, the edges can be uniquely colored by set A of the partition if they lie on a path
    connecting two elements from A.

    Args:
        T: A tree with unique leaf labels.
        partition: A partition of the labels of the trees' leaves.
        lca: Precomputed LCA datastructure for the tree. The default is None in which case it is
            computed internally.

    Returns:
        Returns None if the tree and the partition are not refinement-compatible, and otherwise a
            tuple of a dict and a list. The dict contains nodes v of the tree such that the edge
            (parent(v), v) lies on a path connecting two elements from the same set of the partition
            with the index of this set as value. The list contains the last common ancestor in the
            tree for each set in the partition.

    References:
        1. M. Hellmuth, D. Schaller, P.F. Stadler. Compatibility of partitions with trees,
           hierarchies, and split systems. In: Discrete Applied Mathematics 314 (2022) 265-283.
           doi: https://doi.org/10.1016/j.dam.2022.03.014
    """
    if not isinstance(T, Tree):
        raise TypeError("T must be of type 'Tree'")

    if not isinstance(partition, (list, tuple)):
        raise TypeError("partition must be of type 'list' or 'tuple'")

    if not lca:
        lca = LCA(T)

    # map labels to leaves of the tree
    label_to_leaf = {v.label: v for v in T.leaves()}

    # color each vertex v if the edge (v.parent, v) lies on the path connecting two vertices from
    # the same set
    edge_coloring = {}

    # last common ancestors of the sets in the partition
    lcas = []

    for i, A in enumerate(partition):
        A_iterator = iter(A)
        current_lca = label_to_leaf[next(A_iterator)]
        visited = set()

        for x_label in A_iterator:
            x = label_to_leaf[x_label]
            new_lca = lca(x, current_lca)

            starts = [x]
            if new_lca is not current_lca:
                starts.append(current_lca)
                current_lca = new_lca

            for current in starts:
                while current:
                    if current in edge_coloring and edge_coloring[current] != i:
                        return

                    edge_coloring[current] = i
                    visited.add(current)

                    # parent is always defined since at least the first edge must be added and if
                    # the root is reached, it must equal new_lca
                    current = current.parent
                    if current in visited or current is new_lca:
                        break

        lcas.append(current_lca)

    return edge_coloring, lcas

fitch_orientation(T: Tree, partition: list | tuple, lca: LCA | None = None) -> list[list[str]] | None

(Un)ambiguously present or absent edges in the directed Fitch graph.

The ordered pair (x, y) of two leaves x and y is in the Fitch relation if the path from lca(x, y) to y contains a transfer edge. An undirected Fitch relation ((x, y) if (x, y) or (y, x) in the directed relation) can be represented by a partition (the maximal independent sets). In combination with a tree, this partition determines to some extent the directed Fitch relation.

Parameters:

Name	Type	Description	Default
T	Tree	A tree with unique leaf labels.	required
partition	list | tuple	A partition of the labels of the trees' leaves.	required
lca	LCA | None	Precomputed LCA datastructure for the tree. The default is None in which case it is computed internally.	None

Returns:

Type	Description
list[list[str]] | None	Returns None if the tree and the partition are incompatible. Otherwise, returns a 2-dimensional array that stores at indices i and j corresponding to sets A and B, resp., of the partition whether (A, B) is 'essential', 'forbidden', or 'ambiguous'.

References

1. D. Schaller, M. Hellmuth, P. F. Stadler. Orientation of Fitch Graphs and Reconciliation-Free Inference of Horizontal Gene Transfer in Gene Trees. In: SIAM Journal on Discrete Mathematics, 37(3): 2172-2207 (2023). doi: https://doi.org/10.1137/22M150736X

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def fitch_orientation(
    T: Tree,
    partition: list | tuple,
    lca: LCA | None = None,
) -> list[list[str]] | None:
    """(Un)ambiguously present or absent edges in the directed Fitch graph.

    The ordered pair (x, y) of two leaves x and y is in the Fitch relation if the path from
    lca(x, y) to y contains a transfer edge. An undirected Fitch relation ((x, y) if (x, y) or
    (y, x) in the directed relation) can be represented by a partition (the maximal independent
    sets). In combination with a tree, this partition determines to some extent the directed Fitch
    relation.

    Args:
        T: A tree with unique leaf labels.
        partition: A partition of the labels of the trees' leaves.
        lca: Precomputed LCA datastructure for the tree. The default is None in which case it is
            computed internally.

    Returns:
        Returns None if the tree and the partition are incompatible. Otherwise, returns a
            2-dimensional array that stores at indices i and j corresponding to sets A and B, resp.,
            of the partition whether (A, B) is 'essential', 'forbidden', or 'ambiguous'.

    References:
        1. D. Schaller, M. Hellmuth, P. F. Stadler. Orientation of Fitch Graphs and
           Reconciliation-Free Inference of Horizontal Gene Transfer in Gene Trees. In: SIAM Journal
           on Discrete Mathematics, 37(3): 2172-2207 (2023).
           doi: https://doi.org/10.1137/22M150736X
    """
    if not isinstance(T, Tree):
        raise TypeError("T must be of type 'Tree'")

    if not isinstance(partition, (list, tuple)):
        raise TypeError("partition must be of type 'list' or 'tuple'")

    if not lca:
        lca = LCA(T)

    compatibility = is_compatible(T, partition, lca=lca)
    if not compatibility:
        return

    vertex_coloring, lcas = compatibility

    # compute for each vertex its lowest colored (strict) ancestor
    lowest_colored_ancestor = {}
    for v in T.preorder():
        if not v.parent:
            lowest_colored_ancestor[v] = None
        elif v.parent in vertex_coloring:
            lowest_colored_ancestor[v] = v.parent
        else:
            lowest_colored_ancestor[v] = lowest_colored_ancestor[v.parent]

    n = len(partition)
    matrix = [[None for _ in range(n)] for _ in range(n)]

    for i in range(n):
        for j in range(n):
            if i == j:
                matrix[i][j] = "forbidden"
                continue

            u = lca(lcas[i], lcas[j])  # lca(lca(A), lca(B))
            w = lowest_colored_ancestor[lcas[j]]

            # lca(B) < u and there is a colored vertex w with lca(B) < w <= u
            if lca.ancestor_not_equal(u, lcas[j]) and w and lca(w, u) is u:
                matrix[i][j] = "essential"
            # lca(A) < lca(B)
            elif lca.ancestor_not_equal(lcas[j], lcas[i]):
                matrix[i][j] = "forbidden"
            else:
                matrix[i][j] = "ambiguous"

    return matrix

fitch_orientation_for_refinements(T: Tree, partition: list | tuple, lca: LCA | None = None) -> list[list[str]] | None

(Un)ambiguously present/absent edges in the directed Fitch graph of all refinements.

The ordered pair (x, y) of two leaves x and y is in the Fitch relation if the path from lca(x, y) to y contains a transfer edge. An undirected Fitch relation ((x, y) if (x, y) or (y, x) in the directed relation) can be represented by a partition (the maximal independent sets). In combination with a tree, this partition determines to some extent the directed Fitch relation w.r.t. the possible refinements of the tree.

Parameters:

Name	Type	Description	Default
T	Tree	A tree with unique leaf labels.	required
partition	list | tuple	A partition of the labels of the trees' leaves.	required
lca	LCA | None	Precomputed LCA datastructure for the tree. The default is None in which case it is computed internally.	None

Returns:

Type	Description
list[list[str]] | None	Returns None if the tree and the partition are not refinement-compatible. Otherwise, returns a 2-dimensional array that stores at indices i and j corresponding to sets A and B, resp., of the partition whether (A, B) is 'essential', 'forbidden', or 'ambiguous' for all refinements of the tree that are compatible with the partition.

References

1. D. Schaller, M. Hellmuth, P. F. Stadler. Orientation of Fitch Graphs and Reconciliation-Free Inference of Horizontal Gene Transfer in Gene Trees. In: SIAM Journal on Discrete Mathematics, 37(3): 2172-2207 (2023). doi: https://doi.org/10.1137/22M150736X

Source code in asymmetree/analysis/horizontal_gene_transfer.py

def fitch_orientation_for_refinements(
    T: Tree,
    partition: list | tuple,
    lca: LCA | None = None,
) -> list[list[str]] | None:
    """(Un)ambiguously present/absent edges in the directed Fitch graph of all refinements.

    The ordered pair (x, y) of two leaves x and y is in the Fitch relation if the path from
    lca(x, y) to y contains a transfer edge. An undirected Fitch relation ((x, y) if (x, y) or
    (y, x) in the directed relation) can be represented by a partition (the maximal independent
    sets). In combination with a tree, this partition determines to some extent the directed Fitch
    relation w.r.t. the possible refinements of the tree.

    Args:
        T: A tree with unique leaf labels.
        partition: A partition of the labels of the trees' leaves.
        lca: Precomputed LCA datastructure for the tree. The default is None in which case it is
            computed internally.

    Returns:
        Returns None if the tree and the partition are not refinement-compatible. Otherwise, returns
            a 2-dimensional array that stores at indices i and j corresponding to sets A and B,
            resp., of the partition whether (A, B) is 'essential', 'forbidden', or 'ambiguous' for
            all refinements of the tree that are compatible with the partition.

    References:
        1. D. Schaller, M. Hellmuth, P. F. Stadler. Orientation of Fitch Graphs and
           Reconciliation-Free Inference of Horizontal Gene Transfer in Gene Trees. In: SIAM Journal
           on Discrete Mathematics, 37(3): 2172-2207 (2023).
           doi: https://doi.org/10.1137/22M150736X
    """
    if not isinstance(T, Tree):
        raise TypeError("T must be of type 'Tree'")

    if not isinstance(partition, (list, tuple)):
        raise TypeError("partition must be of type 'list' or 'tuple'")

    if not lca:
        lca = LCA(T)

    r_compatibility = is_refinement_compatible(T, partition, lca=lca)
    if not r_compatibility:
        return

    edge_coloring, lcas = r_compatibility

    # compute for each vertex its lowest colored (strict) ancestor
    lowest_colored_edge = {}
    for v in T.preorder():
        if not v.parent:
            lowest_colored_edge[v] = None
        elif v in edge_coloring:
            lowest_colored_edge[v] = v
        else:
            lowest_colored_edge[v] = lowest_colored_edge[v.parent]

    n = len(partition)
    matrix = [[None for _ in range(n)] for _ in range(n)]

    for i in range(n):
        for j in range(n):
            if i == j:
                matrix[i][j] = "forbidden"
                continue

            u = lca(lcas[i], lcas[j])  # lca(lca(A), lca(B))
            v_A = lcas[i]  # lca(A)
            v_B = lcas[j]  # lca(B)
            lce_B = lowest_colored_edge[v_B]

            if lca.ancestor_not_equal(u, v_B) and lce_B and lca.ancestor_not_equal(u, lce_B):
                matrix[i][j] = "essential"
            elif (
                u.parent
                and lca.ancestor_not_equal(u, v_A)
                and u in edge_coloring
                and edge_coloring[u] == edge_coloring.get(_find_child(u, v_A, lca))
            ):
                matrix[i][j] = "essential"
            elif (
                lca.ancestor_not_equal(v_B, v_A)
                and edge_coloring.get(_find_child(u, v_A, lca)) == j
            ):
                matrix[i][j] = "forbidden"
            else:
                matrix[i][j] = "ambiguous"

    return matrix