1-BIN-301, 2-AIN-501 Metódy v bioinformatike, ZS 2018/19

Úvod · Pravidlá · Termíny a zadania · Prednášky a poznámky · Facebook (oznamy a diskusie) (návod a pravidlá)
Zadania domácich úloh a články na journal club nájdete v časti Termíny a zadania.
Pozrite si ukážkové príklady na skúšku.
Rozpis skupín pre journal club je zverejnený.


CI04

Z MBI
Prejsť na: navigácia, hľadanie

Opakovanie dynamického programovania pre globálne zarovnanie

Uvažujme napríklad skórovanie zhoda +1, nezhoda -1, medzera -1 a vstupné sekvencie X=x_{1}\dots x_{m} a Y=y_{1}\dots y_{n}. Nech s(x,y) je skóre písmen x a y, t.j. 1 ak sa zhodujú a -1 ak nie. Máme rekurenciu:

A[i,j]=\max \left\{A[i-1,j-1]+s(x_{i},y_{j}),A[i-1,j]-1,A[i,j-1]\right\}
  • Ako presne by sme implementovali?
  • Ako spocitame maticu spatnych sipok B?
  • Aka je casova a pamatova zlozitost?

Reprezentácia pomocou grafu

Takéto dynamické programovanie vieme reprezentovať vo forme acyklického orientovaného grafu:

  • vrchol (i,j) pre každé 0\leq i\leq m,0\leq j\leq m, t.j. pre každé políčko dyn. prog. tabuľky
  • hrana z (i-1,j-1) do (i,j) s cenou s(x_{i},y_{j})
  • hrana z (i-1,j) do (i,j) s cenou -1
  • hrana z (i,j-1) do (i,j) s cenou -1
  • súčet súradníc na každej hrane rastie, graf teda nemôže obsahovať cyklus, je acyklický
  • každá cesta z (0,0) do (m,n) zodpovedá zarovnaniu, jej cena je cenou zarovnania (každá hrana jeden stĺpec)
  • optimálne zarovnanie teda zodpovedá ceste s maximálnou cenou

Krátka vsuvka o acyklických orientovaných grafoch

  • Mame dany acyklicky orientovany graf s ohodnotenymi hranami a startovaci vrchol s, koncovi vrchol t a chceme najst cestu s max. cenou z s do t.
  • Hladanie cesty s maximalnou cenou je vo vseobecnosti NP-tazke (podobne na Hamiltonovsku cestu)
  • V acyklickom grafe to vsak vieme riesit efektivne
  • Najskor si graf zotriedime topologicky, t.j. usporiadame vrcholy tak, aby kazda hrana isla z vrcholu z mensim cislom do vrcholu s vacsim cislom. To sa da modifikaciou prehladavania do hlbky v case O(|V|+|E|)
  • Potom pocitame dynamickym programovanim, kde A[u] je dlzka najdlhsej cesty z s do u: A[u]=\max _{{v:v\rightarrow u\in E}}A[v]+c(v\rightarrow u)

pricom na zaciatku nastavime A[s]=0 a na konci mame cenu cesty v A[t].

  • Cas vypoctu je O(|V|+|E|)
  • Vsimnime si, ze tiez dostaneme najdlhsie cesty z s do vsetkych vrcholov.

Ak tento algoritmus nasadime na graf pre globalne zarovnanie, dostavame presne nasu rekurenciu (topologicke triedenie mozno vynechat - poradie zhora dole a zlava doprava je topologicky utriedene). Vyhoda je, ze mozeme modifikaciou grafu ziskavat riesenia roznych pribuznych problemov bez toho, aby sme vzdy vymyslali novu rekurenciu.

Lokálne zarovnanie

  • Zarovnanie moze zacat a skoncit hocikde v matici
  • Pridaj startovaci vrchol s, koncovy vrchol t
  • Pridaj hrany s->(i,j) a (i,j)->t s cenou 0 pre kazde (i,j)
  • Opat ekvivalentne s rekurenciou z prednasky

Variant: chceme zarovnat cely retazec X k nejake casti retazca Y (napr. mapovanie sekvenovacich readov na genom)

  • Iba zmenime hrany z s a hrany do t (ako?)

Afínne skóre medzier

  • Napr. otvorenie medzery o=-3, pokracovanie medzery e=-1
A  -  -  -  T  C  G
A  C  G  C  T  C  C
1 -3 -1 -1  1  1  -1


Nesprávne riešenie pomocou dynamického programovania

Pouzijeme bezne dynamicke programovanie pre globalne zarovnanie, ale v rekurencii zmenime vypocet penalty za medzeru:

A[i,j]=\max \left\{A[i-1,j-1]+s(x_{i},y_{j}),A[i-1,j]+c(i-1,j,hore),A[i,j-1]+c(i,j-1,vlavo)\right\}
  • c(i,j,s) = o, ak v policku A[i,j] mame sipku s
  • c(i,j,s) = e, ak v policku A[i,j] mame inu sipku

Preco toto riesenie nefunguje?

  • Co ak pre policko (i,j) je viac rovnako dobrych rieseni s roznymi sipkami?
  • Co ak pre policko (i,j) je najlepsie riesenie so sipkou napr. sikmo, ale druhe najlepsie je len 1 horsie a ma sipku hore?

Toto je obvykla chyba pri dynamickom programovani:

  • aby bolo dynamicke programovanie spravne, musi platit, ze optimalne riesenie vacsieho podproblemu musi obsahovat optimalne riesenie mensieho podproblemu

Správne riešenie pomocou dynamického programovananie

Riesenie 1:

  • Pridame hrany pre cele suvisle useky medzier so spravnou cenou
  • (i,j)->(i,k) s cenou o+(k-j)e
  • (i,j)->(k,j) s cenou o+(k-i)e
  • Cas O(mn(m+n)), t.j. kubicky
  • pozor, mame aj cesty, ktore nezopodvedaju ziadnemu spravnemu skore, napr. (i.j)->(i+1,j)->(i+2,j) ma cenou 2o, ale ma mat o+e. Nastastie hrana (i,j)->(i+2,j) ma vyssiu cenu, takze ta dlhsia cesta sa nepouzije.

Riesenie 2:

  • ztrojnasobime kazdy vrchol (i,j)_{u},(i,j)_{v},(i,j)_{z}
  • v indexe si pamatame, odkial sme do (i,j) prisli (u=uhlopriecne, v=vodorovne, z=zvislo)
  • ak ideme napr. z (i,j-1)_{v} do (i,j)_{v}, pokracujeme v uz existujucej medzere, takze skore je e
  • ak ideme napr. z (i,j-1)_{u} do (i,j)_{v}, zaciname novu medzeru, takze skore je o
  • ake vsetky hrany teda mozeme mat? Kolko je spolu v grafe hran a vrcholov a aka je zlozitost algoritmu?

Lineárna pamäť: Hirshbergov algoritmus 1975

  • Klasicke dynamicke programovanie potrebuje cas O(nm)
  • Trivialna implementacia tiez pouzije pamat O(mn) - uklada si celu maticu A, pripadne maticu B so sipkami naspat
  • Na vypocet matice A nam z stacia dva riadky tejto matice: riadok i pocitam len pomocou riadku i-1, starsie viem zahodit
  • Ale ak chcem aj vy[isat zarovnanie, stale potrebujem pamat O(mn) na maticu sipok B
  • Hirschbergov algoritmus znizi pamat na O(m+n), zhruba zdvojnasobi cas (stale O(mn))
  • Prejdeme celú maticu a spočítame maticu A. Zapamätáme si, kde moja cesta prejde cez stredný riadok matice
    • Nech B_k[i,j] je najväčší index v riadku k, cez ktorý prechádza najkratšia cesta z (0,0) do (i,j)
  • Ako vieme B_k[i,j] spočítať?
    • ak A[i,j] = A[i-1,j-1]+w(S[i],T[j])$, potom B_k[i,j]=B_k[i-1,j-1].
    • ak A[i,j]=A[i-1,j]+1, potom B_k[i,j]=B_k[i-1,j].
    • ak A[i,j]=A[i,j-1]+1, potom B_k[i,j]=B_k[i,j-1]
    • Toto platí, ak i > k. Pre i=k nastavíme B_k[i,j]=j
  • Ak už poznáme A[i-1,*] a B_k[i-1,*], vieme spočítať A[i,*] a B_k[i,*].
    • Stacia nam teda iba dva riadky matice A a B_k
  • Nech k'=B_k[m,n]. Potom v optimálnom zarovnaní sa S[1..k] zarovná s T[1..k'] a S[k+1..m] s T[k'+1..n].
    • Toto použijeme na rekurzívny algoritmus na výpočet zarovnania:
optA(l1, r1, l2, r2) { // align S[l1..r1] and T[l2..r2]
    if(r1-l1 <= 1 ||  r2-l2 <=1) 
        solve using dynamic programming
    else {
        k=(r-l+1)/2;
        for (i=0; i<=k; i++) 
           compute A[i,*] from A[i-1,*]
        for (i=k+1; i<=r-l+1; i++) 
           compute A[i,*], B_k[i,*] from A[i-1,*], B_k[i-1,*]
        k2=B_k[r1-l1-1,r2-l2-1];
        optA(l1, l1+k-1, l2, l2+k2-1); 
        optA(l1+k, r2, l2+k2, r2); 
    }
}

Casova zlozitost:

  • Označme si N=nm (súčin dĺžky dvoch daných reťazcov).
  • Na hornej úrovni rekurzie spúšťame dynamické programovanie pre celú maticu -- čas bude $cN$.
  • Na druhej urovni mame dva podproblemy, velkosti N1 a N2, pricom N1+N2<=0.5*N (z kazdeho stlpca matice A najviac polovica riadkov pocitana znova)
  • Na tretej urovni mame 4 podproblemy N11, N12, N21, N22, pricom N11+N12 <= 0.5*N1 a N21+N22 <= 0.5*N2 a teda celkovy sucet podproblemov na druhej urvni je najviac N/4.

Na stvrtej urovni je sucet podproblemov najviac N/8 atd, Dostavame geometricky rad cN+cn/2+cN/4+... ktoreho sucet je 2cN

Vypísanie všetkých najlepších riešení

  • Namiesto jednej spatnej sipky si pamatame vsetky, ktore v danom A[i,j] viedli k maximalnej cene
  • Potom mozeme rekurzivne prehladavat a vypisovat vsetky cesty z (m,n) do (0,0) ktore pozostavaju iba zo zapamatanych hran
  • Cas na vypisanie jednej cesty je polynomialny, ale ciest moze byt exponencialne vela!
  • Mozno namiesto toho chceme len pocet takych ciest, alebo vsetky dvojice pismen, ktore mozu byt spolu zarovnane v niektorom optimalnom zarovnani