Needleman-Wunsch-Algorithmus

Python bietet eine Grundausstattung von Funktionen, sprich, sobald man es installiert hat, versteht es eine Reihe bestimmter Befehle. Es gibt jedoch speziellere Befehle, die man in Form so genannter Module zusätzlich zur Grundausrüstung aktivieren kann.
Um den Needleman-Wunsch-Algorithmus in Python umzusetzen, benötigt man das "Numeric"-Modul, denn dieses enthält die benötigten Matrizen. Das Einbinden dieser Matrixklassen an sich hat also noch nichts mit dem Needleman-Wunsch-Algorithumus zu tun, sondern schafft lediglich die erforderlichen Voraussetzungen.

import Numeric

Als nächstes wird die Distanzfunktion festgelegt. Die beiden Variablen, von denen diese abhängt, repräsentieren im Zuge der Berechnung jeden einzelnen Buchstaben ihrer jeweiligen Sequenzen.
Es wird festgelegt, dass die Alinierung zweier ungleicher Zeichen mit einem Strafpunkt belegt wird, die Alinierung gleicher Zeichen hingegen straffrei bleibt.

def d(a,b):

   if a == b
     return 0
  else:
     return 1


Die beiden Sequenzen, die aliniert werden sollen, bekommen nun einen Namen. Damit sind sie zum einen problemlos abrufbar (in diesem Fall z.B. mit print A), zum anderen wird auch der Bezug zu den Variablen klar: Da Zeichenfolgen mit Großbuchstaben, einzelne Zeichen hingegen mit Kleinbuchstaben benannt werden, ist klar, dass sämtliche Buchstaben a dem String A entstammen usw.

A = 'AGC'
B = 'ACGT'

(Das einfache Gleichheitszeichen = meint hier nicht die mathematische Gleichheit, sondern die Zuweisung der Variablen zu einem String. Da Sequenzen "Inhalt" sind, also nicht wie Variablen, Funktionen etc. als "Werkzeug und Zubehör" dienen, werden sie in Anführungszeichen gsetzt.)
Für den Algorithmus ist es wichtig, auch die Längen der Sequenzen zu kennen, denn sie bestimmen die Maße der Matrix. Die Matrix enthält dabei allerdings jedoch je eine Spalte und Zeile mehr als die Strings Zeichen, da jeweils ein Platz für das Gap-Zeichen benötigt wird.

m = len(A)
n = len(B)

Jetzt werden die beiden Matrizen ins Leben gerufen: Zum einen die Distanzmatrix ("D"), in der der direkte Vergleich zwischen den Zeichen und die Bewertung der zugehörigen Schritte stattfindet, zum zweiten diejenige, die die Schritte in der anderen verfolgt und bewertet. Man nennt sie Traceback-Matrix ("T"), sinngemäß also "Spurenzurückverfolgungsmatrix".

D = Numeric.zeros([m+1, n+1])
T = Numeric.zeros([m+1, n+1])

Aus dem Modul "Numeric zeros" ruft man nun die Matrizen auf; m+1 und n+1 sind die Anzahl ihrer Zeilen bzw. Spalten. Python nummeriert die Spalten und Zeilen automatisch durch, beginnt allerdings mit dem Zählen bei Null. Auf dieser Grundlage kann man die einzelnen Felder der Matrix exakt benennen, was abstrakt gehalten ungefähr so aussieht:



Anhand dieser Nummerierung erteilt man nun die Anweisung, wie die erste Zeile und die erste Spalte der Distanz-Matrix befüllt werden sollen. Geht man nur horizontal die erste Zeile entlang, so vergleicht man im ersten Feld das Gapzeichen mit dem Gapzeichen, was eine Distanz von 0 bedeutet. Ab dann jedoch wird die Distanz immer 1 sein, denn keine DNA- oder Proteinsequenz beinhaltet in ihrem Alphabet eine Lücke. Mit jedem Schritt, den man diese Zeile weiter geht, erhält man also einen weiteren Strafpunkt.

for i in range(m+1):
  D[i, 0] = i * d("-", "A")
for j in range(n+1):
  D[0, j] = j * d("-", "A")


Für die Felder der ersten Spalte gilt natürlich das Entsprechende.

Die Traceback-Matrix T dokumentiert den Weg, den man durch die Distanzmatrix nimmt, indem sie "0" für jeden diagonalen, "-1" für jeden horizontalen, und "1" für jeden senkrecht nach unten führenden Schritt speichert.
Wiederum können die erste Zeile und die erste Spalte statisch befüllt werden:

for i in range(m+1):
  T[i,0] = -1
for j in range(n+1):
  ;T[0,j] = 1

Jetzt kommt der entscheidende Punkt: Die dynamische Programmierung.
Indem man für jeden Schritt automatisch die Strafpunkte berechnet, die durch ihn und dadurch, dass er überhaupt möglich wird (Strafpunkte des Ausgangsfeldes werden jeweils zur neu berechneten Distanz addiert), zusammen kommen, muss man den Wert für jedes Feld nur einmal berechnen - die Rechenzeit wächst also proportional zur Sequenzlänge. (vgl.: O-Notation)
Nun führen ja zu jedem Feld drei Wege - von oben, diagonal von links/oben und von links. Um die Strafpunkte, die man für eine bestimmte Alinierung (also für das Einbeziehen des zugehörigen Feldes in den Weg durch die Matrix) erhält, zu berechnen, berechnet man zunächst die Distanz zwischen den beiden zugehörigen Zeichen. Dann addiert man je die Strafpunkte der drei möglichen Ausgangsfelder, erhält also drei Werte, von denen man den kleinsten als denjenigen, der den besten Weg beschreibt, in die Matrix einträgt.

for i in range(1,m+1):
  for j in range(1,n+1):


  d_horiz = D[i-1,j] + d(A[i-1],"-")
  d_diag = D[i-1,j-1] + d(A[i-1],B[j-1])
  d_vert = D[i,j-1] + d("-", B[j-1])


D[i,j] = min([d_horiz, d_diag, d_vert])

if d_diag <= d_horiz and d_diag <= d_vert:
     T[i,j] = 0
  elif d_horiz <= d_vert:
    T[i,j] = -1
  else:
    T[i,j] = +1


Die alinierten, d.h. sinnvoll mit Lücken befüllten Strings entsprechen natürlich nicht mehr den Originalsequenzen, können also auch nicht mehr A und B bezeichnet werden. Üblich wäre z.B. eine Benennung als A' und B', was jedoch nicht funktioniert, weil Python ' und " gleichermaßen als Anführungszeichen erkennt und somit "Input" erwartet, das nachfolgende = folglich als Datenmaterial und nicht als Zuweisungsoperator auffassen würde.
Wir nennen die alinierten Sequenzen daher einfach X und Y:

X = ""
Y = ""


Wir gehen nun in die (zwischenzeitlich ja gefüllte) Tracebackmatrix und verfolgen den Weg zurück, bei dem am wenigsten Strafpunkte zusammen kommen. Startpunkt ist die rechte untere Ecke der Matrix, also das Feld [m,n], Ziel ist das linke obere Eck [0,0].
Der Weg dorthin wird durch die gespeicherten Zahlen bestimmt: Beträgt der Wert eines Feldes 1, so rückt man ein Feld senkrecht nach oben; das bedeutet, dass der Sequenz Y das nun erreichte Zeichen aus B, der Sequenz X hingegen ein Gapzeichen zugeordet wird, da man ja keinen neuen Buchstaben aus A erreicht hat. Für 0 führt der Schritt diagonal nach links oben - beide Strings erhalten ein Zeichen aus der jeweils ursprünglichen Sequenz.

i = m
j = n


while i > 0 or j > 0:
  if T[i,j] == 0:
    X = A[i-1] + X
    Y = B[j-1] + Y
    (i,j) = (i-1, j-1)


  elif T[i,j] == -1:
    X = A[i-1] + X
    Y = "-" + Y
    i = i - 1


  else:
    X = "-" + X
    Y = B[j-1] + Y
    j = j -1



Zuletzt erfolgt der Befehl, die beiden neuen Sequenzen auszugeben:

print X
print Y