【ROSALIND】【练Python，学生信】14 寻找DNA模序
如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面 谢谢配合~
题目：
寻找DNA共有模序（Finding a Shared Motif）
Given: A collection of k (k<100) DNA strings of length at most 1 kbp, each in FASTA format.
所给：不超过100条DNA序列，长度均不超过1kb，以FASTA格式给出。
Return: A longest common substring of the collection. (If multiple solutions exist, you may return any single solution.)
需得：这些DNA序列的最长公共子串（如果有多个，返回任意一个即可）。
测试数据
>Rosalind_1
GATTACA
>Rosalind_2
TAGACCA
>Rosalind_3
ATACA
测试输出
AC
背景
在09 确定DNA子序列出现的位置中，我们练习了从DNA序列中寻找给定子序列位置的方法，但是在实际应用中，我们并不能提前知道哪些是有意义的模序，需要通过比较多条序列以找到公共序列。
本题假设模序不发生突变，我们需要在一组DNA序列中找到每个序列共有的最长子序列。
思路
对本题我的想法是任意取一条DNA序列，将其打成稀碎的碎片，这里“稀碎”指每一个打碎方式都要出现。比如对序列‘AGCT’，打碎后出现的碎片有：‘AGCT’、‘AGC’、‘GCT’、‘AG’、‘GC’、‘CT’，单个碱基的因为没有意义所以不需要列出。如果一个碎片是所有DNA序列的公共子序列，那么它一定在每条序列中都出现，所以我要做的是在剩下的序列中一一检查各碎片是否存在，都在则为公共子序列。
我的代码也是按这个思路实现的。我写了一个打碎序列的函数，在其中用了两个循环以罗列出所有长度大于一的碎片。随后我将第一条序列输入该函数，得到一个储存了所有碎片的列表，为了方便取最长子串和提高比对效率，我用.sort()方法对碎片进行了由长到段的排列。之后就是挨个取碎片与序列比对。用一个标志flag记录碎片在几条序列中出现，若比对完flag与序列条数相等，说明所有序列都有该碎片，即找到了一****序列。
这段代码还有多处可以优化的地方，比如选择最短序列进行打碎、找到一条序列即可停止等。此外网上还有很多更简洁的编码方式和解题思路，在此不作引述，仅展示我自己的原始想法。
代码
def readfasta(lines):
'''读取FASTA文件的函数'''
seq = []
index = []
seqplast = ""
numlines = 0
for i in lines:
if '>' in i:
index.append(i.replace("\n", "").replace(">", ""))
seq.append(seqplast.replace("\n", ""))
seqplast = ""
numlines += 1
else:
seqplast = seqplast + i.replace("\n", "")
numlines += 1
if numlines == len(lines):
seq.append(seqplast.replace("\n", ""))
seq = seq[1:]
return index, seq
def frag(seq):
'''将一条字符串打成长度不等的碎片的函数'''
res = []
i = 0
while i < len(seq):
s_seq = seq[i:]
j = 1
while j < len(s_seq):
res.append(s_seq[:(len(s_seq)-j+1)])
j += 1
i += 1
return res
f = open('input.txt', 'r')
lines = f.readlines()
f.close()
[index, seq] = readfasta(lines)
frags = frag(seq[0])
frags.sort(key=len, reverse=True) # 将碎片按从长到短排列
i = 0
while i < len(frags):
flag = 0
for j in seq:
r = j.count(frags[i]) # .count()方法返回碎片在字符串中出现的次数
if r != 0: # 若碎片出现过，flag值加一
flag += 1
if flag >= len(seq):
print(frags[i]) # 若flag与序列数相同，说明为公共子串，符合要求，输出
break
i += 1
