【ROSALIND】【練Python，學(xué)生信】 47 默慈金數(shù)和RNA二級結(jié)構(gòu)
如果第一次閱讀本系列文檔請先移步閱讀【ROSALIND】【練Python，學(xué)生信】00 寫在前面 ?謝謝配合~
題目：
默慈金數(shù)和RNA二級結(jié)構(gòu)
Given: An RNA string s of length at most 300 bp.
所給：一條RNA序列s，長度不超過300bp。
Return: The total number of noncrossing matchings of basepair edges in the bonding graph of s, modulo 1,000,000.
需得：s所有不交叉匹配情況的數(shù)目，結(jié)果對1,000,000取模。
?
測試數(shù)據(jù)
>Rosalind_57
AUAU
測試輸出
7
?
生物學(xué)背景
? ? ? ? 在33 卡特蘭數(shù)與RNA二級結(jié)構(gòu)中，我們學(xué)會了如何計算不交叉完美匹配的數(shù)量。但在真實世界中，并非所有的堿基都會參與互補(bǔ)配對，就像在一次派對中并非每個人都能找到人握手。所以本題在之前的基礎(chǔ)上又去掉了所有堿基都參與配對的條件，要求計算的是所有情況的總數(shù)。
?
數(shù)學(xué)背景
? ? ? ??默慈金數(shù)（Motzkin number）的介紹可以參考下面這篇博文（https://www.cnblogs.com/yaoyueduzhen/p/5456530.html）。把默慈金數(shù)記為M，第n個數(shù)為Mn，舉例來說，下圖是M5=21的形象展示?？梢钥吹?，默慈金數(shù)把5個結(jié)點之間無配對、一個配對、兩個配對的情況都列舉了出來，總數(shù)就是第5個默慈金數(shù)的大小。
? ? ? ? 如何計算默慈金數(shù)呢？與卡特蘭數(shù)相似，默慈金數(shù)也有遞推公式，我們假設(shè)M0=M1=0，假設(shè)有n個結(jié)點Kn。先考慮第一個結(jié)點K1，它可能參與配對，也可能不參與，如果參與配對，那它可以與后面n-1個結(jié)點中的任意一個配對；假設(shè)第1個結(jié)點和第k（2≤k≤n）個結(jié)點配對，那么所有結(jié)點就被分成了兩部分，k-2個在1和k之間，n-k個在1和k外面，這兩部分之間不可相連，否則就出現(xiàn)交叉了；這樣剩下的結(jié)點就有Mk?2?Mn?k種方式連接。以此類推，我們得到了遞推公式：
? ? ? ??上面是用任意兩者都可配對的結(jié)點做示范，具體到RNA序列本質(zhì)也是一樣的，只是要另外考慮堿基互補(bǔ)配對。如下圖是序列UAGCGUGAUCAC的不交叉配對的其中兩種情況：
?
思路
? ? ? ? 本題默慈金數(shù)計算代碼來自
https://github.com/fernandoBRS/Rosalind-Problems/blob/master/motz.py
? ? ? ? 這道題的思路和卡特蘭數(shù)那道題相似，都是借用這種遞推公式求解可能二級結(jié)構(gòu)的數(shù)量。因此可借用33 卡特蘭數(shù)與RNA二級結(jié)構(gòu)的結(jié)構(gòu)，修改部分判斷條件，把遞推公式替換為默慈金數(shù)即可。
?
?
代碼
memorize = {} # 建立一個字典，存儲已經(jīng)出現(xiàn)過的字符串及不交叉匹配的數(shù)量
memorize[''] = 1 # 如果序列為空，說明只有這一種情況
def ismatch(c1, c2):
??? """判斷是否堿基配對的函數(shù)"""
??? if (c1 == 'A' and c2 == 'U') or (c1 == 'U' and c2 == 'A') or \
??????????? (c1 == 'G' and c2 == 'C') or (c1 == 'C' and c2 == 'G'):
??????? return 1
??? else:
??????? return 0
def noncross(seq):
??? """判斷是否有不交叉的匹配"""
??? if seq in memorize.keys(): # 如果這段序列之前已經(jīng)被分析過了，直接取出結(jié)果即可
??????? return memorize[seq]
??? if len(seq) == 2: # 如果序列只有兩個字符，直接判斷是否相等即可
??????? if ismatch(seq[0], seq[1]):
??????????? return 2
??? # 當(dāng)序列長度大于2時，代入默慈金數(shù)公式
??? num = noncross(seq[1:]) + sum([ismatch(seq[0], seq[k]) * noncross(seq[1: k]) *
?????????????????????????????????? noncross(seq[k + 1:]) for k in range(1, len(seq))]) ???
? ? memorize[seq] = num # 將新出現(xiàn)的字符串存入字典方便直接查詢
??? return num
f = open('input.txt', 'r')
input = f.readlines()
f.close()
index = input[0].replace('\n', '')
input = input[1:]
i = 0
seq = ''
while i < len(input):
??? seq = seq + input[i].replace('\n', '')
??? i += 1
res = noncross(seq)
print(res % 1000000)
?
