【ROSALIND】【练Python，学生信】33 卡特兰数与RNA二级结构

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

题目：

卡特兰数和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，包含数目相同的A和U，以及数目相同的G和C。

Return: The total number of noncrossing perfect matchings of basepair edges in the bonding graph of s, modulo 1,000,000.

需得：由s得到的所有不相交完美匹配的数目，结果对1,000,000取模。

测试数据

>Rosalind_57

AUAU

测试输出

2

生物学背景

        在26 RNA二级结构与完美匹配中，我们初步了解了图论与RNA折叠的关系，并计算了完美匹配。但在实际中，不是所有的完美匹配都可能存在。在用数学方法进行RNA二级结构预测时，通常禁止配对碱基间的相互交叉。从而避免形成一种名为“假结（pseudoknot）”的折叠结构，如下图：

数学背景

        我们假设有一条长为n的RNA序列，根据这条序列可绘制出一个有n个结点的bonding graph。为了避免RNA折叠中出现碱基相互交叉的情况，需要两条边{i, j}，{k,l}满足i<k<j<l这种情况不存在。下图就是26 RNA二级结构与完美匹配例子中唯一符合不交叉要求的两个完美匹配：

        为了解决这个问题，我们首先要引入卡特兰数（catalan）的概念，这是一个经常在各种计数问题中现身的数组。递归公式为：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写一个计算卡特兰数各项的代码，如下所示：

def catalan(n):
    """返回卡特兰数的函数"""

    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]

        卡特兰数与本题有什么关系呢？我们假设有一个RNA序列s，长为2n，要求其所有不交叉完美匹配的数量。既然是完美匹配，第一个字符，即s[0]需要与剩下2n-1个字符中的一个匹配，不妨设其与s[k]配对，同时，为了达成完美匹配，k与0之间必须有偶数个字符。这样我们就把s分成了两部分，第一部分是s[1:k]，第二部分是s[k+1:]，这两部分可以看成两个相对独立的、具有偶数个字符的序列。由于两个序列独立，要求总的不交叉完美匹配数量就是把两个序列分别的不交叉完美匹配数量相乘。这两个序列又可以用同样的思路进行分析，以此类推，直到所有碱基都配对为止。这个过程可以用卡特兰数的递推公式求解。

思路

        本题程序编写思路与求卡特兰数各项的函数相似。下面代码用递归的方法实现，不再详述。需要注意的是，完整计算一个序列的所有匹配会消耗大量的时间，但实际上一个长序列中会有很多重复的短序列，每次都要计算匹配数未免过于浪费。所以建立一个自动储存每个序列的匹配数量，下次遇到时直接取出即可，可以大大缩减程序运行时间。

代码

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 == 1: # 如果这个序列长度是奇数，不可能存在完美匹配
        memorize[seq] = 0
        return 0

    if seq.count('A') != seq.count('U') or seq.count('G') != seq.count('G'): # 如果这个序列中配对的碱基数量不相同，不可能存在完美匹配
        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 # 只需从第偶数个碱基中搜索
    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)