<p>我对Python和生物信息学都是新手,但我正在努力通过罗莎琳德信息在网站上学习一些两者。你用后缀树来做这个。后缀树(见<a href="http://en.wikipedia.org/wiki/Suffix_tree" rel="nofollow noreferrer">http://en.wikipedia.org/wiki/Suffix_tree</a>)是一种神奇的数据结构,它使生物信息学中的一切成为可能。您可以快速找到多个长序列中的公共子字符串。后缀树只需要线性时间,所以长度10000000是可行的。在</p>
<p>所以首先找到序列的反向补码。然后把它们放到后缀树中,找到它们之间的公共子串(长度最小)。在</p>
<p>下面的代码使用这个后缀树实现:<a href="http://www.daimi.au.dk/~mailund/suffix_tree.html" rel="nofollow noreferrer">http://www.daimi.au.dk/~mailund/suffix_tree.html</a>。它是用C语言编写的,带有Python绑定。它不能处理大量的序列,但是两个不是问题。不过,我不能说这是否能处理长度10000000。在</p>
<pre><code>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)
</code></pre>
<p>这就产生了输出</p>
^{pr2}$
<p>它找到每个匹配项两次,每端一次。还有一个反向回文!在</p>