【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序列的最長公共子串（如果有多個，返回任意一個即可）。
?
測試數(shù)據(jù)
>Rosalind_1
GATTACA
>Rosalind_2
TAGACCA
>Rosalind_3
ATACA
測試輸出
AC
?
背景
在09 確定DNA子序列出現(xiàn)的位置中，我們練習了從DNA序列中尋找給定子序列位置的方法，但是在實際應用中，我們并不能提前知道哪些是有意義的模序，需要通過比較多條序列以找到公共序列。
本題假設模序不發(fā)生突變，我們需要在一組DNA序列中找到每個序列共有的最長子序列。
?
思路
對本題我的想法是任意取一條DNA序列，將其打成稀碎的碎片，這里“稀碎”指每一個打碎方式都要出現(xiàn)。比如對序列‘AGCT’，打碎后出現(xiàn)的碎片有：‘AGCT’、‘AGC’、‘GCT’、‘AG’、‘GC’、‘CT’，單個堿基的因為沒有意義所以不需要列出。如果一個碎片是所有DNA序列的公共子序列，那么它一定在每條序列中都出現(xiàn)，所以我要做的是在剩下的序列中一一檢查各碎片是否存在，都在則為公共子序列。
我的代碼也是按這個思路實現(xiàn)的。我寫了一個打碎序列的函數(shù)，在其中用了兩個循環(huán)以羅列出所有長度大于一的碎片。隨后我將第一條序列輸入該函數(shù)，得到一個儲存了所有碎片的列表，為了方便取最長子串和提高比對效率，我用.sort()方法對碎片進行了由長到段的排列。之后就是挨個取碎片與序列比對。用一個標志flag記錄碎片在幾條序列中出現(xiàn)，若比對完flag與序列條數(shù)相等，說明所有序列都有該碎片，即找到了一****序列。
這段代碼還有多處可以優(yōu)化的地方，比如選擇最短序列進行打碎、找到一條序列即可停止等。此外網(wǎng)上還有很多更簡潔的編碼方式和解題思路，在此不作引述，僅展示我自己的原始想法。
?
代碼
def readfasta(lines):
??? '''讀取FASTA文件的函數(shù)'''
??? 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):
??? '''將一條字符串打成長度不等的碎片的函數(shù)'''
??? 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()方法返回碎片在字符串中出現(xiàn)的次數(shù)
??????? if r != 0:? # 若碎片出現(xiàn)過，flag值加一
??????????? flag += 1
??????? if flag >= len(seq):
??????????? print(frags[i])? # 若flag與序列數(shù)相同，說明為公共子串，符合要求，輸出
??????????? break
??? i += 1
