Hirschberg-Algorithmus

Der Hirschberg-Algorithmus berechnet d​as paarweise Sequenzalignment u​nd hat e​inen zur Eingabe linearen Speicherbedarf. Der i​n den 1970er Jahren v​on Dan Hirschberg entwickelte Algorithmus verwendet d​ie Methode d​er Dynamischen Programmierung u​nd das Divide-and-conquer Prinzip.

Allgemeines

Der Hirschberg-Algorithmus ist ein allgemein einsetzbarer und optimaler Algorithmus zum Auffinden eines Sequenzalignment. Der bekannte BLAST-Algorithmus und der FASTA-Algorithmus sind nur suboptimale Heuristiken. Vergleicht man den Hirschberg-Algorithmus mit dem Needleman-Wunsch-Algorithmus, so handelt es sich beim Hirschberg-Algorithmus weniger um einen komplett neuen Algorithmus, sondern eher um eine clevere Strategie, die den Needleman-Wunsch-Algorithmus geschickt einsetzt, um den Speicherverbrauch zu linearisieren, was auch das Besondere an diesem Algorithmus ist: Die Berechnungen für ein Sequenzalignment benötigen nur linear viel Speicherplatz, womit die Platzkomplexität des Algorithmus in O(n) liegt. Zur Berechnung eines Alignments zweier Zeichenketten ${\displaystyle x}$ und ${\displaystyle y}$ mit ${\displaystyle m=|x|}$ und ${\displaystyle n=|y|}$ besitzt der Algorithmus eine Laufzeit von ${\displaystyle \Theta (mn)}$ und einen Speicherverbrauch von ${\displaystyle \Theta (\min\{m,n\})}$. O.B.d.A soll im Folgenden ${\displaystyle n\leq m}$ gelten, so dass der Platzbedarf in ${\displaystyle \Theta (n)}$ liegt.

Anwendung findet d​er Algorithmus z​um Beispiel i​n der Bioinformatik z​um Abgleich verschiedener DNA- o​der Proteinsequenzen.

In einer leicht abgewandelten Form wird Hirschbergs Algorithmus auch dazu verwendet um in einem Graphen parallel Zusammenhangskomponenten mit Aufwand ${\displaystyle \Theta (\log ^{2}n)}$ auf ${\displaystyle \Theta (n^{2})}$ Prozessoren zu berechnen.

Berechnung der Levenshtein-Distanz auf linearem Speicherplatz

Zum Verständnis d​es Hirschberg-Algorithmus i​st es zunächst wichtig z​u verstehen, d​ass sich d​ie Levenshtein-Distanz a​uf linearem Speicherplatz berechnen lässt:

01 ${\displaystyle T_{0}}$ := 0
02 for j in 1..n loop
03       ${\displaystyle T_{j}}$ := ${\displaystyle T_{j-1}}$ + ${\displaystyle Ins(y_{j})}$
04 end loop
05 for i in 1..m loop
06       s := ${\displaystyle T_{0}}$
07       ${\displaystyle T_{0}}$ := ${\displaystyle T_{0}}$ + ${\displaystyle Del(x_{i})}$
08       c := ${\displaystyle T_{0}}$
09       for j in 1..n loop
10             c := ${\displaystyle \min {\begin{cases}s&+Sub(x_{i},y_{j})\\T_{j}&+Del(x_{i})\\c&+Ins(y_{j})\end{cases}}}$
11             s := ${\displaystyle T_{j}}$
12             ${\displaystyle T_{j}}$ := c
13       end loop
14 end loop

In den Zeilen 1–4 wird das eindimensionale Feld ${\displaystyle T}$ mit n Plätzen Speicherbedarf initialisiert. In Zeile 6 wird die Initialisierung des ersten Elements ${\displaystyle T_{0}}$ in ${\displaystyle s}$ gerettet. Danach wird ${\displaystyle T_{0}}$ und ${\displaystyle c}$ mit dem Startwert für die nächste Zeile initialisiert. Die nachfolgende Abbildung zeigt eine Momentaufnahme eines Zeilendurchlaufs. In der inneren Schleife zeigt ${\displaystyle c}$ immer auf das jeweils zuvor berechnete Ergebnis, während ${\displaystyle s}$ das noch benötigte Ergebnis der letzten Zeile sichert. Nach Zeile 14 steht die Levenshtein-Distanz als Ergebnis in ${\displaystyle T_{n}}$.

  ε ${\displaystyle y_{1}}$ ${\displaystyle y_{2}}$ ${\displaystyle y_{3}}$ ${\displaystyle y_{4}}$ ...
ε    0  1  2  3  ...
${\displaystyle x_{1}}$   1
${\displaystyle x_{2}}$
...

s = 0
c = ${\displaystyle T_{0}}$ = 1

Es sollte k​lar sein, d​ass sich d​iese Berechnung a​uch rückwärts durchführen lässt. Dabei w​ird die gedachte Matrix n​icht von l​inks nach rechts u​nd von o​ben nach u​nten durchlaufen, sondern v​on rechts u​nten nach l​inks oben:

01 ${\displaystyle T_{n}}$ := 0
02 for j in n-1..0 loop
03       ${\displaystyle T_{j}}$ := ${\displaystyle T_{j+1}}$ + ${\displaystyle Ins(y_{j+1})}$
04 end loop
05 for i in m-1..0 loop
06       s := ${\displaystyle T_{n}}$
07       ${\displaystyle T_{n}}$ := ${\displaystyle T_{n}}$ + ${\displaystyle Del(x_{i+1})}$
08       c := ${\displaystyle T_{n}}$
09       for j in n-1..0 loop
10             c := ${\displaystyle \min {\begin{cases}s&+Sub(x_{i+1},y_{j+1})\\T_{j}&+Del(x_{i+1})\\c&+Ins(y_{j+1})\end{cases}}}$
11             s := ${\displaystyle T_{j}}$
12             ${\displaystyle T_{j}}$ := c
13       end loop
14 end loop

Berechnung des Alignments auf linearem Speicherplatz

Der Divide-&-Conquer-Algorithmus von Hirschberg berechnet ein Alignment der Zeichenketten ${\displaystyle |x|}$ und ${\displaystyle |y|}$, indem er Vorwärts- und Rückwärtsdurchlauf miteinander kombiniert (Zeilenangaben beziehen sich auf den nachfolgend angegebenen Pseudocode):

1. Wenn ${\displaystyle |x|=1}$ oder ${\displaystyle |y|=1}$ liegt ein triviales Alignment-Problem vor (Zeilen 14 – 22). Ein String bestehend aus nur einem Zeichen muss auf einen anderen String ausgerichtet werden und ein Alignment wird zurückgegeben. Ist ${\displaystyle |x|>1}$ und ${\displaystyle |y|>1}$ geht man über zu Schritt 2.

2. Ein Vorwärtsdurchlauf berechnet ein Alignment von ${\displaystyle y}$ und der ersten Hälfte von ${\displaystyle x}$ (Zeilen 27 – 40). Das Ergebnis des Vorwärtsdurchlaufs ist ein Feld ${\displaystyle T^{\ell }}$, dessen Elemente die Kosten für einen Durchlauf von ${\displaystyle (0,0)}$ bis ${\displaystyle (|x|/2,j)}$ (mit ${\displaystyle 0\leq j\leq n}$) angeben.

3. Ein Rückwärtsdurchlauf berechnet ein Alignment von ${\displaystyle y}$ mit der zweiten Hälfte von ${\displaystyle x}$ (Zeilen 42 – 55). Das Ergebnis ist ein weiteres Feld ${\displaystyle T^{r}}$, dessen Elemente die Kosten für einen Durchlauf von ${\displaystyle (n,m)}$ zurück zu ${\displaystyle (|x|/2,j)}$ angeben.

4. In den Feldelementen ${\displaystyle T^{\ell }(n)}$ und ${\displaystyle T^{r}(0)}$ stehen die beiden Levenshtein-Distanzen für die globalen Alignments von ${\displaystyle y}$ und den jeweiligen Hälften von ${\displaystyle x}$. In den restlichen Elementen von ${\displaystyle T^{\ell }}$ stehen die Distanzen von der ersten ${\displaystyle x}$-Hälfte zu allen Präfixen von ${\displaystyle y}$. Entsprechend enthalten die restlichen Elemente von ${\displaystyle T^{r}}$ die Distanzen von der zweiten ${\displaystyle x}$-Hälfte zu allen Suffixen von ${\displaystyle y}$.

5. Die Levenshtein-Distanz von ${\displaystyle x}$ zu ${\displaystyle y}$ kann nun errechnet werden, indem man entlang der mittleren Zeile der Alignment-Matrix läuft und nach einer kleinsten Summe von korrespondierenden ${\displaystyle T^{\ell }}$- und ${\displaystyle T^{r}}$-Elementen sucht. Ist eine solche minimale Summe gefunden, hat man eine Position in der mittleren Zeile gefunden, in der das optimale Alignment die mittlere Zeile bzw. die Mitte von ${\displaystyle x}$ schneidet. An dieser Stelle wird ${\displaystyle y}$ in zwei Teile zerteilt und damit kann auch das Alignment-Problem in zwei kleinere Alignment-Probleme zerteilt werden (Zeilen 59 – 65).

6. Schritt 1 wird rekursiv auf den beiden Teilen von ${\displaystyle x}$ und ${\displaystyle y}$ aufgerufen. Die beiden rekursiven Aufrufe geben Teil-Alignments zurück, die zu einem einzigen Alignment verknüpft werden. Das Alignment wird zurückgegeben (Zeilen 68 und 69).

01 --
02 -- Der Divide-&-Conquer-Algorithmus von Hirschberg zur
03 -- Berechnung des globalen Alignments auf linearem Speicher.
04 --
05 -- Bei ${\displaystyle m=|x|,n=|y|,n\leq m}$ besitzt der Algorithmus eine Laufzeit von ${\displaystyle \Theta (nm)}$
06 -- und einen Speicherverbrauch von ${\displaystyle \Theta (n)}$.
07 --
08 function HirschbergAlignment(x,y : string) return A is
09        function SubAlignment(${\displaystyle i_{1}}$,${\displaystyle j_{1}}$,${\displaystyle i_{2}}$,${\displaystyle j_{2}}$ : integer) return A is
10                mitte,cut : integer
11                s,c : real
12                ${\displaystyle T^{\ell },T^{r}}$ : array(${\displaystyle j_{1}}$..${\displaystyle j_{2}}$) of real
13        begin
14                if ${\displaystyle i_{1}+1=i_{2}}$ or ${\displaystyle j_{1}=j_{2}}$ then
15                        -- Konstruiere Matrix T für die Teil-Strings
16                        -- x(${\displaystyle i_{1}+1}$..${\displaystyle i_{2}}$) und y(${\displaystyle j_{1}+1}$..${\displaystyle j_{2}}$)
17                        -- Achtung: Nur linearer Speicherplatz erforderlich!
18                        T := ...
19                        -- Berechne triviales Alignment auf Matrix T
20                        -- in linearer Laufzeit
21                        return Alignment(T,x(${\displaystyle i_{1}+1}$..${\displaystyle i_{2}}$),y(${\displaystyle j_{1}+1}$..${\displaystyle j_{2}}$))
22                end if
23
24                mitte := ${\displaystyle (i_{1}+i_{2})/2}$
25                -- finde ausgehend von ${\displaystyle (i_{1},j_{1})}$ den minimalen Pfad
26                -- mit dem Vorwärtsalgorithmus:
27                ${\displaystyle T^{\ell }(j_{1})}$ := 0
28                for j in ${\displaystyle j_{1}+1}$..${\displaystyle j_{2}}$ loop
29                        ${\displaystyle T^{\ell }(j)}$ := ${\displaystyle T^{\ell }(j-1)+Ins(y_{j})}$
30                end loop
31                for i in ${\displaystyle i_{1}+1}$..mitte loop
32                        s := ${\displaystyle T^{\ell }(j_{1})}$
33                        c := ${\displaystyle T^{\ell }(j_{1})+Del(x_{i})}$
34                        ${\displaystyle T^{\ell }(j_{1})}$ := c
35                        for j in ${\displaystyle j_{1}+1}$..${\displaystyle j_{2}}$ loop
36                                c := ${\displaystyle \min {\begin{cases}T^{\ell }(j)&+Del(x_{i})\\s&+Sub(x_{i},y_{j})\\c&+Ins(y_{j})\end{cases}}}$
37                                s := ${\displaystyle T^{\ell }(j)}$
38                                ${\displaystyle T^{\ell }(j)}$ := c
39                        end loop
40                end loop
41                -- finde minimalen score-pfad nach ${\displaystyle (i_{2},j_{2})}$
42                ${\displaystyle T^{r}(j_{2})}$ := 0
43                for j in ${\displaystyle j_{2}-1}$..${\displaystyle j_{1}}$ loop
44                        ${\displaystyle T^{r}(j)}$ := ${\displaystyle T^{r}(j+1)+Ins(y_{j+1})}$
45                end loop
46                for i in ${\displaystyle i_{2}-1}$..mitte loop
47                        s := ${\displaystyle T^{r}(j_{2})}$
48                        c := ${\displaystyle T^{r}(j_{2})+Del(x_{i+1})}$
49                        ${\displaystyle T^{r}(j_{2})}$ := c;
50                        for j in ${\displaystyle j_{2}-1}$..${\displaystyle j_{1}}$ loop
51                                c := ${\displaystyle \min {\begin{cases}T^{r}(j)&+Del(x_{i+1})\\s&+Sub(x_{i+1},y_{j+1})\\c&+Ins(y_{j+1})\end{cases}}}$
52                                s := ${\displaystyle T^{r}(j)}$
53                                ${\displaystyle T^{r}(j)}$ := c
54                        end loop
55                end loop
56                -- finde den Punkt aus ${\displaystyle j_{1}}$..${\displaystyle j_{2}}$ in dem der Minimale Pfad die
57                -- mittlere Zeile schneidet:
58                -- ${\displaystyle cut:=_{def}{\mbox{argmin}}_{j_{1}\leq j\leq j_{2}}(T^{\ell }(j)+T^{r}(j))}$
59                for j in ${\displaystyle j_{1}}$..${\displaystyle j_{2}}$ loop
60                        if j=${\displaystyle j_{1}}$ then
61                                cut := ${\displaystyle j_{1}}$
62                        elsif ${\displaystyle T^{\ell }(j)+T^{r}(j)<T^{\ell }(cut)+T^{r}(cut)}$ then
63                                cut := j
64                        end if
65                end loop
66                -- Alignment entsteht durch Konkatenation von linkem und
67                -- rechtem Teil-Alignment:
68                return SubAlignment(${\displaystyle i_{1}}$,${\displaystyle j_{1}}$,mitte,cut)
69                                ${\displaystyle \star }$ SubAlignment(mitte,cut,${\displaystyle i_{2}}$,${\displaystyle j_{2}}$)
70        end SubAlignment
71        m,n : integer
72 begin
73        m := ${\displaystyle |x|}$; n := ${\displaystyle |y|}$
74        -- Sonderbehandlung: ${\displaystyle x}$ ist der leere String und lässt keine Zerteilung zu:
75        if m=0 then
76                return ${\displaystyle {\begin{pmatrix}-\\y_{1}\end{pmatrix}}\star {\begin{pmatrix}-\\y_{2}\end{pmatrix}}\star \cdots \star {\begin{pmatrix}-\\y_{n}\end{pmatrix}}}$
77        else
78                return SubAlignment(0,0,m,n)
79        end if
80 end HirschbergAlignment

Literatur

• D. S. Hirschberg: A linear space algorithm for computing maximal common subsequences. In: Communications of the ACM. Band 18, Nr. 6, 1975, S. 341–343 (PDF).
• Chao, K.M., Hardison, R.C. and Miller, W.: Recent developments in linear-space alignment methods: a survey. In: Journal of Computional Biology. Nr. 4, 1994, S. 271–291 (PDF).