如何在FASTA序列中找到反向重复模式？

Question

假设我的长序列看起来像：

5’-AGGGTTTCCC**TGACCT**TCACTGC**AGGTCA**TGCA-3

在这个长序列中的两个斜体子序列（在这两个星星内）一起被称为反转重复模式。这两个子序列中的四个字母（如A、T、G、C）的长度和组合将是不同的。但这两个子序列之间有联系。请注意，当您考虑第一个子序列时，它的互补子序列是ACTGGA（根据A与T的组合，G与C的组合），当您反转此互补子序列（即最后一个字母排在第一位）时，它与第二个子序列相匹配。在

在FASTA序列中有大量这样的模式（包含1000万个ATGC字母），我想找到这些模式及其起始和结束位置。在

Answer 1

这是一个强力实现，尽管它对超长序列可能不是很有用。在

def substrings(s, lmin, lmax):
    for i in range(len(s)):
        for l in range(lmin, lmax+1):
            subst = s[i:i+l]
            if len(subst) == l:
                yield i, l, subst

def ivp(s, lmin, lmax):
    mapping = {'A': 'T', 'G': 'C', 'T': 'A', 'C': 'G'}
    for i, l, sub in substrings(s, lmin, lmax):
        try:
            from string import maketrans
        except ImportError: # we're on Python 3
            condition = sub.translate(
                       {ord(k): v for k, v in mapping.items()})[::-1] in s
        else: # Python 2
            condition = sub.translate(maketrans('ATGC', 'TACG'))[::-1] in s
        if condition:
            yield i, l, sub

让我们找出长度为6的“反向重复模式”（及其起始位置和长度）：

^{pr2}$

不过，这并不能检查这两种模式是否重叠。例如，'CTGCAG'匹配自身。在

Answer 2

我对Python和生物信息学都是新手，但我正在努力通过罗莎琳德信息在网站上学习一些两者。你用后缀树来做这个。后缀树（见http://en.wikipedia.org/wiki/Suffix_tree）是一种神奇的数据结构，它使生物信息学中的一切成为可能。您可以快速找到多个长序列中的公共子字符串。后缀树只需要线性时间，所以长度10000000是可行的。在

所以首先找到序列的反向补码。然后把它们放到后缀树中，找到它们之间的公共子串（长度最小）。在

下面的代码使用这个后缀树实现：http://www.daimi.au.dk/~mailund/suffix_tree.html。它是用C语言编写的，带有Python绑定。它不能处理大量的序列，但是两个不是问题。不过，我不能说这是否能处理长度10000000。在

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)

这就产生了输出

^{pr2}$

它找到每个匹配项两次，每端一次。还有一个反向回文！在

Answer 3

我有一个想法来找到一个长序列的倒转回文序列。考虑整个序列的DNA序列的一部分并生成其补体。然后反转这个补码序列的部分。然后对这两部分进行动态比对，并计算其代价（一部分来自实际序列，另一部分来自反向补码序列）。成本将提供一个想法，哪一个路线是最好的。现在，如果最佳对齐成本=阈值成本，则选择该部分并找到公共区域。这个特定部分的两个公共区域将是一个反向重复单元。一旦找到单位，然后在下一个部分重复它，其他方式继续增加部分的长度。有人能实现这个算法吗。也许这将是一个有用的解决办法。在