如果给出一个外群,我如何为一组物种生成所有可能的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"设置为outgroup

理想输出:

((("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"作为outgroup的树是表格{tree, "E"}中tree的集合集合,其中s取自可以从叶子集合构建的树集合{"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功能

现在让我们关注这个trees功能.

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

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

这可以推广到具有两个以上部分的分区,并且product函数from itertools在这里似乎很有用.

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

生成集合的分区

在这里,我做了一个partitions_of_set改编自功能该解决方案:

# 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)

trees功能的实际代码

现在我们可以编写tree函数:

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"}}).