【ROSALIND】【练Python，学生信】37 KMP算法与模序搜索

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

题目：

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

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

所给：长度不超过100kb的一条DNA序列s，以FASTA给出。

Return: The failure array of s.

需得：s的“失败数列”。

测试数据

>Rosalind_87

CAGCATGGTATCACAGCAGAG

测试输出

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

背景

        在之前的题目中，我们曾经接触过从长的DNA序列中搜索短的模序的问题。由于真核生物的基因组太过巨大，用之前的方法势必消耗大量的时间。我们如何改进呢？

        在之前的搜索策略中，我们的做法是用一个窗口在长的序列上滑动。表现为当从长序列的某一位与短序列第一位相同时，从这一位开始依次比较长短序列，比较完成后还要退回到长序列开始那一位的后一位重新开始这个过程，但这中有一些比较是可以避免的。

        举个例子，假设我们要在长序列s=TAGGTACGTACGGCATCACG上搜索t=ACGTACGT，s[6]到s[12]之间的字符与t[1]到t[7]之间的字符相匹配，然而s[13]与t[8]不相同。按传统方法我们要退回s[7]重新进行比较，但显然s[7]=C，s[8]=G，s[9]=T都不同于t[1]=A；而对于s[10]，t[1:4]=t[5:8]=ACGT，而s[13]=G，t[8]=T，使错配再次出现。基于以上事实，我们可以直接去扫描s[14]，不需要向回滑动。上述的思想来自于KMP算法（Knuth-Morris-Pratt algorithm）。

        在KMP算法中有一个重要概念叫failure array。对一个长为n的序列s，前缀是子串s[1:j]，后缀是子串s[k:n]。而failure array是一个长为n的数组P，P[k]是与前缀s[1:k−j+1] 相同的最长子串s[j:k]的长度，其中j不为1，否则P[k]就总与k相等了。另外我们假定P[1]=0。

思路

        我这里还是用两重循环解决这个问题。第一重循环从头开始扫描整个序列，与前缀字符比较。一旦找到与第一个字符相同的字符，从这里进入第二重循环，把其后与前缀相同的字符找出来。

代码

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

[index, seq] = readfasta(lines)
num = [0] * len(seq[0]) # num即存储failure array的列表，初始都设为0

i = 1 # i扫描序列
while i < len(seq[0]):
    tmp = 0 # tmp记录重合字符的范围
    j = 1 # 首先要在序列中找到和第一个字符相同的字符位置
    if seq[0][i:i + j] == seq[0][0:j]: # 找到与第一个字符相同的字符位置
        tmp += 1
        if tmp > num[i]: # 如果tmp比此时列表中i位置的数字大，则替代原来的数字
            num[i] = tmp
        if i == 2: # 处理开始两个字符重复的特殊情况
            if seq[0][0:2] == seq[0][1:3]:
                tmp += 1
                num[i] = tmp
        while i + j < len(seq[0]): # 窗口的最后一个字符滑动到序列末尾前
            j += 1
            if i == 1: # 处理开始两个字符重复的特殊情况
                break
            else:
                if seq[0][0:j] == seq[0][i:i + j]: # 如果后续字符相同，tmp继续累加
                    tmp += 1
                    if tmp > num[i+j-1]:
                        num[i+j-1] = tmp
                else: # 如果后续字母不再相同，终止当前循环
                    break
    i += 1

# print(num)

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

f.close()