Newick格式与进化树结点间的距离(Distances in

2022-10-08  本文已影响0人  iBioinformatics

Given: A collection of n trees (n≤40) in Newick format, with each tree containing at most 200 nodes; each tree Tk is followed by a pair of nodes xk and yk in Tk.

所给:n个(不超过40个)以Newick格式表示的进化树,每棵树包含的结点不超过200个;所给文件的形式为一棵树Tk,后面跟随着这棵树中的一对结点xk和yk。

Return: A collection of n positive integers, for which the kth integer represents the distance between xk and yk in Tk.

需得:n个正整数,分别代表xk和yk在Tk内的距离。

测试数据

(cat)dog;

dog cat

(dog,cat);

dog cat

测试输出

1 2

生物学背景

进化树是在生物学中用来表示物种之间进化关系的图。进化树上每个叶子结点代表一个物种,叶子结点之间的距离可以代表两个物种之间的差异程度。

Newick是最常见的进化树文件格式。举例来说,假如n个叶子结点{v1, v2, … ,vn}都和一个内结点u相连,如下图:

用Newick格式可表示为(v1,v2,…,v3)u。Newick格式中结点名称并非必须,我们可以省略一些或全部结点名称,将这棵树表示为(v1,v2,…,v3),(, ,…,)u或者(, ,…,)。

无根树来说在用Newick表示时可以把任意结点作为根节点,因此可以用多个Newick描述同一棵树,如(C, D, (A, B))和(A, (D, C), B)都表示下图:

对Newick格式更详细地解释可以参考这篇文章:https://blog.csdn.net/weixin_43569478/article/details/108079223

数学背景

对于一棵树中的两个结点,肯定能找到唯一路径将它们相连。因为如果两点间有多条路径,必然会形成循环,这与树的定义不相符。

思路

方法1:善于利用现成的工具包

本方法的代码来自于如下博文:http://saradoesbioinformatics.blogspot.com/2016/09/distances-in-trees.html

在这里博主直接用了Biopython包Phylo模块中的现成方法,这样就无需自己动手对Newick格式进行解析,Phylo模块的介绍如下:https://biopython.org/wiki/Phylo

import sys
from Bio import Phylo # 引入Biopython包中的Phylo模块
import io

f = open('rosalind_nwck.txt','r')
# 将所给的输入序列读入并按条分开
pairs = [i.split('\n') for i in f.read().strip().split('\n\n')]  

for i, line in pairs:
    # 把两个叶子结点分开保存到x和y中
    x,y = line.split() 
    # 用Biopython中Phylo模块方法直接把Newick解析成树结构
    tree = Phylo.read(io.StringIO(i),'newick') 
    clades = tree.find_clades()
    # 逐条给各分支加上长度,这里统一赋值为1
    for clade in clades: 
        clade.branch_length = 1
    sys.stdout.write('%s' % tree.distance(x,y) + ' ') # 直接输出两结点间的距离,sys.stdout.write()输出不会自动换行
sys.stdout.write('\n') 

方法2:理解括号和逗号的含义

借助方法1的代码成功混入讨论区后,我终于得以阅读大佬们的解题思路,下面的代码是我在学习后给出的,但有一些地方需要我们静下心来理解。

首先,要完成这道题并不需要写一个方法把Newick格式解析为树形结构,再数出两个结点的距离。我们需要把’(‘ , ’)’和’,’这些括号与树结构建立起联系,这样只需要关注这些符号的特点就可以计算出结点距离。

从网站给出的例子,我们能对括号和逗号有一个直觉上的认识,即:逗号代表并列,括号代表层级。比如:(A,B),A和B的关系应该如下图左;(A)B应该如下图右。这似乎意味着两个结点间多半个括号,距离应加1,多一个逗号,距离应加上2.

不过,实际情况并非这么简单。我们看下面这个示例:(A,( (()),()((,(,))) ),B),A、B之间似乎远隔千山万水,但实际它们的距离只有2,原因如下图:

可以看到,这里A、B中间的内容都是配对的括号,形成一个完整的分支,与A、B同为一个父结点的不同分支,我们在计算A、B距离时尽可以把这个分支整个忽略掉。在计算A、B间的距离时,我们只需要关注A、B间不配对的括号即可,成功配对的括号全部消去,每多一个不匹配的括号,距离就加上1。

对于逗号来说,我把它看成’)(‘,直接用’)(‘替换掉原文中的’,’,以上文的方法消去配对的括号,只剩下错配的括号,就是最终的结果。如何理解呢?我们画几幅图来说明。

假设一个进化树中没有逗号,在删去其中的配对括号后,变成了A))))B或A((((B这种形式,可以用图表示为:

显然,A、B间的距离为4。

现在我们让情况复杂些,最后一个错配括号后面出现了一个逗号,如A)))),B,现在图可以这样表示:

可以看到,加了一个逗号之后,距离增加了2。那为什么不直接计算多出逗号的数量,然后乘2加到最终结果里呢?我们再看一个例子,这次最后一个错配括号前也有逗号,如A),))),B,用图可以表示为:

可以看到,这种情况A、B间的距离依然是6。事实上,’,’的含义相当于一个右括号和一个左括号。

import re

def nwck(data, target):
        newdata = []
        for num, i in enumerate(re.split(r'([(),])', data)): # 用re包处理字符串,以'(',')'和','将树分割
            if i != '' and i != ';\n': # 除去空字符
                newdata.append(i)
        # print(newdata)

        num1 = num2 = -1
        for num, i in enumerate(newdata): # 将两结点的位置赋给num1和num2
            if i == target[0]:
                num1 = num
            if i == target[1].replace('\n', ''):
                num2 = num
            if num1 != -1 and num2 != -1: # 两结点都找到就可以停止了
                break

        start = min(num1, num2) # 左节点赋值给start
        end = max(num1, num2) # 右节点赋值给end

        left = right = 0
        for i in range(start + 1, end): # 从左节点像右节点遍历
            # print(newdata[i])
            if newdata[i] in ',)': # 碰到右括号含义的字符
                if left > 0: # 如果左边有与其配对的左括号
                    left -= 1 # 删去左括号
                else: # 如果左边没有与其配对的左括号
                    right += 1 # 记录这个错配右括号
            if newdata[i] in ',(': # 碰到左括号含义的字符
                left += 1 # 记录左括号

        result = left + right
        return result


f = open('rosalind_nwck.txt','r')
input = f.readlines()
f.close()

trees = [] # 用trees保存每棵树
pairs = [] # 用pairs保存各对结点
for num, i in enumerate(input): # 根据所给输入数据的格式
    if num % 3 == 0: # 第3k行是进化树
        trees.append(i)
    if num % 3 == 1: # 第3k+1行是结点,第3k+2行是空行
        pairs.append(i.split(' '))

result = [] # 用result存储距离
for i in range(len(trees)):
    result.append(nwck(trees[i], pairs[i])) # 调用newick函数解决问题

for i in result:
    print(i, end=' ')

转自:https://www.bilibili.com/read/cv8524280/

上一篇下一篇

猜你喜欢

热点阅读