2017-08-21 6 views
0

Ich habe in den letzten paar Tagen mit einem bestimmten Bioinformatikproblem zu kämpfen, und ich fragte mich, ob jemand meine Logik oder Fehler in meinem Code finden könnte (oder beides). Die Funktion soll alle k-mers mit einer Hamming-Distanz von höchstens d von allen DNS-Strängen finden. Das ist, was ich versuche zu tun:DNA-Motivaufzählung mit Try/Except und Loops - Python3

  1. Beginnen Sie mit einer Iteration aller möglichen k-mers und vergleichen sie

  2. zu jeder Saite Dies bedeutet, ich eine andere Schleife müssen, die durch die Saiten geht von DNA auch

  3. Ich mache eine Weile Schleife für c+k <= len(DNA[0])-1. c+k ist mein Fenster der Größe k, und ich möchte mindestens ein Fenster in jeder DNA-Kette finden, wo meine Combo eine Hamming-Distanz von dieser Zeichenkette gleich oder kleiner als eine beliebige d hat. Wenn die Hamming-Distanz das Kriterium erfüllt, bricht die while-Schleife ab und erlaubt den nächsten String zu vergleichen. Ist dies nicht der Fall, wird das Fenster geändert, und wenn c+k==len(DNA[0])-1, und die Hamming-Distanz immer noch nicht das Kriterium erfüllt, erstelle ich einen Namensfehler int(a) und die Ausnahme tötet die inner_loop.

aber meine Funktion gibt nichts außer set(), die ich nicht verstehe.

import itertools 
def combination(k): 
    bases=['A','T','G','C'] 
    combo=[''.join(p) for p in itertools.product(bases, repeat=k)] 
    return combo 

def hammingDistance(Pattern, seq): 
     if Pattern == seq: 
       return 0 
     else: 
       dist=0 
       for i in range(len(seq)): 
         if Pattern[i] != seq[i]: 
           dist += 1 
     return dist 

def motif_enumeration(k, d, DNA): 
    combos = combination(k) 
    global pattern 
    for combo in combos: 
     try: 
      inner_loop(k, d, DNA, combo) 
     except: 
      continue 

    return set(pattern) 

def inner_loop(k, d, DNA, combo): 
    global pattern 
    for strings in DNA: 
     inner_loop_two(k, d, DNA, combo, strings) 


def inner_loop_two(k, d, DNA, combo, strings): 
    global pattern 
    c=0 
    while c+k < len(DNA[0]): 
     print(combo, strings[c:c+k], hammingDistance(combo, strings[c:c+k])) 
     if d >= hammingDistance(combo, strings[c:c+k]) and strings == DNA[len(DNA)-1]: 
      #if we've reached the last string and the condition is met, 
      #that means that the combo is suitable for each string of DNA 
      pattern += [combo] 
     elif d >= hammingDistance(combo, strings[c:c+k]): 
      #condition is met for one string, now move onto next 
      break 
     elif d < hammingDistance(combo, strings[c:c+k]) and c+k == len(DNA[0])-1: 
      #Name error causes this inner loop two to crash, thus causing the first inner loop 
      #to pass 
      int(a) 
     elif d < hammingDistance(combo, strings[c:c+k]): 
      #change the window to see if the combo is valid later in the string 
      c += 1 


pattern = [] 
DNA=['ATTTGGC', 'TGCCTTA', 'CGGTATC', 'GAAAATT'] 
print(motif_enumeration(3,1,DNA)) 
print(pattern) 

Ich dachte, dass meine zweite innere Schleife abgestürzt, da dies wäre meine erste innere Schleife verursachen, passieren und dann wäre eine andere Combo in motif_enumeration getestet werden, aber die erste Bedingung in meinem inner_loop_two nie etwas druckt. Ich bemerkte auch, dass, wenn die innere Schleife abstürzt und motif_enumeration weitergeht, es für die äußere und die innere Schleife weitergeht. Hier ist ein Beispiel dafür, was ich meine ...

AAA ATT 2 
AAA TTT 3 
AAA TTG 3 
AAA TGG 3 
AAT ATT 1 
AAT TGC 3 
AAT GCC 3 
AAT CCT 2 
AAT CTT 2 
AAG ATT 2 
AAG TTT 3 
AAG TTG 2 
AAG TGG 2 
AAC ATT 2 
AAC TTT 3 
AAC TTG 3 
AAC TGG 3 
ATA ATT 1 etc... 

Meine erwartete Ausgabe ist pattern=[ATA, ATT, GTT, TTT]

+0

Wer ist nur hit-and-run meine Posts downvoting, wenn Sie das tun, können Sie zumindest versuchen, mir zu helfen und mir sagen, warum meine Frage dumm ist? Ich frage nur, weil ich viel Zeit in dieses Problem gesteckt habe und ehrlich gesagt nicht weiß, wie ich meinen Code reparieren soll. – DrJessop

+1

Ich habe gerade Ihren Code geladen und festgestellt, dass ich keine Ahnung habe, wie ich das anfangen soll. Sie rufen niemals eine Funktion auf, um diese Funktion zu starten. Wie laufen wir das eigentlich? – roganjosh

+0

@roganjosh Ich habe gerade meinen Post bearbeitet. Es ist am Ende – DrJessop

Antwort

2

Die Kernkomponente der Logik ist, dass wir eine Combo in das Muster gesetzt sammeln möchten, wenn die Combo bei Spiele alle Position auf alle der Zielzeichenfolgen. Wir können Pythons all und any Funktionen dazu verwenden. Diese Funktionen arbeiten effizient, weil sie den Test beenden, sobald das Ergebnis entschieden ist.

import itertools 

def combination(k): 
    return (''.join(p) for p in itertools.product('ATCG', repeat=k)) 

def hamming_distance(pattern, seq): 
    return sum(c1 != c2 for c1, c2 in zip(pattern, seq)) 

def window(s, k): 
    for i in range(1 + len(s) - k): 
     yield s[i:i+k] 

def motif_enumeration(k, d, DNA): 
    pattern = set() 
    for combo in combination(k): 
     if all(any(hamming_distance(combo, pat) <= d 
       for pat in window(string, k)) for string in DNA): 
      pattern.add(combo) 
    return pattern 

DNA = ['ATTTGGC', 'TGCCTTA', 'CGGTATC', 'GAAAATT'] 
print(motif_enumeration(3, 1, DNA)) 

Ausgang

{'GTT', 'ATA', 'TTT', 'ATT'} 

Ich habe ein paar andere Änderungen an Ihrem Code. Wir können die Hamming-Distanz effizient berechnen, indem wir einen Generator an die sum-Funktion übergeben. Und wir können Zeit sparen & RAM durch die Umwandlung der Combo Tupel zu Strings mit einem Generator, anstatt sie in eine Liste zu setzen.


Die motif_enumeration Funktion kann noch weiter in eine Menge Verständnis kondensiert sein, aber ich muss zugeben, dass es ziemlich dicht ist, und noch schwieriger als die vorherige Version zu lesen. Es kann jedoch etwas effizienter sein.

def motif_enumeration(k, d, DNA): 
    return {combo for combo in combination(k) 
     if all(any(hamming_distance(combo, pat) <= d 
      for pat in window(string, k)) for string in DNA)} 

Und hier ist ein etwas besser lesbare Version, wo ich motif_enumeration eine Hilfsfunktion in_window auszuführen, um den inneren Test gegeben habe.

+0

Ich verstehe jetzt. Ich habe es verpasst, allen Strings, die der Funktion zugewiesen wurden, gemeinsam zu sein. – roganjosh

+0

Dies ist eine sehr elegante Antwort. Wenn ich vorher schon alles gewusst hätte ... Danke – DrJessop

+0

@DrJessop Mein Vergnügen! Ich habe eine leicht lesbare Version am Ende meiner Antwort hinzugefügt. Generatoren/Comprehensions können zu den besten Zeiten schwer zu lesen sein, und ihre Verschachtelung macht es viel schwieriger. ;) –