python - 如何为给定外群的一组物种生成所有可能的 Newick 树排列？

Question

如何为给定外群的一组物种生成所有可能的 Newick 树排列？

对于那些不知道 Newick 树格式是什么的人，可以在以下位置找到很好的描述: https://en.wikipedia.org/wiki/Newick_format

我想为给定外群的一组物种创建所有可能的 Newick 树排列。我期望处理的叶节点数最有可能是 4、5 或 6 个叶节点。

允许“软”和“硬”多项式。 https://en.wikipedia.org/wiki/Polytomy#Soft_polytomies_vs._hard_polytomies https://biology.stackexchange.com/questions/23667/evidence-discussions-of-hard-polytomy

下图是理想输出，“E”设为外群

理想输出:

((("A","B","C"),("D"),("E"));
((("A","B","D"),("C"),("E"));
((("A","C","D"),("B"),("E"));
((("B","C","D"),("A"),("E"));
((("A","B")("C","D"),("E"));
((("A","C")("B","D"),("E"));
((("B","C")("A","D"),("E"));
(("A","B","C","D"),("E"));
(((("A","B"),"C"),"D"),("E"));

但是，我使用 itertools 的任何可能的解决方案，特别是 itertools.permutations，都遇到了等效输出的问题。我想到的最后一个想法涉及如下所示的等效输出。

等效输出:

((("C","B","A"),("D"),("E"));
((("C","A","B"),("D"),("E"));
((("A","C","B"),("D"),("E"));

这是我的解决方案想法的开始。但是，我现在不太确定除了 itertools 之外还能解决这个问题。

import itertools

def Newick_Permutation_Generator(list_of_species, name_of_outgroup)
    permutations_list =(list(itertools.permutations(["A","B","C","D","E"])))

    for given_permutation in permutations_list:
        process(given_permutation)

Newick_Permutation_Generator(["A","B","C","D","E"], "E")

Answer 1

一棵树作为叶子集合的递归集

让我们暂时搁置 newick 表示，并考虑问题的可能 python 表示。

有根树可以看作是集合的递归层次结构(集合(集合...))叶。集合是无序的，非常适合描述树中的进化枝:{{{"A", "B"}, {"C", "D"}}, "E"} 应该与 {{{"C", "D"}, {"B", "A"}}, "E"} 相同。

如果我们考虑初始叶集 {"A", "B", "C", "D", "E"}，则以“E”为外群的树是{tree, "E"} 形式的一组集合，其中 tree 取自可以从叶子集合构建的树集合 { “A”、“B”、“C”、“D”。我们可以尝试编写一个递归的 trees 函数来生成这组树，我们的总树集将表示如下:

{{tree, "E"} for tree in trees({"A", "B", "C", "D"})}

(在这里，我使用了 set comprehension 符号。)

实际上，python 不允许集合的集合，因为集合的元素必须是“可散列的”(也就是说，python 必须能够计算对象的一些“散列”值才能检查它们是否属于或不到集合)。恰好python集没有这个属性。幸运的是，我们可以使用名为 frozenset 的类似数据结构。 ，它的行为很像一个集合，但不能修改并且是“可散列的”。因此，我们的树集将是:

all_trees = frozenset(
    {frozenset({tree, "E"}) for tree in trees({"A", "B", "C", "D"})})

实现树函数

现在让我们关注 trees 函数。

对于叶子集合的每个可能的分区(分解成一组不相交的子集，包括所有元素)，我们需要为每个部分找到所有可能的树(通过递归调用)分区。对于给定的分区，我们将为每个可能的子树组合创建一棵树。

例如，如果一个分区是 {"A", {"B", "C", "D"}}，我们将考虑所有可能的树，可以从 "A"(实际上，只是叶子 "A" 本身)，以及可以从部分 {"B", "C", “D”(即 trees({"B", "C", "D"}))。然后，该分区的可能树将通过获取所有可能的对来获得，其中一个元素仅来自 "A"，而另一个来自 trees({"B", "C", "D"}).

这可以推广到具有两个以上部分的分区，itertools 中的 product 函数在这里似乎很有用。

因此，我们需要一种方法来生成一组叶子的可能分区。

生成集合的分区

这里我做了一个 partitions_of_set 函数改编自 this solution :

# According to https://stackoverflow.com/a/30134039/1878788:
# The problem is solved recursively:
# If you already have a partition of n-1 elements, how do you use it to partition n elements?
# Either place the n'th element in one of the existing subsets, or add it as a new, singleton subset.
def partitions_of_set(s):
    if len(s) == 1:
        yield frozenset(s)
        return
    # Extract one element from the set
    # https://stackoverflow.com/a/43804050/1878788
    elem, *_ = s
    rest = frozenset(s - {elem})
    for partition in partitions_of_set(rest):
        for subset in partition:
            # Insert the element in the subset
            try:
                augmented_subset = frozenset(subset | frozenset({elem}))
            except TypeError:
                # subset is actually an atomic element
                augmented_subset = frozenset({subset} | frozenset({elem}))
            yield frozenset({augmented_subset}) | (partition - {subset})
        # Case with the element in its own extra subset
        yield frozenset({elem}) | partition

为了检查获得的分区，我们创建了一个函数来使它们更容易显示(这对于稍后制作树的 newick 表示也很有用):

def print_set(f):
    if type(f) not in (set, frozenset):
        return str(f)
    return "(" + ",".join(sorted(map(print_set, f))) + ")"

我们测试分区是否有效:

for partition in partitions_of_set({"A", "B", "C", "D"}):
    print(len(partition), print_set(partition))

输出:

1 ((A,B,C,D))
2 ((A,B,D),C)
2 ((A,C),(B,D))
2 ((B,C,D),A)
3 ((B,D),A,C)
2 ((A,B,C),D)
2 ((A,B),(C,D))
3 ((A,B),C,D)
2 ((A,D),(B,C))
2 ((A,C,D),B)
3 ((A,D),B,C)
3 ((A,C),B,D)
3 ((B,C),A,D)
3 ((C,D),A,B)
4 (A,B,C,D)

树函数的实际代码

现在我们可以编写树函数了:

from itertools import product
def trees(leaves):
    if type(leaves) not in (set, frozenset):
        # It actually is a single leaf
        yield leaves
        # Don't try to yield any more trees
        return
    # Otherwise, we will have to consider all the possible
    # partitions of the set of leaves, and for each partition,
    # construct the possible trees for each part
    for partition in partitions_of_set(leaves):
        # We need to skip the case where the partition
        # has only one subset (the initial set itself),
        # otherwise we will try to build an infinite
        # succession of nodes with just one subtree
        if len(partition) == 1:
            part, *_ = partition
            # Just to be sure the assumption is correct
            assert part == leaves
            continue
        # We recursively apply *tree* to each part
        # and obtain the possible trees by making
        # the product of the sets of possible subtrees.
        for subtree in product(*map(trees, partition)):
            # Using a frozenset guarantees
            # that there will be no duplicates
            yield frozenset(subtree)

测试它:

all_trees = frozenset(
    {frozenset({tree, "E"}) for tree in trees({"A", "B", "C", "D"})})

for tree in all_trees:
    print(print_set(tree) + ";")

输出:

(((B,C),A,D),E);
((((A,B),D),C),E);
((((B,D),A),C),E);
((((C,D),A),B),E);
(((A,D),B,C),E);
((A,B,C,D),E);
((((B,D),C),A),E);
(((A,B,C),D),E);
((((A,C),B),D),E);
((((C,D),B),A),E);
((((B,C),A),D),E);
(((A,B),C,D),E);
(((A,C),(B,D)),E);
(((B,D),A,C),E);
(((C,D),A,B),E);
((((A,B),C),D),E);
((((A,C),D),B),E);
(((A,C,D),B),E);
(((A,D),(B,C)),E);
((((A,D),C),B),E);
((((B,C),D),A),E);
(((A,B),(C,D)),E);
(((A,B,D),C),E);
((((A,D),B),C),E);
(((A,C),B,D),E);
(((B,C,D),A),E);

希望结果是正确的

要正确使用这种方法有点棘手。我花了一些时间才弄清楚如何避免无限递归(当分区为 {{"A", "B", "C", "D"}} 时会发生这种情况)。