【ROSALIND】【练Python，学生信】49 求解最短公共超序列

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

题目：

两个交错的模序（Interleaving Two Motifs）

Given: Two DNA strings s and t.

所给：两条DNA序列s和t。

Return: A shortest common supersequence of s and t. If multiple solutions exist, you may output any one.

需得：一条s和t的最短公共超序列。如果有多个解，给出一个即可。

测试数据

ATCTGAT

TGCATA

测试输出

ATGCATGAT

背景

    在38 求解最长公共子序列中我们计算了两条序列间的最长公共子序列，从而帮助我们找到分散在不同外显子中的模序。本题我们来求解这个反演问题：两条序列分别包含两个模序，如何把它们交织在一起，产生一个同时包含这两个序列的最短超序列。

思路

   这里的思路和38 求解最长公共子序列一样，代码也大部分是原样复制。只是本题的目的发生了改变，需求的是最短公共超序列，而这个最短公共超序列其实就是最长公共子序列加上两个序列中不在公共子序列中的其他字符。所以我在这里记下两个序列公共子序列每个字符的位置，把两个序列其他位置的字符顺次插入公共字符之间，即可得到最短公共超序列其中一种情况，具体方法见代码。

    用动态规划求解最长公共子序列的思路请详见24 求解最长递增子序列。

代码

def LCSQ(s1, s2):
    """用动态规划求解最长公共子序列的函数"""
    m, n = len(s1), len(s2)
    dp = [[["", 0] for j in list(range(n + 1))] for i in list(range(m + 1))]
    for i in list(range(1, m + 1)):
        dp[i][0][0] = s1[i - 1]
    for j in list(range(1, n + 1)):
        dp[0][j][0] = s2[j - 1]

        for i in list(range(1, m + 1)):
            for j in list(range(1, n + 1)):
                if s1[i - 1] == s2[j - 1]:
                    dp[i][j] = ['↖', dp[i - 1][j - 1][1] + 1]
                elif dp[i][j - 1][1] > dp[i - 1][j][1]:
                    dp[i][j] = ['←', dp[i][j - 1][1]]
                else:
                    dp[i][j] = ['↑', dp[i - 1][j][1]]

        s3 = []
        loc1 = []  # loc1和loc2分别记录公共字符在两个序列中的位置
        loc2 = []
        i = m
        j = n
        while i > 0 and j > 0:
            if dp[i][j][0] == '↖':
                s3.append(dp[i][0][0])
                loc1.append(i) # i、j是分别是这个公共字符在两个序列中的位置
                loc2.append(j)
                i -= 1
                j -= 1
            elif dp[i][j][0] == '←':
                j -= 1
            elif dp[i][j][0] == '↑':
                i -= 1
        s3 = s3[::-1]
        return s3, loc1[::-1], loc2[::-1]


f = open('rosalind_scsp.txt', 'r')
lines = f.readlines()
f.close()

seq1 = lines[0].strip()
seq2 = lines[1]
res = []
loc1 = loc2 = []
[res, loc1, loc2] = LCSQ(seq1, seq2)
res = ''.join(res)

num = len(loc1) # num是最长公共子序列的长度
result = [] # result记录最短公共超序列
for i in range(num):
    l1 = loc1[i]
    l2 = loc2[i]
    if i == 0: # 当i是第一个公共字符时
        # 把第一个公共字符前的两序列内容以及第一个公共字符本身写入result
        result = seq1[:l1-1] + seq2[:l2-1] + res[i]
    else: # 当i不是第一个公共字符时
        # 把第i-1个和第i个公共字符之间的内容以及第i个公共字符本身写入result
        result = result + seq1[loc1[i-1]:l1 - 1] + seq2[loc2[i-1]:l2 - 1] + res[i]

# 把公共字符结束一直到两序列结束的内容写入result
result = result + seq1[loc1[i]:] + seq2[loc2[i]:]
print(result)

f = open('output.txt', 'w')
f.write(result)
f.close()