我是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]

它每次匹配两次,每次一次.那里也有反向回文！