Autor
Norbert Blum, Uni Bonn
Matthias Kretschmer, Uni Bonn
Da vor 150 Jahren nahe bei Düsseldorf im Neandertal erstmals Skeletteile des so genannten Neandertalers gefunden wurden, ist 2006 das Jahr des Neandertalers. Seit diesem Fund stellte sich u.a. die Frage, inwiefern der Homo sapiens, von dem wir abstammen, und der Neandertaler miteinander verwandt sind. Lange hielt man es für möglich, dass wir vom Neandertaler abstammen. Mittlerweile haben Wissenschaftler herausgefunden, dass die Entwicklungslinie des Homo neanderthalensis sich vor ca. 315.000 Jahren von jener trennte, die schließlich zum Homo sapiens führte. Dies hat man aufgrund der Unterschiede im Erbgut des Homo sapiens und des Homo neanderthalensis festgestellt. Neue Verfahren ermöglichen signifikante Anteile des Erbguts aus den mehr als 30.000 Jahre alten Knochen zu extrahieren. Dieses erhält man in Form von DNA-Sequenzen. DNA-Sequenzen kann man sich als Baupläne für Tiere und Menschen vorstellen. Während der evolutionären Geschichte haben sich DNA-Sequenzen aufgrund von Mutationen geändert. Hat man DNA-Sequenzen verschiedener Spezies, dann kann man die Ähnlichkeit dieser Sequenzen mit Hilfe des Computers schätzen. Wir werden die Ähnlichkeit von zwei DNA-Sequenzen über die Distanz dieser DNA-Sequenzen ausdrücken. Umso ähnlicher zwei DNA-Sequenzen sind, umso geringer ist die Distanz zwischen diesen.
Wie kann man einen Algorithmus entwickeln, der die Distanz zweier
DNA-Sequenzen schätzt?
Mit dieser Frage wollen wir uns zunächst beschäftigen. Hierzu benötigen
wir als erstes ein mathematisches Modell. D.h., wir müssen DNA-Sequenzen und
Mutationen modellieren, um die Distanz zweier Sequenzen
formal definieren zu können.
Eine DNA-Sequenz ist eine Zeichenkette bestehend aus den Zeichen A, G, C
und T. Formal sagt man, eine DNA-Sequenz ist ein String x über dem
Alphabet
= {AGCT}. So ist zum Beispiel
CAGCGGAAGGTCACGGCCGGGCCTAGCGCCTCAGGGGTG |
ein Ausschnitt der DNA-Sequenz des Huhnes. Eine Mutation führt eine DNA-Sequenz x in eine DNA-Sequenz y über. Wir gehen davon aus, dass es nur drei verschiedene Mutationen (die wir auch Operationen nennen werden) gibt:
- Streichen eines Zeichens,
- Einfügen eines Zeichens und
- Ersetzen eines Zeichens durch ein anderes.
![$ \to$](/~algorithmus/Algorithmen/algo31/pictures/img2.png)
![$ \to$](/~algorithmus/Algorithmen/algo31/pictures/img2.png)
![$ \to$](/~algorithmus/Algorithmen/algo31/pictures/img2.png)
Um ein Maß für die Distanz
zweier DNA-Sequenzen zu erhalten ordnen wir jeder Mutation Kosten zu.
Diese Kosten haben etwas mit der Wahrscheinlichkeit, dass eine korrespondierende
Mutation auftritt, zu tun. Je wahrscheinlicher die Mutation umso geringer die
Kosten. Wir ordnen den drei verschiedenen Mutationen folgende Kosten zu:
- Streichen: Kosten 2
- Einfügen: Kosten 2
- Ersetzen: Kosten 3
![$ \to$](/~algorithmus/Algorithmen/algo31/pictures/img2.png)
![$ \to$](/~algorithmus/Algorithmen/algo31/pictures/img2.png)
![$\displaystyle \sum_{{i=1}}^{t}$](/~algorithmus/Algorithmen/algo31/pictures/img3.png)
Wir wollen die Kosten einer Mutationenfolge verwenden, um die Distanz zweier
DNA-Sequenzen zu definieren. Dabei müssen wir beachten, dass es mehrere
verschiedene Mutationenfolgen geben kann, die eine DNA-Sequenzen x nach einer
DNA-Sequenz y transformieren. Falls zum Beispiel x = AG und y = T, wie oben
ist, dann gibt es unter anderem folgende Möglichkeiten:
-
S1 = A
, G
T; c(S1) = c(A
) + c(G
T) = 2 + 3 = 5
-
S2 = A
T, G
; c(S2) = c(A
T) + c(G
) = 3 + 2 = 5
-
S3 = A
, G
,
T; c(S3) = c(A
) + c(G
) + c(
T) = 2 + 2 + 2 = 6
-
S4 = A
C, G
, C
,
T; c(S4) = c(A
C) + c(G
) + c(C
) + c(
T) = 3 + 2 + 2 + 2 = 9
Wie berechnet man nun die evolutionäre Distanz? D.h., wie berechnen wir
dc(x, y)? Wir haben nur die Operationen Streichen, Einfügen und Ersetzen zur
Verfügung, um x nach y zu transformieren. Die Operationen und das
Kostenmodell sind so definiert, dass eine Folge von mehreren Operationen an
derselben Position kostengünstiger durch eine einzelne Operation ersetzt werden
kann. So kann zum Beispiel des Streichen eines Zeichens a und das Einfügen
eines Zeichens b an derselben Position durch das Ersetzen von a durch b
kostengünstiger durchgeführt werden (Kosten 2 + 2 = 4 für eine Streiche- und eine
Einfüge-Operation gegenüber Kosten 3 für eine Ersetze-Operation). Die
Mutationen an den verschiedenen Positionen beeinflussen sich gegenseitig nicht.
Wir können diese Mutationen somit in beliebiger Reihenfolge durchführen. Wir
gehen nun davon aus, dass die letzte Mutation immer an den letzten Positionen
der beiden Sequenzen stattfindet.
Sei
x = a1a2...am und
y = b1b2...bn. D.h. x besteht aus m und
y aus n Zeichen. Gemäß unseren drei Operationen gibt es drei Möglichkeiten,
x nach y zu transformieren:
- Streichen: Wir streichen am und
a1a2...am-1 wird nach
b1b2...bn
transformiert.
- Einfügen: Wir fügen bn ein
und a1a2...am wird nach
b1b2...bn-1
transformiert.
- Ersetzen: Wir ersetzen am durch bn und a1a2...am-1 wird nach b1b2...bn-1 transformiert.
Wir werden jetzt obige Beobachtung zu einem Algorithmus ausarbeiten. Sei nun
x[1 : i] der Präfix von x, der aus den ersten i Zeichen besteht. Der Präfix
der Länge 0 von x ist der leere String x[1 : 0]. Der Präfix der Länge m
von x ist x selber, da x gerade m Zeichen lang ist. Entsprechend ist
y[1 : j] der Präfix von y, der aus den ersten j Zeichen besteht. Wir können
nun die obige Beobachtung wie folgt umformulieren:
- x[1 : m - 1] wird nach y transformiert und wir streichen am.
- x wird nach y[1 : n - 1] transformiert und wir fügen bn ein.
- x[1 : m - 1] wird nach y[1 : n - 1] transformiert und wir ersetzen am durch bn.
![\begin{displaymath}\begin{cases}d_c(x[1:i-1],y[1:j])+c(a_i\to), & \text{(Stre...
...[1:i-1],y[1:j-1])+c(a_i\to b_j) & \text{(Ersetzen)}
\end{cases}\end{displaymath}](/~algorithmus/Algorithmen/algo31/pictures/img4.png)
Beachte, dass für i = 0 und j > 0 nur eine Einfüge-Operation (d.h. Fall 2) und
für j = 0 und i > 0 nur eine Streiche-Operation (d.h. Fall 1) in Frage kommen.
Das rührt daher, dass wir in einem leeren String kein Zeichen streichen und auch
keines ersetzen können. Um den leeren String zu erhalten, können wir in einem
nichtleeren String nur Zeichen streichen. Im Fall i = 0 und j = 0 müssen wir den
leeren String x[1 : 0] nach den leeren String y[1 : 0] transformieren. Dafür
benötigen wir keine einzige Operation, da beide Strings identisch sind. Somit
ergibt sich
dc(x[1 : 0], y[1 : 0]) = 0.
Wähle zum Beispiel x = AGT und y = CAT. Nehmen wir an, wir hätten schon
dc(x[1 : 1], y[1 : 2]) = dc(A, CA),
dc(x[1 : 2], y[1 : 1]) = dc(AG, C) und
dc(x[1 : 1], y[1 : 1]) = dc(A, C) bestimmt. Dann können wir daraus die evolutionäre
Distanz von
dc(x[1 : 2], y[1 : 2]) = dc(AG, CA) nach obigem Muster bestimmen. Dazu
bestimmen wir das Minimum von
-
dc(x[1 : 1], y[1 : 2]) + c(a2
) = dc(A, CA) + c(G
) = dc(A, CA) + 2,
-
dc(x[1 : 2], y[1 : 1]) + c(
b2) = dc(AG, C) + c(
A) = dc(AG, C) + 2 und
-
dc(x[1 : 1], y[1 : 1]) + c(a2
b2) = dc(A, C) + c(G
A) = dc(A, C) + 3.
Wie können wir nun diese Überlegungen verwenden, um einen Algorithmus zu
formulieren? Die evolutionären Distanzen von Präfixen von x und y werden
mehrfach benötigt. So wird zum Beispiel
dc(x[1 : i - 1], y[1 : j - 1]) zur Berechnung
von
dc(x[1 : i], y[1 : j - 1]),
dc(x[1 : i - 1], y[1 : j]) und
dc(x[1 : i], y[1 : j])
benötigt. Wir möchten die Berechnung von
dc(x[1 : i - 1], y[1 : j - 1]) nur einmal
durchführen. Demzufolge müssen wir diesen Wert speichern, so dass wir ihn
mehrfach verwenden können. Wir verwenden hierzu eine Tabelle, in der wir die
schon berechneten evolutionären Distanzen zwischenspeichern. In dieser Tabelle
verwenden wir die Zelle (i, j) (Zeile i und Spalte j) zur Speicherung der
evolutionären Distanz
dc(x[1 : i], y[1 : j]). Der Vorteil dieser Methode ist, dass
es wesentlich effizienter ist, die Werte direkt aus der Tabelle zu lesen, als
sie jedesmal erneut zu berechnen. Wir benötigen die evolutionären Distanzen für
alle
0i
m und
0
j
n. Also besteht unsere Tabelle aus m + 1
Zeilen und aus n + 1 Spalten.
Sei zum Beispiel x = ATGAACG und y = TCAAT. Die dazugehörige Tabelle wäre, bei
der obigen Kostenfunktion (Streichen und Einfügen kosten 2, Ersetzen kostet
3):
Wir fangen bei der Berechnung der Distanzen mit dc(x[1 : 0], y[1 : 0]) an. Wir wissen schon, dass dc(x[1 : 0], y[1 : 0]) = 0 ist. Also notieren wir uns in der Zelle (0, 0) den Wert 0. Wenn wir die Spalte 0 von oben nach unten durchgehen, können wir nur Streiche-Operation durchführen. Das haben wir oben bereits bemerkt (j ist in diesem Fall gleich 0). Man kann nur Zeichen streichen, da wir x[1 : i] nach den leeren String y[1 : 0] transformieren wollen. Entsprechend können wir in der Zeile 0 nur Einfüge-Operationen durchführen.
Die anderen Tabelleneinträge berechnen sich nun nach der Formel von oben.
Betrachten wir zum Beispiel die Zelle (2, 1). Wir wollen den String x[1 : 2] = AT
nach
den String y[1 : 1] = T transformieren. Wir sehen direkt, dass es mit einer
Streiche-Operation erledigt ist. Wir müssen nur A entfernen und danach T
durch T ersetzen. Das ist auch
intuitiv die Lösung mit den minimalen Kosten. Der Algorithmus muss demnach auch
diese Distanz (eine Streiche-Operation:
Kosten = 2) berechnen. Dies ist
auch der Fall, denn er bildet das Minimum aus den drei Werten, die den drei
Operationen entsprechen:
-
dc(x[1 : 1], y[1 : 1]) + c(a2
) = dc(A, T) + c(a2
) = 3 + 2 = 5 (Streichen von T),
-
dc(x[1 : 2], y[1 : 0]) + c(
b1) = dc(AT, y[1 : 0]) + c(
b1) = 4 + 2 = 6 (Einfügen von T) und
-
dc(x[1 : 1], y[1 : 0]) + c(a2
b1) = dc(A, y[1 : 0]) + c(a2
b1) = 2 + 0 = 2 (Ersetzen von T durch T).
Die kleinen Striche innerhalb der Tabelle sollen verdeutlichen über welche
,,Wege`` die optimale Transformation von x[1 : i] nach y[1 : j] ablaufen kann.
In dem Fall der Zelle (2, 1) sind wir gerade diagonal von (1, 0) über eine
Ersetzen-Operation zu den minimalen Kosten gekommen. Diese Wege sind nicht
notwendigerweise eindeutig, wie man anhand der Tabelle sehen kann. Die Striche
können direkt bei dem Aufbau der Tabelle bestimmt werden, da sie nur angeben,
welcher der drei möglichen Schritte an den letzten Positionen der Präfixe, zu
minimalen Kosten führt.
Die rot markierten Wege, geben die ,,Wege minimaler Kosten`` für die
Transformation von x nach y an. D.h. sie geben die Wege an, über die man mit
minimalen Kosten x nach y transformieren und somit die evolutionäre
Distanz bestimmen kann.
Da wir nicht wissen, welche Tabelleneinträge zu einem ,,Weg minimaler Kosten``
für die Transformation von x nach y beitragen, müssen wir alle bestimmen.
Wir benötigen zur Berechnung des Eintrags von Zelle (i, j) die Werte der Zellen
(i - 1, j), (i, j - 1) und (i - 1, j - 1). Um sicherzustellen, dass diese
evolutionären Distanzen bereits berechnet sind, wenn wir
dc(x[1 : i], y[1 : j])
berechnen wollen, bauen wir die Tabelle zeilenweise von oben nach unten oder
spaltenweise von links nach rechts auf. Wenn wir die Tabelle zeilenweise
aufbauen, müssen wir jede Zeile von links nach rechts durchlaufen. Entsprechend
muss jede Spalte beim spaltenweisen Aufbau von oben nach unten durchlaufen
werden. Am Ende steht dann in Zelle (m, n) das gewünschte Ergebnis: die
evolutionäre Distanz
dc(x[1 : m], y[1 : n]) von x[1 : m] = x und y[1 : n] = y.
Das folgende Applet bietet eine Visualisierung des Verfahrens (das Applet wurde
von Ryan
McKinley erstellt und von Sascha
Meinert und Martin
Nöllenburg modifiziert).
Rekapitulieren wir: Wir haben, beginnend bei dem einfachsten und kleinsten
Teilproblem, die Berechnung von
dc(x[1 : 0], y[1 : 0]), immer größer werdende
Teilprobleme gelöst. Wir haben die Präfixe von x und y immer größer werden
lassen und für diese die evolutionäre Distanz berechnet. Dabei haben wir zur
Bestimmung von
dc(x[1 : i], y[1 : j]) die drei evolutionären Distanzen
dc(x[1 : i - 1], y[1 : j]),
dc(x[1 : i], y[1 : j - 1]) und
dc(x[1 : i - 1], y[1 : j - 1]) benötigt.
Abstrakt ausgedrückt: Wir haben optimale Teillösungen verwendet, um ein größeres
Teilproblem optimal zu lösen.
Die Frage ist, ob das ein allgemeines Prinzip ist, welches ein allgemeines
Programmierverfahren zulässt? Dies ist der Fall. Das Prinzip heißt
Optimalitätsprinzip und sieht folgendermaßen aus:
- Jede Teillösung einer optimalen Lösung, die Lösung eines Teilproblems ist, ist selbst eine optimale Lösung des betreffenden Teilproblems.
Das Programmierverfahren zur Lösung eines Optimierungsproblems, das auf dem Optimalitätsprinzip aufbaut, heißt dynamische Programmierung. Bei der dynamischen Programmierung teilen wir unser Problem in Teilprobleme auf. Die kleinsten Teilprobleme lösen wir direkt. In unserem Fall war dies gerade die Berechnung von dc(x[1 : 0], y[1 : 0]). Die Lösung hierfür war sehr einfach, da immer dc(x[1 : 0], y[1 : 0]) = 0 gilt. Aus den optimalen Lösungen zu kleinen Teilproblemen werden dann optimale Lösungen für größere Teilprobleme berechnet. Dies wird solange wiederholt, bis man eine optimale Lösung für das eigentliche Problem berechnet hat. Die dynamische Programmierung ist ein wichtiges allgemeines Verfahren, das beim Algorithmenentwurf häufig seine Anwendung findet.
Das hier vorgestellte Verfahren kann nicht nur zur Berechnung der evolutionären
Distanz zweier DNA-Sequenzen verwendet werden, sondern auch um die Ähnlichkeit
zweier Wörter der deutschen Sprache zu bestimmen. Dies wird unter anderem bei
Rechtschreibkorrekturprogrammen angewendet. Zum Beispiel wird dem Benutzer jedes
Wort, dass eine Distanz kleiner oder gleich 5 zu einem dem
Rechtschreibkorrekturprogramm unbekannten Wort hat, als Alternative für dieses
unbekannte Wort angeboten.
Wie man heute weiß, stammt der Homo sapiens nicht vom Neandertaler ab. Es wurden derartige Unterschiede im Erbgut des Neandertalers und des Homo sapiens festgestellt, dass daraus geschlossen werden konnte, dass der Homo sapiens kein Nachfahre des Neandertalers ist. Dies hat man mit Computerverfahren, ähnlich dem oben vorgestellten, nachweisen können.