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

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

第一 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 ,with必须比我的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')

在仔细考虑了这个问题之后,我意识到,除非你从整个字符串开始,并不断地删掉结尾,直到找到回文或没有更多的字符串,否则你真的不知道你是否找到了最长的回文。

# 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文件的多处理版本中使用6个工作进程运行了这个脚本( 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