Globales Sequenzalignment mit dem Needleman-Wunsch-Algorithmus in C#

Wenn es darum geht in der Bioinformatik eine funktionelle oder evolutionäre Verwandschaft zwischen Nukleotid- oder Aminosäuresequenzen nachzuweisen, so ist der Needleman-Wunsch-Algorithmus wohl einer der ersten über die man stolpert. Nukleotide und Aminosäuren werden durch Buchstaben dargestellt und müssen so angeordnet werden, dass man Mutationen (dargestellt durch Fehlpaarung), Deletionen oder Insertionen (dargestellt durch Gaps, „-„) ablesen kann. Dabei gilt, mehr gleiche oder ähnliche Elemente in gleicher Reihenfolge, bedeuten mehr Übereinstimmung.

Der Needleman-Wunsch-Algorithmus bedient sich der Methode der Dynamischen Programmierung. Hierbei wird eine Scoring-Matrix initialisiert anhand welcher dann das Alignement selbst erstellt wird. Das ganze Verfahren kann man generell in drei Schritte unterteilen.

  1. Initialisierung der Scoring-Matrix (befülle Matrix mit initialen Werten)
  2. Berechnung der Scoring-Matrix (berechne alle Werte in der Matrix)
  3. Traceback (benutze die Matrix um ein Alignement zu erstellen)

Doch was brauchen wir eigentlich alles für so ein Alignment? Nun zuerst brauchen wir natürlich die zwei Strings die wir vergleichen wollen.

a=„ACTCCTTAA“

b=„ATCCAA“

Nun brauchen wir noch drei weitere Werte. Die Gapkosten (hier f) für die Bewertung einer Insertion oder Deletion und den Match- bzw. den Mismatchscore. Wir nehmen:

f = -1

MatchScore =  1

MismatchScore = -1

Außerdem brauchen wir noch eine Bewertungsfunktion. Diese bekommt zwei Elemente  der Sequenz als Parameter. Falls diese identisch sind liefert die Funktion den Match-Score, wenn nicht dann den Mismatch-Score. Die Bewertungsfunktion ist wie folgt definiert:

Hier steht m für den MatchScore und mm für den MismatchScore.

Hier steht m für den MatchScore und mm für den MismatchScore.

Nun können wir mit der Initialisierung der beginnen. Nehmen wir zunächst eine  Matrix F(m+1 x n+1), wobei m die Länge des Strings a und n die Länge des Strings b ist. Dabei ist F(0,0) = 0.

Schritt 1 der Initialisierung.

Schritt 1 der Initialisierung.

Im zweiten und letzten Schritt der Initialisierung wird wie folgt vorgegangen:

F(i,0) = F(i -1, 0) – 1 und F(0,j) = F(0, j – 1) – 1

Das bedeutet jede Zelle der Ersten Zeile enthält den Wert der vorhergehenden Zelle – 1. Jede Zelle der ersten Spalte enthält den Wert der vorhergehenden Zelle – 1.

matrix2

Nun können wir uns dem zweiten Schritt des Algorithmus widmen, dem Befüllen der Matrix. Dabei wird der Wert in der Zelle F(m+1 x n+1) (also ganz unten rechts) als der Score der Matrix bezeichnet. Die einzelnen Werte werden wie folgt berechnet.

F(i,j) ist die aktuelle Zelle. w(x,y) ist die Bewertungsfunktion. f steht für die Gapkosten.

F(i,j) ist die aktuelle Zelle, w(x,y) ist die Bewertungsfunktion, f steht für die Gapkosten.

Das sieht auf den ersten Blick nicht trivial aus, ist aber recht simpel. Nehmen wir dazu als Beispiel gleich die erste Zelle F(1,1) hier werden die Buchstaben „A“ und „A“ verglichen. Jetzt werden drei Werte berechnet, jeweils aus den Vorgängerzelle. Als erstes gucken wir uns die diagonale Vorgängerzelle an, F(0,0). Hierzu benutzen wir die erste Zeile in der Formel.

(Diagonal)F(i-1,j-1) + w(x,j) = F(0,0) + w(‚A‘,’A‘) = 0 + 1 = 1

Wie oben bereits erwähnt bekommt w(x,y) zwei Buchstaben übergeben und liefert uns entweder den Match- oder Mismatchscore zurück. In diesem Fall klar den Matchscore weil ‚A’=’A‘. Nun nehmen wir uns die obere und linke Vorgängerzelle vor.

(Links) F(i-1,j)+f = F(0,1) + (-1) = -1 + (-1) = -2

(Oben) F(i,j-1)+f = F(1,0)+(-1) = -1 + (-1) = -2

Wie in der Formel beschrieben gilt es jetzt den größten dieser Werte für die aktuelle Zelle zu nehmen also ist F(i,j) = F(1,1) = 1.

Der Wert der Zelle F(1,1) wird aus den drei Vorgängerzellen errechnet.

Der Wert der Zelle F(1,1) wird aus den drei Vorgängerzellen errechnet.

An dieser Stelle möchte ich anmerken, dass es mehrere Möglichkeiten gibt das Traceback auszuführen. Ich mache hier etwas um mir das Traceback zu erleichter. Bei der Berechnung der Matrix, fülle ich gleichzeitig eine zweite Matrix aus. Nur in dieser Matrix speichere ich einfach Buchstaben, quasi eine Matrix die sich merkt welche der Vorgängerzellen den aktuellen Wert geliefert haben. Dabei steht ‚D‘ für diagonal, ‚L‘ für Links und ‚U‘ für Up, also oben :). Da die diagonale Zelle den aktuellen Wert 1 geliefert hat tragen wir in unsere zweite Matrix F(1,1) = ‚D‘ ein. So sollte unsere zweite Matrix zu diesem Zeitpunkt aussehen.

zweitematrix

Die erste Zeile und Spalte werden bereits bei der Initialisierung ausgefüllt und sehen immer so aus.

Wichtig dabei ist auch, dass ich in die allererste Zelle den Wert 0 eintrage. Warum wir das Ganze machen und wie es uns später hilft sehen wir gleich wenn wir zum Traceback dem letzten Schritt des Algorithmus kommen. Hier könnt ihr aber zuerst sehen wie die beiden Matrizen aussehen wenn wir sie komplett ausgefüllt haben.

Der Score der Matrix beträgt 3.

Der Score der Matrix beträgt 3.

Die dazugehörige Traceback-Matrix.

Die dazugehörige Traceback-Matrix.

An dieser Stelle können wir mit dem Traceback beginnen. Dieses wird uns das Alignement liefern. Das funktioniert indem wir beim Score, also rechts unten anfangen und den „am besten bewerteten“ Weg zurück verfolgen. Wie bereits erwähnt gibt es dazu mehrere Ansätze. Wir könnten das Alignement aus der Scoring-Matrix auslesen indem wir berechnen welche der Vorgängerzellen jeweils den aktuellen Wert geliefert hat. Ich habe mich hier entschieden aber eine zweite Matrix, die Traceback Matrix zu benutzen. Hier habe wir in jeder Zelle bereits vermerkt welche Zellen die Werte geliefert haben. Jetzt gehen wir Zelle für Zelle durch bis wir die Zelle F(0,0) erreichen. Dabei gehen wir wie folgt vor:

D → Die diagonale Vorgängerzelle F(i-1,j-1) hat den Wert geliefert → Die beiden Elemente (Buchstaben) werden jeweils Alignement a und b hinzugefügt. Wir bewegen uns auf die Zelle F(i-1,j-1).

U → Die obere Vorgängerzelle F(i,j-1) hat den Wert geliefert → Alignement a wird eine Gap (‚-‚) hinzugefügt. b bekommt den jeweiligen Buchstaben dazu. Wir bewegen uns auf die Zelle F(i,j-1).

L → Die linke Vorgängerzelle F(i-1,j) hat den Wert geliefert → Alignement a wird der entsprechende Buchstabe hinzugefügt. b bekommt eine Gap (‚-‚) dazu. Wir bewegen uns auf die Zelle F(i-1,j).

Fangen wir mit dem Score rechts unten an. Die Zelle hat den Wert ‚D‘. Wir bewegen uns auf die diagonale Zelle und fügen beiden Alignements die Buchstaben hinzu.

‚A‘

‚A‘

Die nächste Zelle hat wieder den Wert ‚D‘. Wir bewegen uns wieder auf die Diagonale Vorgengängerzelle und fügen den Strings die dazugehörigen Buchstaben an.

‚AA‘

‚AA‘

Diese Zelle hat den Wert ‚L‘. Wir bewegen uns auf die linke Zelle. Alignement a bekommt den Buchstaben. b wird eine Gap angehängt.

‚TAA‘

‚-AA‘

Wichtig ist, dass wir die jeweiligen Zeichen immer an den Anfang der Alignements setzten oder aber an das Ende dann aber das Alignement umdrehen. So gehen wir vor bis wir die Zelle F(0,0) erreicht haben und unser Alignement fertig ist.

traceback2

Das Traceback fängt beim Score an und endet in der Zelle (0,0).

Der Algorithmus liefert dieses Alignment.

Der Algorithmus liefert dieses Alignement.

Wichtig hier ist, dass dies nicht die einzige Lösung ist. Abhängig davon welche Werte bei der Erstellung der Matrix zuerst geprüft werden, kann das Traceback variieren. Beispielsweise wenn zwei Zellen den gleichen Wert liefern.

Bei der Implementierung habe ich beschlossen den ganzen Algorithmus in drei Methoden aufzuteilen. Zuerst die einfachste, die Bewertungsfunktion.

//Bekommt die zu vergleichenden Buchstaben, den Match- und Mismatch-Score als Parameter
        public static int ScoreFunction(char a,char b,int matchScore,int mismatchScore)
        {
            //Wenn die Buchstaben gleich sind ist der Match-Score der Rückgabewert.
            //Wenn nicht gib den Mismatch-Score zurück.
            return a == b ? matchScore : mismatchScore;
        }

In der wichtigsten Methode Align(), findet der Hauptteil des Algorithmus statt. Hier findet man die Initialisierung, Bewertung und Traceback…

//Die Align-Methode bekommt die beiden Sequenzen A und B, den Gap-Penalty, Match- und Mismatch-Score als Parameter.
        public static string[] Align(string sequenceA, string sequenceB, int gapPenalty, int matchScore, int mismatchScore)
        {
            //Hier werden beide Matrizen initialisiert. Die Scoring- und Traceback-Matrix.
            #region Initialize
            int[,] matrix = new int[sequenceA.Length + 1, sequenceB.Length + 1];
            char[,] tracebackMatrix = new char[sequenceA.Length + 1, sequenceB.Length + 1];
            matrix[0, 0] = 0;

            //Befülle die erste Reihe der Matrizen.
            for (int i = 1; i < sequenceA.Length + 1;i++)
            {
                matrix[i,0] = matrix[i-1,0] + gapPenalty;
                tracebackMatrix[i,0] = 'L';
            }

            //Befülle die erste Spalte der Matrizen.
            for (int i = 1; i < sequenceB.Length + 1; i++)
            {
                matrix[0, i] = matrix[0 , i - 1] + gapPenalty;
                tracebackMatrix[0,i] = 'U';
            }
            #endregion

            //Hier werden die Werte für die Scoring-Matrix berechnet. Zeitgleich befüllen
            //wir die Traceback-Matrix.
            #region Scoring
            for (int i = 1; i < sequenceA.Length + 1;i++)
            {
                for (int j = 1; j < sequenceB.Length + 1;j++)
                {
                    //Lese und berechne die Werte der Vorgängerzellen.
                    int diagonal = matrix[i - 1, j - 1] + ScoreFunction(sequenceA[i-1],sequenceB[j-1],matchScore,mismatchScore);
                    int links = matrix[i - 1, j] + gapPenalty;
                    int oben = matrix[i, j - 1] + gapPenalty;

                    //Der höchste Wert wird eingetragen.
                    matrix[i, j] = Math.Max(oben,Math.Max(links, diagonal));

                    //Finde raus welche der Zellen den Wert geliefert hat und trage
                    //dies in die Traceback-Matrix ein.
                    if(matrix[i,j] == diagonal && i > 0 && j > 0)
                    {
                        tracebackMatrix[i, j] = 'D';
                    }
                    else if (matrix[i, j] == links)
                    {
                        tracebackMatrix[i, j] = 'L';
                    }
                    else if (matrix[i, j] == oben)
                    {
                        tracebackMatrix[i, j] = 'U';
                    }
                }
            }
            #endregion

            //Hier wird einfach die Traceback-Methode aufgerufen und liefert uns das Alignment.
            #region Traceback
            return TraceBack(tracebackMatrix,sequenceA,sequenceB);
            #endregion

        }

Die Methode Align() macht Gebrauch von der Traceback-Methode. Hier wird der Weg vom Score zur Zelle (0,0) ausgelesen und entsprechend das Alignement erstellt:

//Eine TraceBack-Matrix und die beiden Sequenzen sind die Paramter für das Traceback.
        public static string[] TraceBack(char[,] tracebackMatrix,string sequenzA, string sequenzB)
        {

            int i = tracebackMatrix.GetLength(0) - 1;
            int j = tracebackMatrix.GetLength(1) - 1;

            StringBuilder alignedSeqA = new StringBuilder();
            StringBuilder alignedSeqB = new StringBuilder();

            //Das Traceback wird asugeführt solange wir uns nicht in der Zelle (0,0) befinden.
            while (tracebackMatrix[i, j] != 0)
            {
                switch (tracebackMatrix[i, j])
                {
                    case 'D':
                        alignedSeqA.Append(sequenzA[i - 1]);
                        alignedSeqB.Append(sequenzB[j - 1]);
                        i--;
                        j--;
                        break;
                    case 'U':
                        alignedSeqA.Append("-");
                        alignedSeqB.Append(sequenzB[j - 1]);
                        j--;
                        break;
                    case 'L':
                        alignedSeqA.Append(sequenzA[i - 1]);
                        alignedSeqB.Append("-");
                        i--;
                        break;

                }
            }

            string[] alignments = new string[2];
            //Da wir die Zeichen jeweils rechts an die Alignments angefügt haben,
            //rufen wir hier die Reverse-Methode um die Strings umzudrehen.
            alignments[0] = new string(alignedSeqA.ToString().Reverse().ToArray());
            alignments[1] = new string(alignedSeqB.ToString().Reverse().ToArray());

            //Der Rückgabewert ist ein array welches beide alignments enthält.
            return alignments;

        }

Im Anschluss noch folgende Überlegung, die Traceback-Matrix könnte effizienter sein, wenn man byte[] anstatt char[] verwendet. Dies würde zwar einige Sonderzeichen ausschließen uns aber gleichzeitig Speicherverbrauch sparen. Dies wäre insofern kein Problem da wir eigentlich nur drei Buchstaben brauchen.

Am besten kann man diesen Algorithmus nachvollziehen indem man vorgeht wie bei den meisten anderen, mit Stift und Papier 🙂

Hier haben wir uns das globale Alignment zweier Sequenzen angeguckt. Als nächstes bietet es sich an entweder mit dem lokalen Alignment (Smith-Waterman-Algorithmus) oder gleich dem Multiple-Sequence-Alignment weiter zu machen.

Bis bald.

Advertisements
Dieser Beitrag wurde unter .NET, Algorithmen, Beispiel, Beispiel, Bioinformatik, c#, global alignment, globales alignment, Mathematik, needleman-wunsch, Programmierung, Sequence alignmenet, Sequenzalignment abgelegt und mit , , , , , , , , , verschlagwortet. Setze ein Lesezeichen auf den Permalink.

9 Antworten zu Globales Sequenzalignment mit dem Needleman-Wunsch-Algorithmus in C#

  1. Pingback: Lokales Sequenzalignment mit dem Smith-Waterman Algorithmus in C# | RealityBites

  2. thortroll schreibt:

    Danke :). Gut und verständlich auch für Nicht-Biologen. Habe den Blog mal für die „nicht-Klausur-Zeit“ abgespeichert und gucke dann mal, was es hier noch so gibt. Scheint ja diverse Algorithmen in C-Derivaten zu geben :).

  3. MALMAL schreibt:

    hi, hast du bei deinem beispiel nicht Zeilen und Spalten vertauscht?

    • MALMAL schreibt:

      bzw. nicht nur Zeile und Spalte, sondern auch Links und oben 😉

      • plusnoir schreibt:

        Ist das eine Frage oder eine Feststellung? Der Code ist ewig her und ich habe ihn damals eigentlich getestet. ABer danke für den Hinweis, werde ihn bei nächsten gelegenheit nochmal überprüfen.

      • MALMAL schreibt:

        nicht der Code, aber deine Beispiele mit den Bildern: Hier:
        „Nehmen wir zunächst eine Matrix F(m+1 x n+1), wobei m die Länge des Strings a und n die Länge des Strings b ist. “
        aber du hast a auf den Spalten eingetragen, und b haste in den Zeilen eingetragen, und du sagst aber M(Zeile x Spale) aber das macht dann wenig Sinn weil m die Länge deiner Spalten sind und n die Länge der Zeilen und nicht umgekehrt…

      • plusnoir schreibt:

        Du hast absolut Recht. So habe ich das gemacht ja. Dennoch ist der Code aber an eben so ein Vorgehen angepasst. Oder siehst du hier ein Problem? Dieses Vorgehen zerschießt die Funktionalität nicht.

      • MALMAL schreibt:

        nein es kommt trotzdem das richtige raus, ich fands nur irgendwie verwirrend wenn man die Schritte nach den Formeln nachvollziehen will und dann immer Spalte mit Zeile vertauscht, habs aber auch auf Wikipedia so gesehen… naja ist nicht so wichtig glaube ich, und ich habs auch so verstanden 😉

Kommentar verfassen

Trage deine Daten unten ein oder klicke ein Icon um dich einzuloggen:

WordPress.com-Logo

Du kommentierst mit Deinem WordPress.com-Konto. Abmelden /  Ändern )

Google+ Foto

Du kommentierst mit Deinem Google+-Konto. Abmelden /  Ändern )

Twitter-Bild

Du kommentierst mit Deinem Twitter-Konto. Abmelden /  Ändern )

Facebook-Foto

Du kommentierst mit Deinem Facebook-Konto. Abmelden /  Ändern )

w

Verbinde mit %s