首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >Dendropy:在两个节点中间添加内部节点

Dendropy:在两个节点中间添加内部节点
EN

Stack Overflow用户
提问于 2022-05-19 15:52:25
回答 1查看 57关注 0票数 1

我对DendroPy非常陌生。我想要做的事情似乎很简单,但我不知道如何正确地做,我没有在互联网上找到任何东西。

我想在现有的树状树状树结构中,在两个节点之间添加一个节点。

代码语言:javascript
复制
from dendropy import Tree, Taxon, Node

t1 = Tree.get_from_string(
    "(((sp1: 0.35, sp2: 0.15):0.75, sp3:1): 0.5, (sp4: 0.5, sp5: 0.05)MRCA_sp4&sp5: 1)root;",
    "newick", rooting='force-rooted'
)
t1.print_plot()
mrca = t1.mrca(taxon_labels=["sp4", "sp5"])
print(mrca.description())

树和MRCA节点描述

正确地找到了sp4和sp5的MRCA。现在,我尝试在MRCA和root之间添加一个节点,使用下面的代码:

代码语言:javascript
复制
def add_node_midway_between_2_nodes(lowernode, taxon_label=None, node_label=None):
    newtaxon = Taxon(label=taxon_label)
    newnode = Node(taxon=newtaxon, label=node_label)
    newnode.parent_node = lowernode.parent_node
    newnode.edge_length = lowernode.edge_length/2
    lowernode.parent_node = newnode
    lowernode.edge_length = newnode.edge_length
    return newnode

node = add_node_midway_between_2_nodes(mrca, node_label="midway between root and MRCA sp4&sp5")
t1.print_plot()
str_t1 = t1.as_string(schema='newick')
print(str_t1)

根与MRCASP4&SP5之间有节点的树

&R根;

看了看情节和字符串,它似乎奏效了。但是,当我再次尝试计算sp4 et sp5的MRCA时,它不再找到“MRCASP4&SP5”,而是找到根节点。

代码语言:javascript
复制
mrca = t1.mrca(taxon_labels=["sp4", "sp5"])
print(mrca.description())

Output =根节点的描述

从parent_node到sp5,我仍然可以找到“MRCASP4&SP5”。

出于绝望,我尝试使用字符串str_t1重做树,但也没有工作,甚至给出了另一个结果(同样不正确):节点“根和MRCA sp4&sp5之间的中间”。

代码语言:javascript
复制
t1 = Tree.get_from_string(
    str_t1,
    "newick", rooting='force-rooted'
)
mrca = t1.mrca(taxon_labels=["sp4", "sp5"])
print(mrca.description())

Output =节点“根到MRCASP4和SP5之间的中途”的描述

那么,在两个节点中间添加一个节点的干净方法是什么,这不会在之后创建奇怪的事件?

非常感谢

EN

回答 1

Stack Overflow用户

发布于 2022-09-18 17:36:05

您的代码主要是工作的,但是您应该按照update_bipartitions中的建议,正确地应用树拓扑的任何更改。所以,在你的例子中,它看起来是这样的:

代码语言:javascript
复制
def add_node_midway_between_2_nodes(lowernode, taxon_label=None, node_label=None):
    newtaxon = Taxon(label=taxon_label)
    newnode = Node(taxon=newtaxon, label=node_label)
    newnode.parent_node = lowernode.parent_node
    newnode.edge_length = lowernode.edge_length/2
    lowernode.parent_node = newnode
    lowernode.edge_length = newnode.edge_length
    return newnode

node = add_node_midway_between_2_nodes(
    mrca, node_label="midway between root and MRCA sp4&sp5"
)
t1.update_taxon_namespace()
t1.update_bipartitions(
    suppress_unifurcations=False, suppress_storage=True
)  # suppress_storage is optional, I just do not want to create a bipartitions list
t1.print_plot()
str_t1 = t1.as_string(schema='newick')
print(str_t1)

NB!

更新taxa命名空间应该在更新bu分区之前,因为后者必须使用正确的TaxonNamespace。否则,你还是会有奇怪的行为。

然而,最好采用内建Node方法进行精细树的重建。例如,我会以这样的方式重写函数:

代码语言:javascript
复制
def insert_new_node_posterior(
    node: Node,
    *,
    taxon_label: Optional[str] = None,
    node_label: Optional[str] = None,
    edge_length: Real,
    # If it was product or at least reusable in the future code,
    # I would add more arguments for proportion specification,
    # using height, distance from root &c.
) -> Node:
    parent = node.parent_node
    if not parent:
        raise Exception("You cannot insert a node in posterior to the root.")

    new_taxon = Taxon(label=taxon_label)
    i = parent.child_nodes().index(node)
    parent.remove_child(node)
    intermediate_node = parent.insert_new_child(
        index=i, taxon=new_taxon, label=node_label, edge_length=edge_length
    )
    node.edge_length -= edge_length
    intermediate_node.add_child(node)
    return intermediate_node


node = insert_new_node_posterior(
    mrca,
    node_label="midway between root and MRCA sp4&sp5",
    edge_length=mrca.edge_length / 2
)
t1.update_taxon_namespace()
t1.update_bipartitions(
    suppress_unifurcations=False, suppress_storage=True
)

t1.print_plot()
str_t1 = t1.as_string(schema='newick')
print(str_t1)

尽管如此,Tree.mrca仍然显示了一个不正确的节点:

代码语言:javascript
复制
Node object at 0x1be8b959460<Node object at 0x1be8b959460: 'midway between root and MRCA sp4&sp5' (<Taxon 0x1be8cdceeb0 'None'>)>
    [Edge]
        Edge object at 0x1be8b959400 (1917897249792, Length=0.5)
    [Taxon]
         Taxon object at 0x1be8cdceeb0: <Unnamed Taxon>
    [Parent]
        Node object at 0x1be8c3b2b80<Node object at 0x1be8c3b2b80: 'root' (None)>
    [Children]
        [0] Node object at 0x1be8be4d6a0<Node object at 0x1be8be4d6a0: 'MRCA sp4&sp5' (None)>

不过,在这种情况下,这并不是一个bug,因为这只是方法的一个特性。由于这一点在源代码

代码语言:javascript
复制
if cms:
    # for at least one taxon cm has 1 and bipartition has 1
    if cms == leafset_bitmask:
        # curr_node has all of the 1's that bipartition has
        if cm == leafset_bitmask:
            return curr_node  # Vovin's comment: Since there is a unifurcation,
                              # it returns the current node
                              # instead of the next iteration
        last_match = curr_node
        nd_source = iter(curr_node.child_nodes())
    else:
        # we have reached a child that has some, but not all of the
        #   required taxa as descendants, so we return the last_match
        return last_match

例如,如果我们将一个子节点添加到新节点中,那么它工作得很好:

代码语言:javascript
复制
node = insert_new_node_posterior(
    mrca,
    node_label="midway between root and MRCA sp4&sp5",
    edge_length=mrca.edge_length / 2
)
node.new_child(label="sp6", edge_length=1, taxon=Taxon(label="sp6"))
t1.update_taxon_namespace()
t1.update_bipartitions(
    suppress_unifurcations=False, suppress_storage=True
)
mrca = t1.mrca(taxon_labels=["sp4", "sp5"])
print(mrca.description())
代码语言:javascript
复制
Node object at 0x1be8c1e22b0<Node object at 0x1be8c1e22b0: 'MRCA sp4&sp5' (None)>
    [Edge]
        Edge object at 0x1be8c1e2340 (1917906199360, Length=0.5)
    [Taxon]
        None
    [Parent]
        Node object at 0x1be8ba8da30<Node object at 0x1be8ba8da30: 'midway between root and MRCA sp4&sp5' (<Taxon 0x1be8ba8d370 'None'>)>
    [Children]
        [0] Node object at 0x1be8c1e2910<Node object at 0x1be8c1e2910: 'None' (<Taxon 0x1be8c1e2640 'sp4'>)>
        [1] Node object at 0x1be8c1e28e0<Node object at 0x1be8c1e28e0: 'None' (<Taxon 0x1be8c1e2e50 'sp5'>)>
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/72307591

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档