Lokales Sequenzalignment mit dem Smith-Waterman Algorithmus in C#

Wir haben festgestellt, dass der Needleman-Wunsch Algorithmus sich gut dazu eignet ein globales Alignement zweier Sequenzen zu berechnen. Dabei werden beide Sequenzen über ihre gesamte Länge aligniert. Das bedeutet je weniger sich die beiden Sequenzen ähneln, desto mehr Lücken (gaps) wird das resultierende Alignement aufweisen. Dieses Verfahren eignet sich besonders wenn beide Sequenzen ungefähr dieselbe Länge haben und möglicherweise ähnlich sind.

Was ist jedoch wenn wir nur eine teilweise Übereinstimmung der Sequenzen vermuten? Wenn wir annehmen, dass eine Übereinstimmung nur in kleineren Regionen einer längeren Sequenz vorkommt benutzen wir den Smith-Waterman Algorithmus.

Bezüglich der Implementierung unterscheiden sich die beiden Verfahren nur in wenigen Punkten. Auch der Smith-Waterman Algorithmus bedient sich der dynamischen Programmierung. Das bedeutet wir können die Berechnung eines lokalen Alignements in drei Schritte unterteilen, genau wie beim Needleman-Wunsch Algorithmus.

  1. Initialisierung einer Matrix
  2. Berechnung des Score der Matrix
  3. Traceback

Doch worin genau sich die Verfahren unterscheiden sieht man wohl am besten an einem Beispiel. Nehmen wir folgende zwei Sequenzen:

A=‘MISSISSIPPI‘

B=‘ISSP‘

Auch hier brauchen wir Gapkosten f, einen Match- und Mismatchscore als Eingabeparameter für den Algorithmus. Wenn man Elemente der Sequenzen vergleicht, liefert eine Bewertungsfunktion w(x, y) den Matchscore wenn die Elemente identisch sind und entsprechend den Mismatchscore wenn sie sich unterscheiden.

f = -1

MatchScore = 2

MismatchScore = -1

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.

Die Unterschiede beginnen direkt bei der Initialisierung der Matrix. Beim Needleman-Wunsch Algorithmus werden die erste Zeile mit i * (-1) initialisiert und die erste Spalte mit j * (-1). Beim Smith-Waterman hingegen werden einfach sämtliche Initialisierungsfelder auf 0 gesetzt.

initSW

Ein weiterer Unterschied zu einem globalen Alignement besteht in der Berechnung des Scores. Zwar liefern auch hier die jeweiligen Vorgängerzellen die Werte für die aktuelle Zelle, jedoch dürfen wir diesmal die Werte nicht kleiner als 0 werden lassen. Aus diesem Grund werden Werte in der Score Matrix wie folgt berechnet.

scoreSW

Nehmen wir die Zelle (1,1) hier wird das Element „M“ mit dem Element „I“ verglichen. Also gucken wir uns die Vorgängerzellen an.

(Diagonal) F(i-1, j-1) + w(x, j) = F(0, 0) + w(“M”,”I”) = 0 – 1 = -1

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

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

(Vierte Option)    F(i, j) = 0

Da wir uns bei der Berechnung stets für den größtmöglichen Wert entscheiden müssen nehmen wir in diesem Fall natürlich die vierte Option, nämlich 0.

firstCellSW

Wir werden auch in diesem Fall zeitgleich eine Matrix erstellen die uns beim Traceback helfen soll. Der Unterschied zu Needleman-Wunsch besteht darin, dass wir immer wenn der vierte Fall eintritt auch in die Traceback-Matrix eine 0 eintragen. Da auch bei ihrer Initialisierung 0 in die Zellen eingetragen wird, ist die Traceback-Matrix zu diesem Zeitpunkt identisch mit der Scoring-Matrix. So berechnen wir die beiden Matrizen und das Ergebnis sieht wie folgt aus.

scoreMatrixSW

Die Scoring Matrix.

Anders als bei den Needleman-Wunsch Algorithmus liegt der Score der Matrix in der Zelle mit dem höchsten Wert, in diesem Fall 7. Dieser bestimmt auch den Startpunkt des Traceback.

tracebackmatrixSW

Die Traceback Matrix.

Die Einträge in der Traceback Matrix weisen jeweils darauf hin welche Vorgängerzelle den aktuellen Wert geliefert hat. Dabei steht D für diagonal, L für Links und U für „Up“ also oben. Wir beginnen mit dem Score der Matrix und schreiben solange die Alignements bis wir auf eine Zelle mit dem Wert 0 treffen, diese bildet nämlich das Abbruchkriterium für das Traceback. Das mitschreiben der Alignements erfolgt ansonsten genau wie im Needleman-Wunsch Verfahren.

finalTracebackSW

Das Traceback, wobei eine 0 das Abbruchkriterium bildet.

alignmentSW

Das resultierende Alignement.

Wie man sieht wird im Gegensatz zum globalen Alignement nicht über die ganze Länge aligniert sondern nur auf einem begrenzten, also lokalen Bereich.

Auch hier werden wir die ScoreFunction-Methode brauchen. Darin hat sich aber nichts im Vergleich zum NW-Verfahren geändert.

//Bekommt die zu vergleichenden Buchstaben, den Match- und Mismatch-Score als Parameter
        private 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;
        }

Die Methode Align() hat sich dagegen an mehreren Stellen geändert, beispielsweise bei der Initialisierung. Außerdem benutzen wir jetzt eine Variable score[2] um uns die Zelle mit dem Höchsten Wert zu merken. Auch dürfen wir negative Werte nicht mehr zulassen. Außerdem kann die Traceback-Matrix jetzt auch eine 0 enthalten.

//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] = 0;
                tracebackMatrix[i, 0] = '0';
            }

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

            //Hier werden die Werte für die Scoring-Matrix berechnet. Zeitgleich befüllen
            //wir die Traceback-Matrix.
            #region Scoring

            //Diese Variable benutzen wir um uns die Zelle mit dem Höchsten Wert zu merken.
            //Diese wird uns beim Traceback als Startpunkt dienen.
            int[] score = new int[2];

            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(Math.Max(oben, Math.Max(links, diagonal)),0);                     //Hier wird geprüft ob die aktuelle Zelle als Score für die matrix taugt.                     if(matrix[i,j] > matrix[score[0],score[1]])
                    {
                        score[0] = i;
                        score[1] = j;
                    }

                    //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';
                    }
                    else if (matrix[i, j] == 0)
                    {
                        tracebackMatrix[i, j] = '0';
                    }
                }
            }
            #endregion

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

        }

Die einzige Veränderung in der Methode Traceback() ist der Score der mit als Parameter übergeben werden muss und als Startzelle fungiert.

//Eine TraceBack-Matrix und die beiden Sequenzen sind die Paramter für das Traceback.
        //Außerdem wird hier die Position der Score-zelle an die Methode übergeben.
        private static string[] TraceBack(char[,] tracebackMatrix,int[] score ,string sequenzA, string sequenzB)
        {
            //Lege die Startposition anhand des Scores fest
            int i = score[0];
            int j = score[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;

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

7 Antworten zu Lokales Sequenzalignment mit dem Smith-Waterman Algorithmus in C#

  1. pupsi schreibt:

    saucooler blog
    weiter so

  2. CB schreibt:

    Eine perfekte Erklärung, hat mir sehr beim Verständnis der Aligning-Algorithmen (auch Needleman-Wunsch) geholfen. Ab jetzt hat dieser Blog ein Lesezeichen.

  3. lucyinthesky42 schreibt:

    Die Erklärung des Smith Waterman Algorithmus ist super! Habe etliche Videos und Uni-PDFs durchgestöbert, aber erst die Seite hier hat mir das Vorgehen klar gemacht, danke dafür! 🙂

  4. Jens schreibt:

    Wirklich gute Erklärung. Und vor allem gut dargestellt mit gut gewähltem Beispiel. Top!
    Ein Frage: In der Scoring-Matrix zum Beispiel – müsste da F(3,2) nicht 4 sein? F(2,1) = 2 und wir haben für die Bewertungsfunktion w(‚S‘,’S‘) = 2.

    (Oder aber ich habe es doch nicht so gut verstanden, wie ich gehofft habe – und wäre für eine Erklärung dankbar :-))

    • plusnoir schreibt:

      Zuerstmal freut es mich natürlich sehr wenn ich dir helfen konnte. Zu deinem zweiten Punkt werde ich mir das natürlich noch (Pen & Paper) angucken. Auf den ersten Blick muss ich aber sagen, du hast recht, ich habe wohl einen Fehler gemacht! Setzen 1! Werde ich korrigieren. Oder habe ich einfach eine Kontrolle für den unaufmerksamen Leser eingebaut? Muahaha…

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