我是Python和生物信息学的新手,但我正在通过rosalind.info网站学习两种方法.您可以使用后缀树执行此操作.后缀树(见
http://en.wikipedia.org/wiki/Suffix_tree)是神奇的数据结构,它使生物信息学中的所有事情成为可能.您可以快速找到多个长序列中的公共子串.后缀树只需要线性时间,因此长度10,000,000是可行的.
所以首先找到序列的反向补码.然后将两者放入后缀树中,找到它们之间的公共子串(最小长度).
下面的代码使用了这个后缀树实现:http://www.daimi.au.dk/~mailund/suffix_tree.html.它是用C语言编写的,带有Python绑定.它不会处理大量的序列,但是两个没有问题.但是我不能说这是否会处理长度10,000,000.
from suffix_tree import GeneralisedSuffixTree
baseComplement = { 'A' : 'T', 'C' : 'G', 'G' : 'C', 'T' : 'A' }
def revc(seq):
return "".join([baseComplement[base] for base in seq[::-1]])
data = "AGGGTTTCCCTGACCTTCACTGCAGGTCATGCA"
# revc TGCATGACCTGCAGTGAAGGTCAGGGAAACCCT
# 012345678901234567890123456789012
# 1 2 3
minlength = 6
n = len(data)
tree = GeneralisedSuffixTree([data, revc(data)])
for shared in tree.sharedSubstrings(minlength):
#print shared
_, start, stop = shared[0]
seq = data[start:stop]
_, rstart, rstop = shared[1]
rseq = data[n-rstop:n-rstart]
print "Match: {0} at [{1}:{2}] and {3} at [{4}:{5}]".format(seq, start, stop, rseq, n-rstop, n-rstart)
这会产生输出
Match: AGGTCA at [23:29] and TGACCT at [10:16]
Match: TGACCT at [10:16] and AGGTCA at [23:29]
Match: CTGCAG at [19:25] and CTGCAG at [19:25]
它每次匹配两次,每次一次.那里也有反向回文!