Использование pyparsing Forward() с несколькими правилами для поиска палиндромной последовательности

Question

Я пытаюсь найти палиндромные последовательности в ДНК с помощью рекурсии. Я делаю это, потому что невозможно узнать точную длину палиндромной последовательности в ДНК. Я справился с этой проблемой в своей голове и на бумаге, но код ниже все еще не приносит ответа, который я хочу. Я новичок в 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 означает «палиндром AT, за которым следует палиндром TA, а затем палиндром GC, а затем палиндром CG». Вместо этого давайте использовать |.

В-четвертых, вы вызываете locatedExpr, должно быть новее, чем моя копия пирапса. Я включил его, и я немного изменил его использование.

Вот модифицированная программа:

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

Обратите внимание, что результат является минимальной начальной палиндром, а не вся строка. Хотя я не думаю, что это ваша цель, я надеюсь, что мои изменения приблизит вас.

EDIT:

Если ваша цель состоит в том, чтобы определить, является ли строка палиндромом (контрастирует с «поиск палиндром в большей строки»), то эта программа может быть намного проще в использовании:

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)] 

EDIT

Я просто побежал это с помощью 6 рабочих процессов в многопроцессорную версию этого скрипта в отношении этого файла FASTA (http://toxodb.org/common/downloads/Current_Release/Tgondii/fasta/ToxoDB-28_Tgondii_ESTs.fasta) после первой обработки формы FASTA t для преобразования каждой завернутой многострочной последовательности в одну строку. В 147151 последовательностях сценарий обнаружил более 1000 палиндромов длиной 20 или более нуклеотидов в длину за 22 минуты. Вот некоторые из палиндромов образца:

ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAT 
CATATATATATATATATATATATATATATATATATATATATATATATATATATATATATG 
CCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGG 
CCTCGTGCCGAATTCGGCACGAGGCCTCGTGCCGAATTCGGCACGAGG 
CTCGTGCCGAATTCGGCACGAGCTCGTGCCGAATTCGGCACGAG 
AAAAAAAAAAAAAAAAAACTCGAGTTTTTTTTTTTTTTTTTT 
GGCACGAGGCCTCGTGCCGAATTCGGCACGAGGCCTCGTGCC 
TTTATATATAAATATTTATATATAAATATTTATATATAAA