Needleman-Wunschov algorithm

Iz MaFiRaWiki

Ta članek ali del članka je v delu. Veseli bomo, če ga boste dopolnili in popravili.

Kaj pomeni to opozorilo?

Vsebina

Algoritem

Needleman–Wunschov algoritem primerja dva niza in ugotovi, kako podobna sta si. Običajno se uporablja v bioinformatiki za primerjavo dveh nizov aminokislin ali nukleotidov. Algoritem sta leta 1970 predlagala Saul Needleman in Christian Wunsch v članku A general method applicable to the search for similarities in the amino acid sequence of two proteins [1].

Needleman–Wunsch algoritem je primer dinamičnega programiranja, ki nam vrne primerjavo nizov z maksimalno oceno in je prvi algoritem, ki s pomočjo dinamičnega programiranja primerja dva biološka niza.

Ocene za primerjavo aminokislin je definirana v matriki podobnosti S. Tako člen Si,j predstavlja oceno za i-to in j-to aminokislino. Pri iskanju maksimalne ocene si pomagamo tudi s penali d za manjkajočo aminokislino.

Če pogledamo na primeru, kjer je matrika podobnosti enaka

    A  G  C  T 
 A 10 -1 -3 -4 
 G -1  7 -5 -3 
 C -3 -5  9  0 
 T -4 -3  0  8 

potem ima poravnava

 AGACTAGTTAC
 CGA---GACGT 

s penalom d = − 5, naslednjo oceno

 S_{A,C} + S_{G,G} + S_{A,A} +  3\cdot d + S_{G,G} + S_{T,A} + S_{T,C} + S_{A,G} + S_{C,T}=  (-3) + 7 + 10 - 3 \cdot 5 + 7  + (-4) + 0 + (-1) + 0 = 1

Da lahko poiščemo poravnavo z maksimalno oceno, moram najprej ustvariti dvodimenzionalno tabelo (matriko). To tabelo običajno imenujemo F matrika in člen (i,j) običajno označimo Fi,j. V matriki F mora za vsa znak v nizu A obstajati stolpec in za vsak znak v nizu B obstajati vrstica. Iz tega sledi, da za izvajanje algoritma potrebujemo O(nm) časa in O(nm) pomnilnika. Možna je tudi izboljšava algoritma za katero potrebujemo samo O(n + m) pomnilnika, vendar se pri temu poveča čas izvajanja.

Med izvajanjem algoritma bo vsak člen Fi,j predstavljal optimalno oceno za poravnavo prvih i znakov niza A in prvih j znakov niza B.

Psevdo koda za izračun matrike F (indeksi se v tabeli in nizih začnejo z 0):


 for i=0 to length(A)-1
   F(i,0) <- d*i
 for j=0 to length(B)-1
   F(0,j) <- d*j
 for i=1 to length(A)
   for j = 1 to length(B)
     Choice1 <- F(i-1,j-1) + S(A(i), B(j))
     Choice2 <- F(i-1, j) + d
     Choice3 <- F(i, j-1) + d
     F(i,j) <- max(Choice1, Choice2, Choice3)

Ko je matrika F izračunana, se v matriki v spodnjem desnem kotu nahaja maksimalna ocena za vsako primerjavo niza. Da izračunamo katera poravnava dobi to oceno, potem pogledamo vrednost v spodnjem desnem kotu in primerjamo z vrednostjo, ki se nahajajo levo, vrednostjo, ki se nahaj zgoraj, in vrednostjo, ki se nahaja levo zgoraj. Če je vrednost enaka diagonalni vrednosti, potem sta znaka poravnana. Če je vrednost enaka zgornjemu sosedu, potem v niz A vstavimo oznako za manjkajoči znak in če je vrednost enaka levemu sosedu potem v niz B vstavimo oznako za manjkajoči znak.

Psevdo koda za poravnavo z maksimalno oceno

  AlignmentA <- ""
  AlignmentB <- ""
  i <- length(A) - 1
  j <- length(B) - 1
  while (i > 0 AND j > 0)
  {
  	Score <- F(i,j)
  	ScoreDiag <- F(i - 1, j - 1)
  	ScoreUp <- F(i, j - 1)
  	ScoreLeft <- F(i - 1, j)
  	if (Score == ScoreDiag + S(A(i), B(j)))
  	{
  	  AlignmentA <- A(i) + AlignmentA
  	  AlignmentB <- B(j) + AlignmentB
  	  i <- i - 1
  	  j <- j - 1
  	}
  	else if (Score == ScoreLeft + d)
  	{
  	  AlignmentA <- A(i) + AlignmentA
  	  AlignmentB <- "-" + AlignmentB
  	  i <- i - 1
  	}
  	otherwise (Score == ScoreUp + d)
  	{
  	  AlignmentA <- "-" + AlignmentA
  	  AlignmentB <- B(j) + AlignmentB
  	  j <- j - 1
  	}
  }
  while (i > 0)
  {
  	AlignmentA <- A(i) + AlignmentA
  	AlignmentB <- "-" + AlignmentB
  	i <- i - 1
  }  
  while (j > 0)
  {
  	AlignmentA <- "-" + AlignmentA
  	AlignmentB <- B(j) + AlignmentB
  	j <- j - 1
  }

Uporaba algoritma na konkretnih podatkih

V tem primeru bomo poiskali poravnavo za sledeča dva niza aminokislin

 GAATTCAGTTA 
 GAATCGA 


Matrika podobnosti S za ta primer je definirana kot

  • Si,j = 2 če je i-ta aminokislina iz prvega niza enaka j-ti aminokislini iz drugega niza.
  • Si,j = − 1 sicer

oziroma če ponazorimo z matriko

    A  G  T  C 
 A  2 -1 -1 -1 
 G -1  2 -1 -1 
 T -1 -1  2 -1 
 C -1 -1 -1  2 


Penal za manjkajočo aminokislino je enak d = − 2.

Inicijalizacija matrike F

Označimo z M velikost prvega niza in z N velikost drugega niza.

Najprej ustvarimo matriko z M + 1 stolpci in N + 1 vrsticami (indeksi v tabeli se začne z 0). Nato izračunamo člene za ničto vrstico in ničti stolpec, kjer je vrednost členov F0,j enaka d * j in vrednost členov Fi,0 enaka d * i.

Inicijalizirana matrika.

Preostali členi v matriki F

Da izračunamo preostale člene Fi,j, moramo za vsak i,j poznati vrednost členov, ki se nahajo levo, zgoraj in diagonalno levo zgoraj od člena Fi,j, . Torej potrebujemo vrednosti Fi,j − 1, Fi − 1,j in Fi − 1,j − 1.

Za vsak člen Fi,j priredimo vrednost, ki je enaka

 F{i,j} = \max      \begin{bmatrix}     F_{i-1,j-1} + S_{i,j} \\     F_{i-1,j} + d \\     F_{i,j-1} + d \\     \end{bmatrix} 

Če uporabimo to informacijo, lahko izračunamo vrednost F1,1. Torej izračunati moramo

F_{1,1} = \max      \begin{bmatrix}     F_{1-1,1-1} + S_{G,G}  \\     F_{1-1,1} + d     \\     F_{1,1-1} + d      \\     \end{bmatrix}      = \max      \begin{bmatrix}        0 +   2 \\      -2 + (-2) \\      -2 + (-2) \\     \end{bmatrix}

Tako F1,1 dobi vrednost 2 (maksimum smo dobili s pomočjo diagonalnega soseda).

Na sliki predstavlja puščica smer soseda s pomočjo katerega smo dobili maksimum.

Sedaj izračunamo vrednost F2,1 in vstavimo vrednost v matriko.

F_{2,1} = \max      \begin{bmatrix}     F_{2-1,1-1} + S_{A,G}  \\     F_{2-1,1} + d      \\     F_{2,1-1} + d      \\     \end{bmatrix}      = \max      \begin{bmatrix}     -2 + (-1) \\     2 + (-2)  \\     -4 + (-2) \\     \end{bmatrix}

Tu smo maksimum dobili s pomočjo zgornjega soseda.

Na enak način izračunamo tudi vrednost za F3,1.

F_{3,1} = \max      \begin{bmatrix}     F_{3-1,1-1} + S_{A,G} \\     F_{3-1,1} + d     \\      F_{2,1-1} + d     \\      \end{bmatrix}      = \max      \begin{bmatrix}     -4 + (-1) \\     0 + (-2)  \\     -6 + (-2) \\     \end{bmatrix}

Tudi tokrat smo maksimum dobili s pomočjo zgornjega soseda.

V določenih primerih lahko dobimo enako maksimalno vrednost od dveh ali vseh treh sosedov na enkrat. Na sliki smo označili vse sosede s pomočjo katerih dobimo maksimalno vrednost.

Maksimalno vredost so dobili od zgornjegasoseda in diagonalnega soseda.

Če označimo še smeri v ničti vrstici in ničtem stolpcu ter izračunamo vse člene, potem ima matrika F na koncu sledeče vrednosti.


Končna oblika matrike F.

Iskanje poravnav

Ko smo izračunali vse vrednosti matrike F, smo dobili tudi maksimalno oceno za to poravnavo. Ta ocena se nahaja v spodnjem desnem kotu in je v našem primeru enaka 6.

Z označevanjem 'predhodnikov' Fi,j-a, smo si olajšali iskanje vseh poravnav s maksimalno oceno. Vse kar moram narediti je, da se postavimo v spodnji desni kot in sledimo puščicam, ki na pove kateri člen moram izbrati na naslednjem koraku.

Začnemo v spodnjem desnem kotu.

Diagonalna puščica nam pove, da sta aminokislini poravnani. Zato ju samo prepišemo in sledimo naslednji puščici. Zato je poravnava do sedaj enaka

 A 
 A 

Slika 7: ačnemo v spodnjem desnem kotu.

V naslednjem koraku vidimo, da je puščica obrnjena v levo, kar pomeni da moramo v drugi niz označiti, da manjka aminokislina. To bomo označili z _ (podčrtajem). Podoben razmislek velja, če imamo opravka z navzgor obrnjeno puščico. V tem primeru označimo manjkajočo aminokislino v prvem nizu.

 TA 
 _A 


Če preskočimo nekaj korakov dobimo sledečo sliko

Po nekaj korakih.

kar nam da sledečo poravnavo

 CAGTTA 
 C_G__A 

Ko pridemo do člena, ki ima dve ali tri puščice, lahko izberemo katero poravnavo bomo uporabili. Tako stanje nam da več različnih poravnav, ki pa imajo enako maksimalno oceno. V našem primeru imamo samo dve taki poravnavi. Prva je

 GAATTCAGTTA 
 GAAT_C_G__A 

Prva možna poravnava

in druga poravana

 GAATTCAGTTA 
 GAA_TC_G__A 

Druga možna poravnava.

Da smo lažje prikazali delovanje algoritma, smo v ta namen risali puščice. Tak način dela nam sicer pomaga, da hitreje pridemo do iskane rešitve toda za to smo morali zavzeti pomnilnik velikosti O(mn). Boljši način je, če uporabimo predlagani algoritem za iskanje poravnave.

Osebna orodja