如果第一次閱讀本系列文檔請先移步閱讀【ROSALIND】【練Python，學(xué)生信】00 寫在前面 ?謝謝配合~

題目：

Newick格式與進(jìn)化樹結(jié)點(diǎn)間的距離（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個(gè)（不超過40個(gè)）以Newick格式表示的進(jìn)化樹，每棵樹包含的結(jié)點(diǎn)不超過200個(gè)；所給文件的形式為一棵樹Tk，后面跟隨著這棵樹中的一對結(jié)點(diǎn)xk和yk。

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

需得：n個(gè)正整數(shù)，分別代表xk和yk在Tk內(nèi)的距離。

?

測試數(shù)據(jù)

(cat)dog;

dog cat

?

(dog,cat);

dog cat

測試輸出

1 2

?

生物學(xué)背景

? ? ? ? 進(jìn)化樹是在生物學(xué)中用來表示物種之間進(jìn)化關(guān)系的圖。進(jìn)化樹上每個(gè)葉子結(jié)點(diǎn)代表一個(gè)物種，葉子結(jié)點(diǎn)之間的距離可以代表兩個(gè)物種之間的差異程度。

? ? ? ??Newick是最常見的進(jìn)化樹文件格式。舉例來說，假如n個(gè)葉子結(jié)點(diǎn){v1, v2, … ,vn}都和一個(gè)內(nèi)結(jié)點(diǎn)u相連，如下圖：

? ? ? ??用Newick格式可表示為(v1,v2,…,v3)u。Newick格式中結(jié)點(diǎn)名稱并非必須，我們可以省略一些或全部結(jié)點(diǎn)名稱，將這棵樹表示為(v1,v2,…,v3)，(, ,…,)u或者(, ,…,)。

? ? ? ??無根樹來說在用Newick表示時(shí)可以把任意結(jié)點(diǎn)作為根節(jié)點(diǎn)，因此可以用多個(gè)Newick描述同一棵樹，如(C, D, (A, B))和(A, (D, C), B)都表示下圖：

? ? ? ??對Newick格式更詳細(xì)地解釋可以參考這篇文章：https://blog.csdn.net/weixin_43569478/article/details/108079223

?

數(shù)學(xué)背景

? ? ? ??對于一棵樹中的兩個(gè)結(jié)點(diǎn)，肯定能找到唯一路徑將它們相連。因?yàn)槿绻麅牲c(diǎn)間有多條路徑，必然會形成循環(huán)，這與樹的定義不相符。

?

?

思路

? ? ? ??我在做本題時(shí)遇到了極大的困難，最終還是找了別人的代碼作弊通過。看了解釋后發(fā)現(xiàn)，我的思路大體方向是對的，只是很遺憾有關(guān)鍵點(diǎn)沒有想明白，導(dǎo)致長時(shí)間地被卡住，遲遲沒法更新這個(gè)系列。作為補(bǔ)償，這里給大家分享通過本題的兩種方法，我認(rèn)為都是值得掌握的。

?

方法1：善于利用現(xiàn)成的工具包

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

? ? ? ??在這里博主直接用了Biopython包Phylo模塊中的現(xiàn)成方法，這樣就無需自己動手對Newick格式進(jìn)行解析，Phylo模塊的介紹如下：https://biopython.org/wiki/Phylo

? ? ? ??盡管已經(jīng)找到現(xiàn)成的代碼，我還是在運(yùn)行時(shí)遇到了問題，在此記錄一下以防大家重復(fù)踩雷。情況是這樣的，我運(yùn)行Python時(shí)使用的IDE是PyCharm，安裝新模塊時(shí)就直接在File-Settings中下載安裝。然而這次安裝成功后，我用import Bio卻顯示錯(cuò)誤，似乎這個(gè)包沒有安裝上。經(jīng)過搜索后，我終于在下面這個(gè)網(wǎng)址：https://www.cnpython.com/qa/72896找到了解決方案，3樓的網(wǎng)友與我遇到的問題相同，按照他的方法把安裝文件夾路徑修改，并從“bio”重命名為“Bio”后，代碼正常運(yùn)行并獲得了正確結(jié)果。

?

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() # 把兩個(gè)葉子結(jié)點(diǎn)分開保存到x和y中
??? tree = Phylo.read(io.StringIO(i),'newick') # 用Biopython中Phylo模塊方法直接把Newick解析成樹結(jié)構(gòu)
??? clades = tree.find_clades()
??? for clade in clades: # 逐條給各分支加上長度，這里統(tǒng)一賦值為1
??????? clade.branch_length = 1
??? sys.stdout.write('%s' % tree.distance(x,y) + ' ') # 直接輸出兩結(jié)點(diǎn)間的距離，sys.stdout.write()輸出不會自動換行
sys.stdout.write('\n')

方法2：理解括號和逗號的含義

? ? ? ??借助方法1的代碼成功混入討論區(qū)后，我終于得以閱讀大佬們的解題思路，下面的代碼是我在學(xué)習(xí)后給出的，但有一些地方需要我們靜下心來理解。

? ? ? ??首先，要完成這道題并不需要寫一個(gè)方法把Newick格式解析為樹形結(jié)構(gòu)，再數(shù)出兩個(gè)結(jié)點(diǎn)的距離。我們需要把’(‘ , ’)’和’,’這些括號與樹結(jié)構(gòu)建立起聯(lián)系，這樣只需要關(guān)注這些符號的特點(diǎn)就可以計(jì)算出結(jié)點(diǎn)距離。

? ? ? ??從網(wǎng)站給出的例子，我們能對括號和逗號有一個(gè)直覺上的認(rèn)識，即：逗號代表并列，括號代表層級。比如：(A,B)，A和B的關(guān)系應(yīng)該如下圖左；(A)B應(yīng)該如下圖右。這似乎意味著兩個(gè)結(jié)點(diǎn)間多半個(gè)括號，距離應(yīng)加1，多一個(gè)逗號，距離應(yīng)加上2.

? ? ? ??不過，實(shí)際情況并非這么簡單。我們看下面這個(gè)示例：(A,( (()),()((,(,))) ),B)，A、B之間似乎遠(yuǎn)隔千山萬水，但實(shí)際它們的距離只有2，原因如下圖：

? ? ? ??可以看到，這里A、B中間的內(nèi)容都是配對的括號，形成一個(gè)完整的分支，與A、B同為一個(gè)父結(jié)點(diǎn)的不同分支，我們在計(jì)算A、B距離時(shí)盡可以把這個(gè)分支整個(gè)忽略掉。在計(jì)算A、B間的距離時(shí)，我們只需要關(guān)注A、B間不配對的括號即可，成功配對的括號全部消去，每多一個(gè)不匹配的括號，距離就加上1。

? ? ? ??對于逗號來說，我把它看成’)(‘，直接用’)(‘替換掉原文中的’,’，以上文的方法消去配對的括號，只剩下錯(cuò)配的括號，就是最終的結(jié)果。如何理解呢？我們畫幾幅圖來說明。

? ? ? ??假設(shè)一個(gè)進(jìn)化樹中沒有逗號，在刪去其中的配對括號后，變成了A))))B或A((((B這種形式，可以用圖表示為：

? ? ? ??顯然，A、B間的距離為4。

? ? ? ??現(xiàn)在我們讓情況復(fù)雜些，最后一個(gè)錯(cuò)配括號后面出現(xiàn)了一個(gè)逗號，如A)))),B，現(xiàn)在圖可以這樣表示：

? ? ? ??可以看到，加了一個(gè)逗號之后，距離增加了2。那為什么不直接計(jì)算多出逗號的數(shù)量，然后乘2加到最終結(jié)果里呢？我們再看一個(gè)例子，這次最后一個(gè)錯(cuò)配括號前也有逗號，如A),))),B，用圖可以表示為：

? ? ? ??可以看到，這種情況A、B間的距離依然是6。事實(shí)上，’,’的含義相當(dāng)于一個(gè)右括號和一個(gè)左括號。

?

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): # 將兩結(jié)點(diǎn)的位置賦給num1和num2
??????????? if i == target[0]:
??????????????? num1 = num
??????????? if i == target[1].replace('\n', ''):
??????????????? num2 = num
??????????? if num1 != -1 and num2 != -1: # 兩結(jié)點(diǎn)都找到就可以停止了
??????????????? break

??????? start = min(num1, num2) # 左節(jié)點(diǎn)賦值給start
??????? end = max(num1, num2) # 右節(jié)點(diǎn)賦值給end

??????? left = right = 0
??????? for i in range(start + 1, end): # 從左節(jié)點(diǎn)像右節(jié)點(diǎn)遍歷
??????????? # print(newdata[i])
??????????? if newdata[i] in ',)': # 碰到右括號含義的字符
??????????????? if left > 0: # 如果左邊有與其配對的左括號
??????????????????? left -= 1 # 刪去左括號
??????????????? else: # 如果左邊沒有與其配對的左括號
??????????????????? right += 1 # 記錄這個(gè)錯(cuò)配右括號
??????????? 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保存各對結(jié)點(diǎn)
for num, i in enumerate(input): # 根據(jù)所給輸入數(shù)據(jù)的格式
??? if num % 3 == 0: # 第3k行是進(jìn)化樹
??????? trees.append(i)
??? if num % 3 == 1: # 第3k+1行是結(jié)點(diǎn)，第3k+2行是空行
??????? pairs.append(i.split(' '))

result = [] # 用result存儲距離
for i in range(len(trees)):
??? result.append(nwck(trees[i], pairs[i])) # 調(diào)用newick函數(shù)解決問題

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