如何为给定外群的一组物种生成所有可能的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.排列，遇到了等效输出的问题。我想到的最后一个想法涉及到如下所示的等效输出。在

等效输出：

^{pr2}$

这是我的解决方案的开始。不过，除了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

这是个很难回答的问题！这是我的旅程。在

第一个观察是outgroup始终是一个连接到newick字符串末尾的节点。让我们把其余的物种称为内群，并尝试生成这些物种的所有排列。然后简单地添加outgroup。在

from itertools import permutations

def ingroup_generator(species, n):
    for perm in permutations(species, n):
        yield tuple([tuple(perm), tuple(s for s in species if s not in perm)])

def format_newick(s, outgroup=''):
    return '(' + ', '.join('({})'.format(', '.join(p)) for p in s) + ',({}));'.format(outgroup)

species = ["A","B","C","D","E"]
outgroup = "E"
ingroup = [s for s in species if s != outgroup]

itertools_newicks= []
for n in range(1, len(ingroup)):
    for p in ingroup_generator(ingroup, n):
        itertools_newicks.append(format_newick(p, outgroup))

for newick in itertools_newicks:
    print newick

这将返回40个newick字符串：

^{pr2}$

其中有些是重复的，但我们稍后将删除这些重复项。在

作为bli noted in the comments，(((("A","B"),"C"),"D"),("E"));及其变体也应被视为有效解。 comments on BioStar向我指出了正确的方向，这与生成二叉树的所有可能的分组是一样的。我找到了一个漂亮的Python implementation in this StackOverflow answer by rici：

# A very simple representation for Nodes. Leaves are anything which is not a Node.
class Node(object):
  def __init__(self, left, right):
    self.left = left
    self.right = right

  def __repr__(self):
    return '(%s, %s)' % (self.left, self.right)

# Given a tree and a label, yields every possible augmentation of the tree by
# adding a new node with the label as a child "above" some existing Node or Leaf.
def add_leaf(tree, label):
  yield Node(label, tree)
  if isinstance(tree, Node):
    for left in add_leaf(tree.left, label):
      yield Node(left, tree.right)
    for right in add_leaf(tree.right, label):
      yield Node(tree.left, right)

# Given a list of labels, yield each rooted, unordered full binary tree with
# the specified labels.
def enum_unordered(labels):
  if len(labels) == 1:
    yield labels[0]
  else:
    for tree in enum_unordered(labels[1:]):
      for new_tree in add_leaf(tree, labels[0]):
        yield new_tree

然后

enum_newicks= []
for t in enum_unordered(ingroup):
    enum_newicks.append('({},({}));'.format(t, outgroup))

for newick in enum_newicks:
    print newick

生成以下15个新图标：

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

所以现在我们已经有了40+15=55个可能的newick字符串，我们必须删除重复项。在

我尝试的第一个死胡同是为每个newick字符串创建一个规范化的表示，这样我就可以将它们用作字典中的键。其想法是递归地对所有节点中的字符串进行排序。但首先我必须捕获所有（嵌套）节点。我不能使用正则表达式，因为nested structures are by definition not regular。在

因此，我使用了pyparsing包，并提出了这个：

from pyparsing import nestedExpr 

def sort_newick(t):
    if isinstance(t, str):
        return sorted(t)
    else:
        if all(isinstance(c, str) for c in t):
            return sorted(t)
        if all(isinstance(l, list) for l in t):
            return [sort_newick(l) for l in sorted(t, key=lambda k: sorted(k))]
        else:
            return [sort_newick(l) for l in t]


def canonical_newick(n):
    n = n.replace(',', '')
    p = nestedExpr().parseString(n).asList()
    s = sort_newick(p)
    return str(s)

这让我屈服了

from collections import defaultdict

all_newicks = itertools_newicks + enum_newicks

d = defaultdict(list)
for newick in all_newicks:
    d[canonical_newick(newick)].append(newick)

for canonical, newicks in d.items():
    print canonical
    for newick in newicks:
        print '\t', newick
    print

有22个键的字典：

[[[['A'], [['C'], ['B', 'D']]], ['E']]]
    ((A, (C, (B, D))),(E));

[[[['B'], [['A'], ['C', 'D']]], ['E']]]
    ((B, (A, (C, D))),(E));

[[[['B'], [['A', 'C'], ['D']]], ['E']]]
    ((B, ((A, C), D)),(E));

[[['A', 'C', 'D'], ['B'], ['E']]]
    ((B), (A, C, D),(E));
    ((A, C, D), (B),(E));
    ((A, D, C), (B),(E));
    ((C, A, D), (B),(E));
    ((C, D, A), (B),(E));
    ((D, A, C), (B),(E));
    ((D, C, A), (B),(E));

[[['A', 'B'], ['C', 'D'], ['E']]]
    ((A, B), (C, D),(E));
    ((B, A), (C, D),(E));
    ((C, D), (A, B),(E));
    ((D, C), (A, B),(E));

[[[[['A'], ['B', 'C']], ['D']], ['E']]]
    (((A, (B, C)), D),(E));

[[[['A', 'C'], ['B', 'D']], ['E']]]
    (((A, C), (B, D)),(E));

[[['A'], ['B', 'C', 'D'], ['E']]]
    ((A), (B, C, D),(E));
    ((B, C, D), (A),(E));
    ((B, D, C), (A),(E));
    ((C, B, D), (A),(E));
    ((C, D, B), (A),(E));
    ((D, B, C), (A),(E));
    ((D, C, B), (A),(E));

[[[['A', 'D'], ['B', 'C']], ['E']]]
    (((B, C), (A, D)),(E));

[[['A', 'B', 'C'], ['D'], ['E']]]
    ((D), (A, B, C),(E));
    ((A, B, C), (D),(E));
    ((A, C, B), (D),(E));
    ((B, A, C), (D),(E));
    ((B, C, A), (D),(E));
    ((C, A, B), (D),(E));
    ((C, B, A), (D),(E));

[[['A', 'C'], ['B', 'D'], ['E']]]
    ((A, C), (B, D),(E));
    ((B, D), (A, C),(E));
    ((C, A), (B, D),(E));
    ((D, B), (A, C),(E));

[[['A', 'B', 'D'], ['C'], ['E']]]
    ((C), (A, B, D),(E));
    ((A, B, D), (C),(E));
    ((A, D, B), (C),(E));
    ((B, A, D), (C),(E));
    ((B, D, A), (C),(E));
    ((D, A, B), (C),(E));
    ((D, B, A), (C),(E));

[[[['A'], [['B'], ['C', 'D']]], ['E']]]
    ((A, (B, (C, D))),(E));

[[[[['A', 'B'], ['C']], ['D']], ['E']]]
    ((((A, B), C), D),(E));

[[[[['B'], ['A', 'C']], ['D']], ['E']]]
    (((B, (A, C)), D),(E));

[[[['C'], [['B'], ['A', 'D']]], ['E']]]
    ((C, (B, (A, D))),(E));

[[[['C'], [['A', 'B'], ['D']]], ['E']]]
    ((C, ((A, B), D)),(E));

[[[['A'], [['B', 'C'], ['D']]], ['E']]]
    ((A, ((B, C), D)),(E));

[[[['A', 'B'], ['C', 'D']], ['E']]]
    (((A, B), (C, D)),(E));

[[[['B'], [['C'], ['A', 'D']]], ['E']]]
    ((B, (C, (A, D))),(E));

[[[['C'], [['A'], ['B', 'D']]], ['E']]]
    ((C, (A, (B, D))),(E));

[[['A', 'D'], ['B', 'C'], ['E']]]
    ((A, D), (B, C),(E));
    ((B, C), (A, D),(E));
    ((C, B), (A, D),(E));
    ((D, A), (B, C),(E));

但仔细检查发现了一些问题。让我们看看newick '(((A, B), (C, D)),(E));和((D, C), (A, B),(E));的例子。在我们的字典d中，它们具有不同的规范密钥，分别是[[[['A', 'B'], ['C', 'D']], ['E']]]和{}。但事实上，这些是重复的树。我们可以通过观察两棵树之间的Robinson-Foulds distance来证实这一点。如果它是零，树是相同的。在

我们使用ete3 toolkit package中的robinson_foulds函数

from ete3 import Tree

tree1 = Tree('(((A, B), (C, D)),(E));')
tree2 = Tree('((D, C), (A, B),(E));')

rf, max_parts, common_attrs, edges1, edges2, discard_t1, discard_t2 = tree1.robinson_foulds(tree2, unrooted_trees=True)
    print rf # returns 0

好的，那么Robinson Foulds是一个更好的方法来检查newick树的等式，而不是我的标准树方法。让我们将所有newick字符串包装在一个自定义的MyTree对象中，其中相等被定义为Robinson-Foulds距离为零：

class MyTree(Tree):

    def __init__(self, *args, **kwargs):
        super(MyTree, self).__init__(*args, **kwargs)

    def __eq__(self, other):
        rf = self.robinson_foulds(other, unrooted_trees=True)
        return not bool(rf[0])

trees = [MyTree(newick) for newick in all_newicks]

如果我们还可以定义一个__hash__()函数，为重复树返回相同的值，那么set(trees)将自动删除所有重复项。在

不幸的是，I haven't been able to find a good way to define ^{}，但是有了__eq__，我可以make use of ^{}：

unique_trees = [trees[i] for i in range(len(trees)) if i == trees.index(trees[i])]
unique_newicks = [tree.write(format=9) for tree in unique_trees]
for unique_newick in unique_newicks:
    print unique_newick

所以，我们的旅程到此结束了。我不能完全证明这是正确的解决方案，但我很有信心，以下19个newick都是可能的不同排列：

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

如果我们将每个newick与所有其他newick成对比较，就可以确认这个列表中不再有重复项了

from itertools import product

for n1, n2 in product(unique_newicks, repeat=2):
    if n1 != n2:
        mt1 = MyTree(n1)
        mt2 = MyTree(n2)
        assert mt1 != mt2