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

題目：

更快速地模序搜索策略（Speeding Up Motif Finding）

Given: A DNA string s (of length at most 100 kbp) in FASTA format.

所給：長(zhǎng)度不超過(guò)100kb的一條DNA序列s，以FASTA給出。

Return: The failure array of s.

需得：s的“失敗數(shù)列”。

?

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

>Rosalind_87

CAGCATGGTATCACAGCAGAG

測(cè)試輸出

0 0 0 1 2 0 0 0 0 0 0 1 2 1 2 3 4 5 3 0 0

?

背景

????????在之前的題目中，我們?cè)?jīng)接觸過(guò)從長(zhǎng)的DNA序列中搜索短的模序的問(wèn)題。由于真核生物的基因組太過(guò)巨大，用之前的方法勢(shì)必消耗大量的時(shí)間。我們?nèi)绾胃倪M(jìn)呢？

????????在之前的搜索策略中，我們的做法是用一個(gè)窗口在長(zhǎng)的序列上滑動(dòng)。表現(xiàn)為當(dāng)從長(zhǎng)序列的某一位與短序列第一位相同時(shí)，從這一位開(kāi)始依次比較長(zhǎng)短序列，比較完成后還要退回到長(zhǎng)序列開(kāi)始那一位的后一位重新開(kāi)始這個(gè)過(guò)程，但這中有一些比較是可以避免的。

????????舉個(gè)例子，假設(shè)我們要在長(zhǎng)序列s=TAGGTACGTACGGCATCACG上搜索t=ACGTACGT，s[6]到s[12]之間的字符與t[1]到t[7]之間的字符相匹配，然而s[13]與t[8]不相同。按傳統(tǒng)方法我們要退回s[7]重新進(jìn)行比較，但顯然s[7]=C，s[8]=G，s[9]=T都不同于t[1]=A；而對(duì)于s[10]，t[1:4]=t[5:8]=ACGT，而s[13]=G，t[8]=T，使錯(cuò)配再次出現(xiàn)?；谝陨鲜聦?shí)，我們可以直接去掃描s[14]，不需要向回滑動(dòng)。上述的思想來(lái)自于KMP算法（Knuth-Morris-Pratt algorithm）。

????????在KMP算法中有一個(gè)重要概念叫failure array。對(duì)一個(gè)長(zhǎng)為n的序列s，前綴是子串s[1:j]，后綴是子串s[k:n]。而failure array是一個(gè)長(zhǎng)為n的數(shù)組P，P[k]是與前綴s[1:k?j+1] 相同的最長(zhǎng)子串s[j:k]的長(zhǎng)度，其中j不為1，否則P[k]就總與k相等了。另外我們假定P[1]=0。

?

思路

????????我這里還是用兩重循環(huán)解決這個(gè)問(wèn)題。第一重循環(huán)從頭開(kāi)始掃描整個(gè)序列，與前綴字符比較。一旦找到與第一個(gè)字符相同的字符，從這里進(jìn)入第二重循環(huán)，把其后與前綴相同的字符找出來(lái)。

?

代碼

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

[index, seq] = readfasta(lines)
num = [0] * len(seq[0]) # num即存儲(chǔ)failure array的列表，初始都設(shè)為0

i = 1 # i掃描序列
while i < len(seq[0]):
??? tmp = 0 # tmp記錄重合字符的范圍
??? j = 1 # 首先要在序列中找到和第一個(gè)字符相同的字符位置
??? if seq[0][i:i + j] == seq[0][0:j]: # 找到與第一個(gè)字符相同的字符位置
??????? tmp += 1
??????? if tmp > num[i]: # 如果tmp比此時(shí)列表中i位置的數(shù)字大，則替代原來(lái)的數(shù)字
??????????? num[i] = tmp
??????? if i == 2: # 處理開(kāi)始兩個(gè)字符重復(fù)的特殊情況
??????????? if seq[0][0:2] == seq[0][1:3]:
??????????????? tmp += 1
??????????????? num[i] = tmp
??????? while i + j < len(seq[0]): # 窗口的最后一個(gè)字符滑動(dòng)到序列末尾前
??????????? j += 1
??????????? if i == 1: # 處理開(kāi)始兩個(gè)字符重復(fù)的特殊情況
??????????????? break
??????????? else:
??????????????? if seq[0][0:j] == seq[0][i:i + j]: # 如果后續(xù)字符相同，tmp繼續(xù)累加
??????????????????? tmp += 1
??????????????????? if tmp > num[i+j-1]:
??????????????????????? num[i+j-1] = tmp
??????????????? else: # 如果后續(xù)字母不再相同，終止當(dāng)前循環(huán)
??????????????????? break
??? i += 1

# print(num)

f = open('output.txt', 'a')
for i in num:
??? f.write(str(i) + ' ')

f.close()