如果第一次閱讀本系列文檔請(qǐng)先移步閱讀【ROSALIND】【練Python，學(xué)生信】00 寫在前面 ?謝謝配合~

題目：

兩個(gè)交錯(cuò)的模序（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的最短公共超序列。如果有多個(gè)解，給出一個(gè)即可。

?

測(cè)試數(shù)據(jù)

ATCTGAT

TGCATA

測(cè)試輸出

ATGCATGAT

?

背景

? ? 在38 求解最長(zhǎng)公共子序列中我們計(jì)算了兩條序列間的最長(zhǎng)公共子序列，從而幫助我們找到分散在不同外顯子中的模序。本題我們來(lái)求解這個(gè)反演問(wèn)題：兩條序列分別包含兩個(gè)模序，如何把它們交織在一起，產(chǎn)生一個(gè)同時(shí)包含這兩個(gè)序列的最短超序列。

?

思路

? ?這里的思路和38 求解最長(zhǎng)公共子序列一樣，代碼也大部分是原樣復(fù)制。只是本題的目的發(fā)生了改變，需求的是最短公共超序列，而這個(gè)最短公共超序列其實(shí)就是最長(zhǎng)公共子序列加上兩個(gè)序列中不在公共子序列中的其他字符。所以我在這里記下兩個(gè)序列公共子序列每個(gè)字符的位置，把兩個(gè)序列其他位置的字符順次插入公共字符之間，即可得到最短公共超序列其中一種情況，具體方法見代碼。

? ? 用動(dòng)態(tài)規(guī)劃求解最長(zhǎng)公共子序列的思路請(qǐng)?jiān)斠?a target="_blank" >24 求解最長(zhǎng)遞增子序列。

?

代碼

def LCSQ(s1, s2):
??? """用動(dòng)態(tài)規(guī)劃求解最長(zhǎng)公共子序列的函數(shù)"""
??? 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分別記錄公共字符在兩個(gè)序列中的位置
??????? 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是分別是這個(gè)公共字符在兩個(gè)序列中的位置
??????????????? 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是最長(zhǎng)公共子序列的長(zhǎng)度
result = [] # result記錄最短公共超序列
for i in range(num):
??? l1 = loc1[i]
??? l2 = loc2[i]
??? if i == 0: # 當(dāng)i是第一個(gè)公共字符時(shí)
??????? # 把第一個(gè)公共字符前的兩序列內(nèi)容以及第一個(gè)公共字符本身寫入result
??????? result = seq1[:l1-1] + seq2[:l2-1] + res[i]
??? else: # 當(dāng)i不是第一個(gè)公共字符時(shí)
??????? # 把第i-1個(gè)和第i個(gè)公共字符之間的內(nèi)容以及第i個(gè)公共字符本身寫入result
??????? result = result + seq1[loc1[i-1]:l1 - 1] + seq2[loc2[i-1]:l2 - 1] + res[i]

# 把公共字符結(jié)束一直到兩序列結(jié)束的內(nèi)容寫入result
result = result + seq1[loc1[i]:] + seq2[loc2[i]:]
print(result)

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

?