使用具有多个规则的pyparsing-Forward()来查找回文序列

Question

我正试图通过递归在DNA中找到回文序列。我这么做是因为不可能知道DNA中回文序列的确切长度。我已经在脑子里和纸上解决了这个问题，但下面的代码仍然没有找到我想要的答案。因此，我设置代码的方式可能有误。欢迎任何帮助。

stem = Forward()
atRule = Word("A") + ZeroOrMore(stem) + Word("T")
taRule = Word("T") + ZeroOrMore(stem) + Word("A")
gcRule = Word("G") + ZeroOrMore(stem) + Word("C")
cgRule = Word("C") + ZeroOrMore(stem) + Word("G")
stem << locatedExpr(Combine(atRule + taRule + gcRule + cgRule))
print(stem.parseString("AAAGGGCCCTTTAAAGGGCCCTTT"))

Answer 1

我想您对一些事情感到困惑。即便如此，我也无法让任何解析器匹配尽可能长的回文，我认为这就是你的目标。

首先，Word("A")匹配一个或多个A。类似地，Word("T")匹配一个或多个T。所以：AAAAT将被匹配为回文。相反，让我们做Literal("A") + ... + Literal("T")

其次，ZeroOrMore(stem)意味着你可以有多个内部回文。这将匹配："A AT TA T"，这不是回文。相反，让我们执行Optional(stem)。

第三，+运算符表示串联，而不是交替。atRule + taRule + gcRule + cgRule的意思是"一个AT回文，后面是TA回文，然后是GC回文，再后面是CG回文。"相反，让我们使用|。

第四，您调用locatedExpr，它必须比我的pyparsing副本更新。我已经包含了它，并且我稍微改变了它的用途。

这是修改后的程序：

from pyparsing import *
def locatedExpr(expr):
    locator = Empty().setParseAction(lambda s,l,t: l)
    return Group(locator("locn_start") + expr("value") + locator.copy().leaveWhitespace()("locn_end"))
stem = Forward()
atRule = Literal("A") + Optional(stem) + Literal("T")
taRule = Literal("T") + Optional(stem) + Literal("A")
gcRule = Literal("G") + Optional(stem) + Literal("C")
cgRule = Literal("C") + Optional(stem) + Literal("G")
stem << Combine(atRule | taRule | gcRule | cgRule)
lstem = locatedExpr(stem)
print(lstem.parseString('AT'))
print(lstem.parseString('ATAT'))
print(lstem.parseString("AAAGGGCCCTTTAAAGGGCCCTTT"))

结果如下：

[[0, 'AT', 2]]
[[0, 'AT', 2]]
[[0, 'AAAGGGCCCTTT', 12]]

注意，结果是最小的初始回文，而不是整个字符串。虽然我不认为这是你的目标，但我希望我的改变能让你更接近。

编辑：

如果你的目标是确定一个字符串是否是回文（与"在更大的字符串中搜索回文"相比），那么这个程序可能更容易使用：

def DNA_complement(s):
    d = {'A':'T', 'T':'A', 'C':'G', 'G':'C'}
    return ''.join(d.get(ch,'') for ch in s)
def DNA_reversed_complement(s):
    return DNA_complement(reversed(s))
def DNA_palindrome(s):
    return s == DNA_reversed_complement(s)
print DNA_palindrome('AT')
print DNA_palindrome('ATAT')
print DNA_palindrome('AAAGGGCCCTTTAAAGGGCCCTTT')

Answer 2

在仔细思考这个问题后，我意识到你真的不知道你是否找到了最长的回文，除非你从整个字符串开始，并不断剪掉结尾，直到你找到一个回文或没有剩下的字符串。

# sample is a bit longer then the original, I added
# some other characters to the beginning of the string
sample = "AATTAAAAAGGGCCCTTTAAAGGGCCCTTTAAAGGGCCCTTT"
nucleotide_map = {'A':'T', 'T':'A', 'G':'C', 'C':'G'}
# simple function to test if a string is a palindrome
def is_palindrome(st, start, end):
    # unlike conventional palindromes, genetic 
    # palindromes must be of even length
    if (end <= start) or ((end-start) % 2 == 0):
        return False
    s,e = start,end
    while e > s:
        if nucleotide_map[st[s]] != st[e]:
            return False
        s += 1
        e -= 1
    return True

def find_longest_palindrome(st, loc):
    s,e = loc, len(st)-1
    while e > s:
        if is_palindrome(st, s, e):
            return st[s:e+1], e+1
        e -= 1
    return '',loc

现在在sample中找到回文只是：

loc = 0
FIND_OVERLAPPING = False
while loc <= len(sample):
    pal, tmploc = find_longest_palindrome(sample, loc)
    if pal:
        print(pal, loc, len(pal))
        # advance to next location to look for palindromes
        if FIND_OVERLAPPING:
            loc += (tmploc-loc)/2
        else:
            loc = tmploc
    else:
        loc += 1

我不确定我是否对寻找重叠的回文有正确的看法。不要只想增加loc，否则你只会得到所有退化的情况-也就是说，对于AAGAATTCTT的回文，你还会得到AGATTCT、GAATTC、AATT和AT。

以下是一种将其缝合到pyparsing解析器中的方法：

from pyparsing import *
def getAllPalindromes(s,loc,t):
    FIND_OVERLAPPING = False
    ret = []
    sample = t[0]
    while loc <= len(sample):
        pal, tmploc = find_longest_palindrome(sample, loc)
        if pal:
            ret.append((pal, loc))
            # advance to next location to look for palindromes
            if FIND_OVERLAPPING:
                loc += (tmploc-loc)/2
            else:
                loc = tmploc
        else:
            loc += 1
    return ret
palin = Word('AGCT').setParseAction(getAllPalindromes)
palin.parseString(sample).pprint()

给出结果：

[('AATT', 0), ('AAAGGGCCCTTTAAAGGGCCCTTTAAAGGGCCCTTT', 6)]

编辑

我只是在这个脚本的多处理版本中使用6个工作进程来运行这个FASTA文件(http://toxodb.org/common/downloads/Current_Release/Tgondii/fasta/ToxoDB-28_Tgondii_ESTs.fasta)，在首先处理FASTA格式以将每个封装的多行序列转换为单个字符串之后。在147151个序列中，该脚本在22分钟内发现了超过1000个长度为20个或更多核苷酸的回文。以下是一些示例回文：

ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAT
CATATATATATATATATATATATATATATATATATATATATATATATATATATATATATG
CCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGG
CCTCGTGCCGAATTCGGCACGAGGCCTCGTGCCGAATTCGGCACGAGG
CTCGTGCCGAATTCGGCACGAGCTCGTGCCGAATTCGGCACGAG
AAAAAAAAAAAAAAAAAACTCGAGTTTTTTTTTTTTTTTTTT
GGCACGAGGCCTCGTGCCGAATTCGGCACGAGGCCTCGTGCC
TTTATATATAAATATTTATATATAAATATTTATATATAAA