From Tomasz Żok

Main: ProgramowanieDynamiczneWBioinformatyceStrukturalnej

Wstęp

• Programowanie dynamiczne to algorytmiczny sposób rozwiązywania problemów
• Opiera się na podziale problemu na mniejsze części łatwiejsze do rozwiązania
• W przeciwieństwie do metody dziel-i-rządź mniejsze podproblemy nie są rozłączne

Algorytm Needlemana-Wunscha

• Służy do znalezienia uliniowienia (ang. alignment) różnej długości sekwencji s i t
• Rezultatem działania algorytmu są sekwencje s' oraz t' , które spełniają następujące cechy:
  • Są równej długości
  • Oprócz alfabetu używanego w oryginalnych sekwencjach s i t mogą zawierać specjalny znak pusty (ang. gap) oznaczany myślnikiem '-'
  • Na każdej pozycji s' i t' para znaków może być:
    • równa np. AA (ang. match )
    • różna np. AU (ang. mismatch )
    • zawierająca znak pusty np. -U (ang. indel , insertion/deletion)
• Istnieją różne schematy punktowania poszczególnych sytuacji; dla uproszczenia przyjmijmy, że match = +1, mismatch = -1, indel = -1
• Programowanie dynamiczne polega na tym, żeby dla coraz większych podciągów s i t podejmować decyzje o przyznawaniu match , mismatch lub indel

Algorytm krok po kroku:

• Przygotuj macierz w taki sposób, żeby każdy znak sekwencji s był w osobnej kolumnie, a każdy znak t w osobnym wierszu; dodaj na początku jeden dodatkowy wierz i kolumnę

• Wypełnij dodatkowy wiersz i kolumnę malejącymi wartościami

• Od teraz macierz można wypełniać tylko dla tych pól, które mają znaną wartość dla sąsiada z lewej, z góry i po ukosie między nimi
• Każde pole macierzy wypełniamy wartością maksymalną z:
  • Komórki po lewej minus 1 ( indel )
  • Komórki na górze minus 1 ( indel )
  • Komórki po ukosie minus 1 jeżeli jest mismatch
  • Komórki po ukosie plus 1 jeżeli jest match
• Zatem wypełniamy cały pierwszy wiersz. Na przecięciu P i P pojawia się 1, bo mamy match , a po ukosie wartość 0. Można zauważyć, że wszystkie dalsze wartości w wierszu stosują krok "komórka po lewej minus 1 (indel)"

• Analogicznie wypełniamy pozostałą część macierzy

• Kluczowe jest teraz odczytanie wyniku z macierzy. Musimy prześledzić wstecz (ang. traceback) skąd wzięła się każda wartość zaczynając od ostatniego pola w macierzy w prawym dolnym rogu
• Śledząc wstecz musimy sprawdzić czy aktualna wartość powstała z operacji indel od komórki z góry lub operacji indel od komórki z lewej, a może przez match lub mismatch
• Czasami jest to niejednoznaczne. Już w pierwszym kroku wstecz nie wiadomo czy wybrać  indel z komórki u góry czy mismatch z komórki po ukosie. W praktyce można przyjąć stałą kolejność sprawdzanych wyborów np. najpierw indel z góry, potem indel z lewej i na końcu mismatch
• Macierz po wykonaniu śledzenia wstecz wygląda następująco

• Aby wypisać wynik należy przy każdym kroku śledzenia tworzyć sekwencje s' oraz t' :
  1. W górę, zatem s' = -, t' = N
  2. W górę, zatem s' = --, t' = IN
  3. Na ukos, zatem s' = Y--, t' = EIN
  4. Na ukos, zatem s' = TY--, t' = TEIN
  5. Na ukos, zatem s' = TTY--, t' = TTEIN
  6. W lewo, zatem s' = ETTY--, t' = -TTEIN
  7. Na ukos, zatem s' = RETTY--, t' = R-TTEIN
  8. Na ukos, zatem s' = PRETTY--, t' = PR-TTEIN
• Ostateczny alignment

$ echo -e 'PRETTY\nPRTTEIN' | ./needleman-wunsch.py
PRETTY--
PR-TTEIN

Algorytm Nussinov

• W 1980 Ruth Nussinov podała algorytm przewidujący strukturę drugorzędową RNA
• Algorytm ten wykorzystuje programowanie dynamiczne by dla zadanej sekwencji znaleźć maksymalną liczbę par kanonicznych, które da się przypisać

Algorytm krok po kroku

• Budujemy macierz kwadratową gdzie wiersze i kolumny reprezentują pojedyncze znaki sekwencji

• Wypełniamy przekątną oraz jej przesunięcie w dół zerami

• W tym algorytmie wypełniamy macierz podążając od przekątnej aż do prawego górnego rogu
• Dla każdej komórki M[i,j] macierzy wyznaczamy maksimum z:
  • M[i + 1][j - 1] + s(i, j) , gdzie s(i, j) to 1 gdy na przecięciu macierzy mamy parę kanoniczną G-C lub A-U i 0 w innym wypadku
  • M[i][k] + M[k + 1][j] dla każdego i ≤ k < j
• Wizualnie prezentuje się to tak, że aby wyznaczyć ostateczny wynik, musimy albo wybrać z albo sumy pól oznaczonych jako t lub u lub v lub w lub x lub y. Oczywiście w toku rozwiązywania problemu wszystkie puste pola macierzy zostaną wypełnione.

• Wypełniamy zatem pierwszą serię pól. Na początku tylko podstawowe dwie pary A-U możemy dopisać

• Kolejna iteracja. W pierwszym wierszu ( i = 1, j = 3 ) mamy G-U, którego nie zaliczamy jako pary, ale dla k = 1 , mamy M[1,1] + M[2,3] = 0 + 1 = 1 , dlatego w macierzy wpisujemy 1

• Kolejna iteracja. W drugim wierszu ( i = 2, j = 5 ) mamy wartość 2, ponieważ zachodzi M[2][3] + M[4][5] = 1 + 1 = 2

• Kolejna iteracja. W pierwszym wierszu ( i = 1, j = 5 ) pojawia się 2, ponieważ dla k = 3 zachodzi M[1][3] + M[4][5] = 1 + 1 = 2

• Kolejna iteracja

• Ostatnia iteracja

• Wynik odczytujemy z lewego górnego narożnika. Jest to maksymalna liczba par, jakie mogą się znaleźć w strukturze o zadanej sekwencji
• Śledzenie wsteczne ponownie wymaga by na każdym kroku odczytać skąd wzięła się dana wartość w macierzy
  1. Wartość  M[1][7] to M[1][6] + M[7][7] (3 + 0)
  2. Wartość  M[1][6] to M[2][5] + 1 (para G-C)
  3. Wartość M[2][5] to M[2][3] + M[4][5] (1 + 1)
  4. Wartość M[2][3] to M[3][2] + 1 (para A-U)
  5. Wartość  M[4][5] to M[5][4] + 1 (para A-U)
• A zatem zestaw par to (1,6), (2,3) i (4,5), a struktura dot-bracket:

   
  GAUUACA
  (()()). 

$ echo GAUUACA | ./nussinov.py 
GAUUACA
(()()).