【ROSALIND】【练Python，学生信】 47 默慈金数和RNA二级结构
如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面 谢谢配合~
题目：
默慈金数和RNA二级结构
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所有不交叉匹配情况的数目，结果对1,000,000取模。
测试数据
>Rosalind_57
AUAU
测试输出
7
生物学背景
在33 卡特兰数与RNA二级结构中，我们学会了如何计算不交叉完美匹配的数量。但在真实世界中，并非所有的碱基都会参与互补配对，就像在一次派对中并非每个人都能找到人握手。所以本题在之前的基础上又去掉了所有碱基都参与配对的条件，要求计算的是所有情况的总数。
数学背景
默慈金数（Motzkin number）的介绍可以参考下面这篇博文（https://www.cnblogs.com/yaoyueduzhen/p/5456530.html）。把默慈金数记为M，第n个数为Mn，举例来说，下图是M5=21的形象展示。可以看到，默慈金数把5个结点之间无配对、一个配对、两个配对的情况都列举了出来，总数就是第5个默慈金数的大小。
如何计算默慈金数呢？与卡特兰数相似，默慈金数也有递推公式，我们假设M0=M1=0，假设有n个结点Kn。先考虑第一个结点K1，它可能参与配对，也可能不参与，如果参与配对，那它可以与后面n-1个结点中的任意一个配对；假设第1个结点和第k（2≤k≤n）个结点配对，那么所有结点就被分成了两部分，k-2个在1和k之间，n-k个在1和k外面，这两部分之间不可相连，否则就出现交叉了；这样剩下的结点就有Mk−2⋅Mn−k种方式连接。以此类推，我们得到了递推公式：
上面是用任意两者都可配对的结点做示范，具体到RNA序列本质也是一样的，只是要另外考虑碱基互补配对。如下图是序列UAGCGUGAUCAC的不交叉配对的其中两种情况：
思路
本题默慈金数计算代码来自
https://github.com/fernandoBRS/Rosalind-Problems/blob/master/motz.py
这道题的思路和卡特兰数那道题相似，都是借用这种递推公式求解可能二级结构的数量。因此可借用33 卡特兰数与RNA二级结构的结构，修改部分判断条件，把递推公式替换为默慈金数即可。
代码
memorize = {} # 建立一个字典，存储已经出现过的字符串及不交叉匹配的数量
memorize[''] = 1 # 如果序列为空，说明只有这一种情况
def ismatch(c1, c2):
"""判断是否碱基配对的函数"""
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: # 如果序列只有两个字符，直接判断是否相等即可
if ismatch(seq[0], seq[1]):
return 2
# 当序列长度大于2时，代入默慈金数公式
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 # 将新出现的字符串存入字典方便直接查询
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)
