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

Question

我正在尝试使用递归查找 DNA 中的回文序列。我这样做是因为不可能知道 DNA 中回文序列的确切长度。我已经在头脑中和纸上解决了这个问题，但下面的代码仍然没有给出我想要的答案。我是 pyparsing 和使用 CFG 的新手。因此，我设置代码的方式可能有误。欢迎任何帮助。

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 表示 "an AT palindrome followed by a TA palindrome, followed by a GC palindrome, followed by a CG palindrome." 相反，让我们使用 |.

第四，你调用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]]

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

编辑：

如果你的目标是判断一个字符串是否为回文（与"search for a palindrome in a larger string"对比），那么这个程序可能更容易使用：

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 的回文，你还会得到 AGAATTCT、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)]

编辑

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

ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAT
CATATATATATATATATATATATATATATATATATATATATATATATATATATATATATG
CCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGG
CCTCGTGCCGAATTCGGCACGAGGCCTCGTGCCGAATTCGGCACGAGG
CTCGTGCCGAATTCGGCACGAGCTCGTGCCGAATTCGGCACGAG
AAAAAAAAAAAAAAAAAACTCGAGTTTTTTTTTTTTTTTTTT
GGCACGAGGCCTCGTGCCGAATTCGGCACGAGGCCTCGTGCC
TTTATATATAAATATTTATATATAAATATTTATATATAAA

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

Using pyparsing Forward() with multiple rules to find palindromic sequence

python

pyparsing

context-free-grammar