【ROSALIND】【练Python，学生信】48 Newick格式与进化树

如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面  谢谢配合~

题目：

Newick格式与进化树结点间的距离（Distances in Trees）

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

        尽管已经找到现成的代码，我还是在运行时遇到了问题，在此记录一下以防大家重复踩雷。情况是这样的，我运行Python时使用的IDE是PyCharm，安装新模块时就直接在File-Settings中下载安装。然而这次安装成功后，我用import Bio却显示错误，似乎这个包没有安装上。经过搜索后，我终于在下面这个网址：https://www.cnpython.com/qa/72896找到了解决方案，3楼的网友与我遇到的问题相同，按照他的方法把安装文件夹路径修改，并从“bio”重命名为“Bio”后，代码正常运行并获得了正确结果。

 

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 = line.split() # 把两个叶子结点分开保存到x和y中
    tree = Phylo.read(io.StringIO(i),'newick') # 用Biopython中Phylo模块方法直接把Newick解析成树结构
    clades = tree.find_clades()
    for clade in clades: # 逐条给各分支加上长度，这里统一赋值为1
        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=' ')