2016-04-08 13 views
1

回帰を使用してDNA内の回文配列を検索しようとしています。 DNA内のパリンドローム配列の正確な長さを知ることは不可能なので、私はこれをやっています。私は頭と紙の問題を解決しましたが、下のコードは、私が望む答えがまだ出ていません。私は、CFGをpyparsingおよび使用するのが初めてです。したがって、私はコードを設定する方法に間違っている可能性があります。どんな助けも歓迎されるでしょう。回文を使用して回文配列を見つけるための複数の規則でForward()を使用する

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")) 
+0

私はこのスタイルのパリンドロームに精通していません - DNA分析は異なる定義を使用していますか?他のフィールドでは、回文は、文字列全体が逆の場合に文字が同じである単語または文です。あなたのサンプルでは、​​唯一の回文は、 "AAA"、 "GGG"、 "CCC"などのグループであり、実際に回文が進むにつれてかなり退化しています。 (私は、 "AAGAACCCCGGTTGGCCCCAAGAA"のようなものを期待していました。これは、同じ前方と後方を読み取ります。)また、なぜあなたは 'locatedExpr'を使用していますか?一致するロジックは変更されず、マッチの開始/終了位置にのみ注釈が付けられます。 – PaulMcG

+1

"*遺伝学の文脈におけるパリンドロームの意味は、単語や文章で使用される定義とは若干異なります。" - https://en.wikipedia.org/wiki/Palindromic_sequence –

+0

@user、あなたのサンプルコードは不完全です。エラーを示す短い、完全な、スタンドアロンのプログラムを制作できますか?詳細は、[mcve]を参照してください。 –

答えて

2

あなたはいくつかのことについて混乱していると思います。それでも、私はあなたの目標だと思う、最も長い可能性のあるパリンドロームに合ったパーサを手に入れることはできません。

まず、Word("A")は、1つ以上のAと一致します。同様にWord("T")は1つ以上のTと一致します。だから:AAAATは回文として一致します。代わりに、やりましょうLiteral("A") + ... + Literal("T")

次に、ZeroOrMore(stem)は、複数の内部のパリンドロームを持つことができることを意味します。それは一致します: "A AT TA T"、これは回文ではありません。代わりに、Optional(stem)をしましょう。

第3に、+演算子は、連結ではなく、連結を表します。 atRule + taRule + gcRule + cgRuleは、「ATパリンドロームに続いてTAパリンドローム、次にGCパリンドローム、続いてCGパリンドローム」を意味する。代わりに、|を使用しましょう。

第4に、locatedExprと呼んでください。これは、私のpyparsingのコピーよりも新しいものである必要があります。私はそれを含めて、私はその使用をわずかに変更しました。ここで

は修正プログラムです:結果は、最小限の初期回文、および全体ではなく文字列であることを

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

は注意:ここで

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

結果です。私はこれがあなたの目標だとは思っていませんが、私の変更があなたをより近づけることを願っています。

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') 
+0

あなたの変更によって私は最終的な答えに近づき、文字列全体に一致させたいと言って書いています。どうもありがとうございました。 defExpr(expr)が何をしているのかを私に説明することができたら、それはすばらしいでしょう。私は、私はこれらのクラスとメソッドのいくつかを使用していないので、pyparsingの新しいです。 – esseldm

+0

pyparsingモジュールの私のコピー(Ubuntu 14.04のpyparsing 2.0.1)には公式バージョンの 'locatedExpr'がありませんでしたので、私は[元のソース](https://pythonhosted.org/pyparsing/pyparsing)からそれを盗まなければなりませんでした-pysrc.html#locatedExpr)。 'pyparsing'モジュールから既に' locatedExpr'を取得しているので、私の定義を無視することができます。 –

+0

私はあなたが何を意味するかを見ます。再度、感謝します。 – esseldm

1

もう少し考えてみたら、文字列全体から始まっていない限り、最長の回文を見つけたかどうか分からず、回文を見つけるか、もう文字列が残っていない。今sampleに回文を見つける

# 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 

だけです:

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

は結果を与える:私はこの使用して6ワーカープロセスを実行した

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

EDIT

をこのFASTAファイル(http://toxodb.org/common/downloads/Current_Release/Tgondii/fasta/ToxoDB-28_Tgondii_ESTs.fasta)に対するこのスクリプトのマルチプロセッシング・バージョンは、最初にFASTA formaを処理した後各ラップされた複数行のシーケンスを単一の文字列に変換します。 147151配列では、スクリプトは、20分以上の長さのパリンドロームを22分で検出した。

ATATATATATATATATATATATATATATATATATATATATATATATATATATATATATAT 
CATATATATATATATATATATATATATATATATATATATATATATATATATATATATATG 
CCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGG 
CCTCGTGCCGAATTCGGCACGAGGCCTCGTGCCGAATTCGGCACGAGG 
CTCGTGCCGAATTCGGCACGAGCTCGTGCCGAATTCGGCACGAG 
AAAAAAAAAAAAAAAAAACTCGAGTTTTTTTTTTTTTTTTTT 
GGCACGAGGCCTCGTGCCGAATTCGGCACGAGGCCTCGTGCC 
TTTATATATAAATATTTATATATAAATATTTATATATAAA 
関連する問題