DNA-Alignment (Sequenzalignment)

Idee

Da DNA-Sequenzen in der Biologie experimentell gewonnen werden, sind sie häufig mit Meßfehlern behaftet, daher macht es Sinn zwei DNA-Sequenzen nicht "genau" zu vergleichen. Deshalb versucht man bei einem Sequenzalignment zwei DNA-Sequenzen so untereinander zu schreiben, daß möglichst viele Übereinstimmungen gefunden werden. Dabei ist es erlaubt Platzhalter in die DNA-Sequenzen einzufügen, damit die Übereinstimmungen "größer" werden. Als Maß für die Übereinstimmung, werden Punkte für ein "Match", "Mismatch" oder ein "Gap" verteilt. Die Werte dafür sind abhängig von dem verwendeten Algorithmus.

Ein weiterer wichtiger Anwendungsfall für das Sequenzalignment ist die Bestimmung von Verwandtschaftsbeziehungen von zwei Organismen. Will man beispielsweise den Verwandtschaftsgrad zwischen zwei DNA-Sequenzen bestimmen, kann ein Sequenzalignment durchgeführt werden.

Beispiel

Zur Verdeutlichung folgt ein kurzes Beispiel von einem Sequenzalignment.

Sequenz 1: KKCGTTTWCGT

Sequenz 2: KKTTCGTCCTTWCCT

Man beachte dabei, daß in der zweiten Sequenz an der dritten und vierten Stelle ein "T" und an der achten und neunten Stelle ein "C" eingefügt wurden. Außerdem ist an der 14. Stelle, statt eines "G", ein "C" zu sehen. Sonst sind die Sequenzen identisch. Also sollte ein guter Algorithmus für das Sequenzalignment diese beiden "Inserts" und die eine "Mutation" finden.

Ein optimales Alignment für die beiden Sequenzen sieht somit folgendermaßen aus:

KK--CGT--TTWCGT || ||| |||| | KKTTCGTCCTTWCCT

Dabei ist zu beachten, daß das gefundene Alignment auch so aussehen könnte:

KK--CGT--TTWCG-T || ||| |||| | KKTTCGTCCTTWC-CT

Dies ist abhängig davon, ob ein "Mismatch" oder zwei "Gaps" höher bewertet werden.

Blosum62

Wie eben schon angedeutet, gibt es verschiedene Möglichkeiten ein "Match", "Mismatch" oder "Gap" zu bewerten. Diese Bewertung nennt man Scoring-Modelle. Ein sehr einfaches Scoring-Modell ist das folgende:

  • Match +1
  • Mismatch -1
  • Gap 0

Eine weitere Möglichkeit ein "Match" oder "Mismatch" zu bewerten, ist die BLOSUM62-Matrix, bei der neben den 20 Aminosäuren auch die Symbole für Mehrdeutigkeit (B, Z, X) und ein Stop-Codon (*) enthalten sind. Diese Einträge fehlen bei der BLOSUM50-Matrix. Die BLOSON62-Matrix ist in der folgenden Tabelle zu sehen. Die roten Einträge symbolisieren ein "Match". Die restlichen Werte stehen für ein "Mismatch", wobei die lila Einträge verdeutlichen sollen, daß es auch positive Bewertungen für ein "Mismatch" geben kann.

A R N D C Q E G H I L K M F P S T W Y V B Z X *
A 4 -1 -2 -2 0 -1 -1 0 -2 -1 -1 -1 -1 -2 -1 1 0 -3 -2 0 -2 -1 0 -4
R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -1 -3 -2 -1 -1 -3 -2 -3 -1 0 -1 -4
N -2 0 6 1 -3 0 0 0 1 -3 -3 0 -2 -3 -2 1 0 -4 -2 -3 3 0 -1 -4
D -2 -2 1 6 -3 0 2 -1 -1 -3 -4 -1 -3 -3 -1 0 -1 -4 -3 -3 4 1 -1 -4
C 0 -3 -3 -3 9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4
Q -1 1 0 0 -3 5 2 -2 0 -3 -2 1 0 -3 -1 0 -1 -2 -1 -2 0 3 -1 -4
E -1 0 0 2 -4 2 5 -2 0 -3 -3 1 -2 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4
G 0 -2 0 -1 -3 -2 -2 6 -2 -4 -4 -2 -3 -3 -2 0 -2 -2 -3 -3 -1 -2 -1 -4
H -2 0 1 -1 -3 0 0 -2 8 -3 -3 -1 -2 -1 -2 -1 -2 -2 2 -3 0 0 -1 -4
I -1 -3 -3 -3 -1 -3 -3 -4 -3 4 2 -3 1 0 -3 -2 -1 -3 -1 3 -3 -3 -1 -4
L -1 -2 -3 -4 -1 -2 -3 -4 -3 2 4 -2 2 0 -3 -2 -1 -2 -1 1 -4 -3 -1 -4
K -1 2 0 -1 -3 1 1 -2 -1 -3 -2 5 -1 -3 -1 0 -1 -3 -2 -2 0 1 -1 -4
M -1 -1 -2 -3 -1 0 -2 -3 -2 1 2 -1 5 0 -2 -1 -1 -1 -1 1 -3 -1 -1 -4
F -2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 6 -4 -2 -2 1 3 -1 -3 -3 -1 -4
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4 7 -1 -1 -4 -3 -2 -2 -1 -2 -4
S 1 -1 1 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4 1 -3 -2 -2 0 0 0 -4
T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5 -2 -2 0 -1 -1 0 -4
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1 1 -4 -3 -2 11 2 -3 -4 -3 -2 -4
Y -2 -2 -2 -3 -2 -1 -2 -3 2 -1 -1 -2 -1 3 -3 -2 -2 2 7 -1 -3 -2 -1 -4
V 0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4 -3 -2 -1 -4
B -2 -1 3 4 -3 0 1 -1 0 -3 -4 0 -3 -3 -2 0 -1 -4 -3 -3 4 1 -1 -4
Z -1 0 0 1 -3 3 4 -2 0 -3 -3 1 -1 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4
X 0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 0 0 -2 -1 -1 -1 -1 -1 -4
* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 1

Benutzt man das einfache Scoring-Modell, um das Alignment aus dem Beispiel zu bewerten, so erhält man folgende Werte:

Wert 1: +1+1+0+0+1+1+1+0+0+1+1+1+1-1+1 = 9

Wert 2: +1+1+0+0+1+1+1+0+0+1+1+1+1+0+0+1 = 10

Benutzt man dagegen die BLOSUM62-Matrix erhält man folgende Werte:

Wert 1: +5+5+0+0+9+6+5+0+0+5+5+11+9-3+5 = 62

Wert 2: +5+5+0+0+9+6+5+0+0+5+5+11+9+0+0+5 = 65

Dabei wird auch deutlich, daß viel von dem gewählten "Penalty" für ein "Gap" abhängt. Wird der Wert nicht, wie hier, auf 0 gesetzt, sondern beispielsweise auf -5, erhält man ein anderes Ergebnis. D.h. beim ersten Alignment ergibt sich unter Umständen ein höherer Wert, als bei dem zweiten Alignment.

Algorithmen

Für das Sequenzalignment existieren verschiedene Verfahren. Das hier vorgestellte Programm untersützt folgende Algorithmen:

  1. Einfacher Dotplot-Algorithmus
  2. Needleman-Wunsch-Algorithmus (globales Alignment)
  3. Smith-Waterman-Algorithmus (lokales Alignment)

Außerdem existieren noch folgende Verfahren, auf die hier nicht weiter eingegangen wird:

  • Hirschberg-Algorithmus (globales Alignment auf linearem Speicherplatz)
  • FASTA (heuristische Verfahren für paarweises Alignment)
  • BLAST (heuristische Verfahren für paarweises Alignment)

zu 1.)

Beim Dotplot-Algorithmus wird eine sogenannte F-Matrix aufgestellt. Wobei die Matrix folgendermaßen initialisiert wird: F(i,0)=F(0,i)=F(0,0)=0 für alle i={1,..,L}. Die restlichen Werte werden rekursiv berechnet. Dabei kommt folgende Vorschrift zum Einsatz:

Wenn x(i)=y(i): Dann setze F(i,j) = max (F(i-1,j), F(i,j-1), F(i-1,j-1)+1). Wenn x(i)!=y(i): Dann setze F(i,j) = max (F(i-1,j), F(i,j-1), F(i-1,j-1)).

Wurde die F-Matrix mit der Vorschrift aufgestellt, dann wird das Alignment dadurch bestimmt, daß man von rechts unten nach links oben läuft. Dabei wird immer der Eintrag mit dem höchsten Wert genommen.

Als C-Programmcode kann das Aufstellen der F-Matrix folgendermaßen ausgedrückt werden:

int F[lengthA][lengthB]; for(int i=0;i<lengthA;i++) { F[i][0] = 0; } for(int i=0;i<lengthB;i++) { F[0][i] = 0; } for(int i=1;i<lengthA;i++) { for(int j=1;j<lengthB;j++) { if(seqA[i]==seqB[j]) { F[i][j] = max(F[i-1][j], F[i][j-1], F[i-1][j-1]+1); } else { F[i][j] = max(F[i-1][j], F[i][j-1], F[i-1][j-1]); } } }

Das Alignment aus der F-Matrix kann mit dem Code bestimmt werden:

int i = lengthA-1; int j = lengthB-1; while(i!=0 && j!=0) { int dirUP = F[i-1][j]; int dirLEFT = F[i][j-1]; int dirDIAGONAL = F[i-1][j-1]; int maximum = max(dirUP, dirLEFT, dirDIAGONAL); if(dirDIAGONAL==maximum) { if(seqA[i-1]==seqB[j-1]) { cout << seqA[i-1] << "--" << seqB[j-1] << endl; } else { cout << seqA[i-1] << " " << seqB[j-1] << endl; } i--; j--; } else if(dirUP==maximum) { cout << seqA[i-1] << " |" << endl; i--; } else { cout << "| " << seqB[j-1] << endl; j--; } }

Dabei wird das Ergebnis in umgekehrter Richtung und untereinander ausgegeben. Außerdem ist zu beachten, daß die Richtung "Diagonal" die höchste Priorität hat. Gefolgt von der Richtung "Hoch" und "Links".

Die F-Matrix für das Beispiel sieht dann so aus:

K K T T C G T C C T T W C C T 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 K 0 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 C 0 1 1 1 2 3 3 3 3 3 3 3 3 3 3 3 G 0 1 2 2 2 3 4 4 4 4 4 4 4 4 4 4 T 0 1 2 3 3 3 4 4 4 5 5 5 5 5 5 5 T 0 1 2 3 3 3 4 4 4 5 6 6 6 6 6 6 T 0 1 2 3 3 3 4 4 4 5 6 7 7 7 7 7 W 0 1 2 3 4 4 4 5 5 5 6 7 8 8 8 8 C 0 1 2 3 4 5 5 5 5 5 6 7 8 8 8 8 G 0 1 2 3 4 5 6 6 6 6 6 7 8 8 9 9 T 0 1 2 3 4 5 6 6 6 6 6 7 8 8 9 10

Die F-Matrix wird von unten rechts nach oben links durchlaufen, wobei immer der Eintrag mit dem höchsten Wert genommen wird. So erhält man folgendes Ergebnis:

K--KCG--TTTWCGT | || |||| | KKTTCGTCCTTWCCT

Man sieht sofort, daß das gefundene Alignment nicht optimal ist.

zu 2.)

Beim Needleman-Wunsch-Algorithmus wird die F-Matrix folgendermaßen initialisiert: F(i,0)=-i*d, F(0,j)=-j*d, i={1,..,L(x)}, j={1,..,L(y)}. Außerdem gibt es bei diesem Algorithmus eine Z-Matrix in der die Richtungen gespeichert werden, in das das Alignment später zusammengebaut wird. Die Vorschrift für das Aufstellen der F-Matrix und Z-Matrix sieht folgendermaßen aus:

F(i,j) = max( F(i-1,j-1)+s(x(i),y(j)), F(i-1,j)-d, F(i,j-1)-d), wobei d der Penalty ist und die Funktion s auf die BLOSUM62-Matrix zugreift.

Eine mögliche Implementierung zur Bestimmung von F- und Z-Matrix sieht folgendermaßen aus:

for(int i=1;i<lengthA;i++) { for(int j=1;j<lengthB;j++) { int f1 = F[i-1][j-1]+s(seqA[i-1],seqB[j-1]); int f2 = F[i-1][j]-d; int f3 = F[i][j-1]-d; int maximum = max(f1, f2, f3); if(maximum==f1) { Z[i][j] = DIAGONAL; } else if(maximum==f2) { Z[i][j] = UP; } else { Z[i][j] = LEFT; } F[i][j] = maximum; } }

Das Alignment wird mit folgendem Code bestimmt:

int i = lengthA-1; int j = lengthB-1; while(i!=0 || j!=0) { int direction = Z[i][j]; char base1 = seqA[i-1]; char base2 = seqB[j-1]; if(direction==UP) { i--; cout << base1 << " |" << endl; } else if(direction==DIAGONAL) { i--; j--; if(base1==base2) { cout << base1 << "--" << base2 << endl; } else { cout << base1 << " " << base2 << endl; } } else { j--; cout << "| " << base2 << endl; } }

Dabei wird das Ergebnis auch hier in umgekehrter Richtung und untereinander ausgegeben. Außerdem ist zu beachten, daß die Richtung "Hoch" die höchste Priorität hat. Gefolgt von der Richtung "Diagonal" und "Links".

Benutzt man den Needleman-Wunsch-Algorithmus mit einem Penalty von 5, erhält man folgende F-Matrix:

K K T T C G T C C T T W C C T 0 -5 -10 -15 -20 -25 -30 -35 -40 -45 -50 -55 -60 -65 -70 -75 K -5 5 0 -5 -10 -15 -20 -25 -30 -35 -40 -45 -50 -55 -60 -65 K -10 0 10 5 0 -5 -10 -15 -20 -25 -30 -35 -40 -45 -50 -55 C -15 -5 5 9 4 9 4 -1 -6 -11 -16 -21 -26 -31 -36 -41 G -20 -10 0 4 7 4 15 10 5 0 -5 -10 -15 -20 -25 -30 T -25 -15 -5 5 9 6 10 20 15 10 5 0 -5 -10 -15 -20 T -30 -20 -10 0 10 8 5 15 19 14 15 10 5 0 -5 -10 T -35 -25 -15 -5 5 9 6 10 14 18 19 20 15 10 5 0 W -40 -30 -20 -10 0 4 7 5 9 13 16 17 31 26 21 16 C -45 -35 -25 -15 -5 9 4 6 14 18 13 15 26 40 35 30 G -50 -40 -30 -20 -10 4 15 10 9 13 16 11 21 35 37 33 T -55 -45 -35 -25 -15 -1 10 20 15 10 18 21 16 30 34 42

Außerdem ergibt sich die Z-Matrix:

K K T T C G T C C T T W C C T 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 K 2 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 K 2 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 C 2 2 2 1 1 1 3 3 1 1 3 3 3 1 1 3 G 2 2 2 2 1 2 1 3 3 3 3 3 3 3 3 3 T 2 2 2 1 1 1 2 1 3 3 1 1 3 3 3 1 T 2 2 2 1 1 1 2 1 1 1 1 1 3 3 3 1 T 2 2 2 1 1 1 1 1 1 1 1 1 3 3 3 1 W 2 2 2 2 2 2 1 2 2 2 1 1 1 3 3 3 C 2 2 2 2 2 1 3 1 1 1 3 1 2 1 1 3 G 2 2 2 2 2 2 1 3 2 2 1 1 2 2 1 1 T 2 2 2 1 1 2 2 1 3 3 1 1 2 2 1 1

Wobei die Zahlen folgende Bedeutung haben:

  • 1: Diagonal gehen
  • 2: Nach oben gehen
  • 3: Nach links gehen

Mit der F- und der Z-Matrix erhält man das Alignment:

KK--CGT--TTWCGT || ||| |||| | KKTTCGTCCTTWCCT

Man sieht sofort, daß die "Inserts" und die einzelne "Mutation" gefunden wurden.

zu 3.

Benutzt man den Smith-Waterman-Algorithmus muß dabei auch eine F- und Z-Matrix bestimmt werden. Initialisiert wird die F-Matrix nach folgender Regel:

F(i,0)=F(0,j)=F(0,0)=0, i={1,..,L(x)}, j={1,..,L(y)}

Die Vorschrift für das Erstellen der F- und Z-Matrix sieht folgendermaßen aus:

F(i,j)= max( 0, F(i-1,j-1)+s(x(i),y(j)), F(i-1,j)-d, F(i,j-1)-d), wobei d das Penalty und s die Funktion für den Zugriff auf die BLOSUM62-Matrix ist.

Zu beachten ist, daß durch die 0 in der Maximum-Funktion dafür gesorgt wird, daß kein Wert in der F-Matrix kleiner als Null ist. Gleichzeitig bedeutet dies, daß in der Z-Matrix keine Richtung angegeben wird.

Der Programmcode für die F- und Z-Matrix kann so aussehen:

for(int i=1;i<lengthA;i++) { for(int j=1;j<lengthB;j++) { int f1 = F[i-1][j-1]+s(seqA[i-1],seqB[j-1]); int f2 = F[i-1][j]-d; int f3 = F[i][j-1]-d; int maximum = max(f1, f2, f3); if(maximum<0) { maximum = 0; } if(maximum==f2) { Z[i][j] = UP; } else if(maximum==f3){ Z[i][j] = LEFT; } else if(maximum==f1) { Z[i][j] = DIAGONAL; } else { Z[i][j] = NOTHING; } F[i][j] = maximum; } }

Der Code für die Ausgabe des Alignment kann folgendermaßen aussehen:

while(i!=0 || j!=0) { int direction = Z[i][j]; char base1 = seqA[i-1]; char base2 = seqB[j-1]; if(direction==UP) { i--; cout << base1 << " |" << endl; } else if(direction==DIAGONAL) { i--; j--; if(base1==base2) { cout << base1 << "--" << base2 << endl; } else { cout << base1 << " " << base2 << endl; } } else { j--; cout << "| " << base2 << endl; } }

Dabei wird das Ergebnis auch hier in umgekehrter Richtung und untereinander ausgegeben. Außerdem ist zu beachten, daß die Richtung "Hoch" die höchste Priorität hat. Gefolgt von der Richtung "Diagonal" und "Links".

Benutzt man den Smith-Waterman-Algorithmus mit einem Penalty von 5, erhält man folgende F-Matrix:

K K T T C G T C C T T W C C T 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K 0 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 K 0 5 10 5 0 0 0 0 0 0 0 0 0 0 0 0 C 0 0 5 9 4 9 4 0 9 9 4 0 0 9 9 4 G 0 0 0 4 7 4 15 10 5 6 7 2 0 4 6 7 T 0 0 0 5 9 6 10 20 15 10 11 12 7 2 3 11 T 0 0 0 5 10 8 5 15 19 14 15 16 11 6 1 8 T 0 0 0 5 10 9 6 10 14 18 19 20 15 10 5 6 W 0 0 0 0 5 8 7 5 9 13 16 17 31 26 21 16 C 0 0 0 0 0 14 9 6 14 18 13 15 26 40 35 30 G 0 0 0 0 0 9 20 15 10 13 16 11 21 35 37 33 T 0 0 0 5 5 4 15 25 20 15 18 21 16 30 34 42

Und die Z-Matrix:

K K T T C G T C C T T W C C T 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 K 2 1 1 3 0 0 0 0 0 0 0 0 0 0 0 0 K 2 1 1 3 3 0 0 0 0 0 0 0 0 0 0 0 C 2 2 2 1 3 1 3 0 1 1 3 0 0 1 1 3 G 2 0 2 2 1 2 1 3 3 1 1 3 0 2 1 1 T 2 0 0 1 1 1 2 1 3 3 1 1 3 3 1 1 T 2 0 0 1 1 1 2 2 1 3 1 1 3 3 3 1 T 2 0 0 1 1 1 1 2 2 1 1 1 3 3 3 1 W 2 0 0 2 2 1 1 2 2 2 1 1 1 3 3 3 C 2 0 0 0 2 1 3 1 1 1 3 1 2 1 3 3 G 2 0 0 0 0 2 1 3 3 2 1 3 2 2 1 1 T 2 0 0 1 1 2 2 1 3 3 1 1 2 2 1 1

Als Alignment ergibt sich auch:

KK--CGT--TTWCGT || ||| |||| | KKTTCGTCCTTWCCT

Also hat auch dieser Algorithmus das optimale Alignment gefunden.

Jetzt kann jeder selbst mal mit verschiedenen Sequenzen und Penalties herumspielen und überprüfen, welcher Algorithmus für welche Sequenzen bessere Ergebnisse liefert.

Das Programm

Das Programm ist ein kommandozeilenbasiertes Tool, das mit verschiedenen Parametern gesteuert werden kann.

Der Aufruf sieht folgendermaßen aus:

./DNAAlignment -switch outputFile sequenceFile1 sequenceFile2 [penalty]

Wobei als switch folgende Werte gesetzt werden können:

  • -simple: Einfaches Alignment
  • -nw: Needleman-Wunsch-Algorithmus
  • -sw: Smith-Waterman-Algorithmus

Wenn der switch nw oder sw gesetzt wird, muß auch ein Wert für den "Penalty" angegeben werden.

Mit dem Paramenter outputFile gibt man eine Datei an, in die das Ergebnis des Alignments geschrieben werden soll. Achtung: Es wird nicht überprüft, ob die Datei vorhanden ist. Sie wird einfach überschrieben!

Die Parameter sequenceFile1 und sequenceFile2 geben die beiden Dateien an, in denen sich die Eingabesequenzen befinden. Die Eingabesequenzen müssen sich in der ersten Zeile der Datei befinden, d.h. die Sequenz darf nicht durch Zeilenumbrüche oder ähnliches getrennt werden!

Das Programm sollte ausschließlich für Übungszwecke eingesetzt werden, da es nicht dafür erstellt wurde, um große Sequenzen zu verarbeiten. Dies liegt daran, daß die beschriebenen Algorithmen recht speicherintensiv sind. Will man ein reales Sequenzalignment durchführen, dann sollte BLAST oder FASTA verwendet werden.

Download

Das Programm wurde mit dem GCC unter Cygwin übersetzt und benötigt deshalb die Cygwin-Dll, für die Ausführung. Diese Dll ist Bestandteil des ZIP-Pakets. Download: DNA-Alignment