【ROSALIND】【練Python,學生信】33 卡特蘭數(shù)與RNA二級結構

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

題目:
卡特蘭數(shù)和RNA二級結構(Catalan Numbers and RNA Secondary Structures)
Given: An RNA string s having the same number of occurrences of 'A' as 'U' and the same number of occurrences of 'C' as 'G'. The length of the string is at most 300 bp.
所給:一個不超過300bp的RNA序列s,包含數(shù)目相同的A和U,以及數(shù)目相同的G和C。
Return: The total number of noncrossing perfect matchings of basepair edges in the bonding graph of s, modulo 1,000,000.
需得:由s得到的所有不相交完美匹配的數(shù)目,結果對1,000,000取模。
?
測試數(shù)據(jù)
>Rosalind_57
AUAU
測試輸出
2
?
生物學背景
????????在26 RNA二級結構與完美匹配中,我們初步了解了圖論與RNA折疊的關系,并計算了完美匹配。但在實際中,不是所有的完美匹配都可能存在。在用數(shù)學方法進行RNA二級結構預測時,通常禁止配對堿基間的相互交叉。從而避免形成一種名為“假結(pseudoknot)”的折疊結構,如下圖:

數(shù)學背景
????????我們假設有一條長為n的RNA序列,根據(jù)這條序列可繪制出一個有n個結點的bonding graph。為了避免RNA折疊中出現(xiàn)堿基相互交叉的情況,需要兩條邊{i, j},{k,l}滿足i<k<j<l這種情況不存在。下圖就是26 RNA二級結構與完美匹配例子中唯一符合不交叉要求的兩個完美匹配:

????????為了解決這個問題,我們首先要引入卡特蘭數(shù)(catalan)的概念,這是一個經常在各種計數(shù)問題中現(xiàn)身的數(shù)組。遞歸公式為:h(n)= h(0)*h(n-1)+h(1)*h(n-2) + ... + h(n-1)h(0) (其中n>=2,h(0) = h(1) = 1)。我們可以用Python寫一個計算卡特蘭數(shù)各項的代碼,如下所示:
def catalan(n):
??? """返回卡特蘭數(shù)的函數(shù)"""
??? if n <= 1:
??????? return 1
??? else:
??????? h = []
??????? h.append(1)
??????? h.append(1)
??????? i = 2
??????? while i <= n:
??????????? j = 0
??????????? tmp = 0
??????????? while j < i:
??????????????? tmp += h[j] * h[i - 1 - j]
??????????????? j += 1
??????????? h.append(tmp)
????????? ??i += 1
??????? return h[n]
?
????????卡特蘭數(shù)與本題有什么關系呢?我們假設有一個RNA序列s,長為2n,要求其所有不交叉完美匹配的數(shù)量。既然是完美匹配,第一個字符,即s[0]需要與剩下2n-1個字符中的一個匹配,不妨設其與s[k]配對,同時,為了達成完美匹配,k與0之間必須有偶數(shù)個字符。這樣我們就把s分成了兩部分,第一部分是s[1:k],第二部分是s[k+1:],這兩部分可以看成兩個相對獨立的、具有偶數(shù)個字符的序列。由于兩個序列獨立,要求總的不交叉完美匹配數(shù)量就是把兩個序列分別的不交叉完美匹配數(shù)量相乘。這兩個序列又可以用同樣的思路進行分析,以此類推,直到所有堿基都配對為止。這個過程可以用卡特蘭數(shù)的遞推公式求解。
?
思路
????????本題程序編寫思路與求卡特蘭數(shù)各項的函數(shù)相似。下面代碼用遞歸的方法實現(xiàn),不再詳述。需要注意的是,完整計算一個序列的所有匹配會消耗大量的時間,但實際上一個長序列中會有很多重復的短序列,每次都要計算匹配數(shù)未免過于浪費。所以建立一個自動儲存每個序列的匹配數(shù)量,下次遇到時直接取出即可,可以大大縮減程序運行時間。
?
代碼
memorize = {} # 建立一個字典,存儲已經出現(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(): # 如果這段序列之前已經被分析過了,直接取出結果即可
??????? return memorize[seq]
??? if len(seq) % 2 == 1: # 如果這個序列長度是奇數(shù),不可能存在完美匹配
??????? memorize[seq] = 0
??????? return 0
??? if seq.count('A') != seq.count('U') or seq.count('G') != seq.count('G'): # 如果這個序列中配對的堿基數(shù)量不相同,不可能存在完美匹配
??????? memorize[seq] = 0
??????? return 0
??? i = 1
??? num = 0
??? while i < len(seq): # 在序列中搜索所有可以與第一個堿基配對的堿基
??????? if ismatch(seq[0], seq[i]) == 1: # 如果第i個堿基配對
??????????? num += (noncross(seq[1:i]) * noncross(seq[i+1:])) # 去檢驗被第k個堿基分出的兩個序列
??????? i += 2 # 只需從第偶數(shù)個堿基中搜索
??? memorize[seq] = num # 記錄這個序列的不交叉完美匹配結果
??? 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)
?