Hirschberg-Algorithmus

Der Hirschberg-Algorithmus berechnet das paarweise Sequenzalignment und hat einen zur Eingabe linearen Speicherbedarf. Der in den 1970er Jahren von Dan Hirschberg entwickelte Algorithmus verwendet die Methode der Dynamischen Programmierung und das Divide-and-conquer Prinzip.

Allgemeines[Bearbeiten | Quelltext bearbeiten]

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 der Algorithmus zum Beispiel in der Bioinformatik zum Abgleich verschiedener DNA- oder 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[Bearbeiten | Quelltext bearbeiten]

Zum Verständnis des Hirschberg-Algorithmus ist es zunächst wichtig zu verstehen, dass sich die Levenshtein-Distanz auf 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 klar sein, dass sich diese Berechnung auch rückwärts durchführen lässt. Dabei wird die gedachte Matrix nicht von links nach rechts und von oben nach unten durchlaufen, sondern von rechts unten nach links 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[Bearbeiten | Quelltext bearbeiten]

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[Bearbeiten | Quelltext bearbeiten]

• 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 (englisch, uci.edu [PDF]). 
• Chao, K.M., Hardison, R.C. and Miller, W.: Recent developments in linear-space alignment methods: a survey. In: Journal of Computational Biology. Nr. 4, 1994, S. 271–291 (englisch, edu.tw [PDF]).