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ý.


CI09

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

Hľadanie motívov zadefinovaných pravdepodobnostnou maticou

  • Mame danych n sekvencii S=(S_{1}\dots S_{n}), kazda dlzky m, dlzku motivu L, nulova hypoteza q (frekvencie nukleotidov v genome)
  • Hladame motiv vo forme pravdepodobnostneho profilu dlzky L a jeho vyskyt v kazdej sekvencii
  • Nech W[a,i] je pravdepodobnost, ze na pozicii i motivu bude baza a, W cela matica
  • o_{i} je pozicia vyskytu v sekvencii S_{i}, O=(o_{1}\dots o_{n}) su vsetky vyskyty
  • \Pr(S|W,O) je jednoduchý súčin, kde pre pozície v oknách použijeme pravdepodobnosti z W, pre pozície mimo okna použijeme q
    • \Pr(S_{i}|W,o_{i})=\prod _{{j=1}}^{{L}}W[S_{i}[j+o_{i}-1],j]\prod _{{j=1}}^{{o_{i}-1}}q[S_{i}[j]]\prod _{{j=o_{i}+L}}^{m}q[S_{i}[j]]
    • \Pr(S|W,O)=\prod _{{i=1}}^{n}\Pr(S_{i}|W,o_{i})
  • Hľadáme W a O, ktoré maximalizujú tuto vierohodnosť Pr(S|W,O)
    • Nepozname efektivny algoritmus, ktory by vedel vzdy najst maximum
    • Dali by sa skusat vsetky moznosti O, pre dane O je najlepsie W frekvencie z dat
    • Naopak ak pozname W, vieme najst najlepsie O
      • v kazdej sekvencii i skusame vsetky pozicie o_{i} a zvolime tu, ktora ma najvyssiu hodnotu Pr(S_{i}|W,o_{i})

EM algoritmus

  • Iterativne zlepsuje W, pricom berie vsetky O vahovane podla ich pravdepodobnosti vzhladom na W z minuleho kola
  • Videli sme na prednaske, tu je trochu prepisany:
  • Inicializácia: priraď každej pozícii j v sekvencii S_{i} nejaké skóre p_{{i,j}}
  • Iteruj:
    • Spočítaj W zo všetkých možných výskytov v S_{1},\dots ,S_{k} váhovaných podľa p_{{i,j}}
    • Prepočítaj všetky skóre p_{{i,j}} tak, aby zodpovedali pomerom pravdepodobností výskytu W na pozícii j v S_{i}, t.j. p_{{i,j}} je umerne Pr(S_{i}|W,o_{i}=j), pricom hodnoty normalizujeme tak, aby sucet v riadku bol 1

Gibbsovo vzorkovanie (Gibbs sampling)

  • Inicializácia: Vezmi náhodné pozície výskytov O
  • Iteruj:
    • Spočítaj W z výskytov O
    • Vyber náhodne jednu sekvenciu S_{i}
    • Pre každú možnú pozíciu j v S_{i} spočítaj skóre p_{{i,j}} (ako v EM) výskytu W na tejto pozícii
    • Zvoľ o_{i} náhodne s váhovaním podľa s_{{i,j}}
  • Takto dostavame postupnost vzoriek O^{{(0)}},O^{{(1)}},....
  • Za sebou iduce vzorky sa podobaju (lisia sa len v jednej zlozke o_{i}) nie su teda nezavisle
  • Pre kazdu vzorku O^{{(t)}} najdeme najlepsie W^{{(t)}} a spocitame vierohodnost \Pr(S|W^{{(t)}},O^{{(t)}}). Nakoniec vyberieme O a W, kde bola vierohodnost najvyssia.
  • Tento algoritmus (s malymi obmenami) bol pouzity v clanku Lawrence, Charles E., et al. (1993) "Detecting subtle sequence signals: a Gibbs sampling strategy for multiple alignment." Science.
    • V clanku v kazdej iteracii maticu W rataju zo vsetkych sekvencii okrem S_{i}
    • Obcas robia krok, kde nahodne skusaju posunut vsetky vyskyty o jedna dolava alebo doprava
    • Tento algoritmus nie je uplne matematicky korektne Gibbsovo vzorkovanie (nema ani poradne zadefinovane rozdelenie, z ktoreho vzorkuje). Na spodku stranky pre informaciu uvadzame algoritmus Gibbsovho vzorkovanie pre hladanie motivov z ineho clanku.

Vzorkovanie z pravdepodobnostného modelu vo všeobecnosti

  • Majme pravdepodobnostny model, kde D su nejake pozorovane data a X nezname nahodne premenne (napr pre nas D su sekvencie S a X su vyskyty O, pripadne aj matica W)
  • mozeme hladat X pre ktore je vierohodnost Pr(D|X) najvyssia
  • alebo mozeme nahodne vzorkovat rozne X z Pr(X|D)

Pouzitie vzoriek

  • spomedzi ziskanych vzoriek zvolime tu, pre ktoru je vierohodnost Pr(D|X) najvacsia (iny pristup k maximalizovaniu vierohodnosti)
  • ale vzorky nam daju aj informaciu o tom, aka je velka neurcitost v odhade X
    • mozeme odhadovat stredne hodnoty a odchylky roznych velicin
    • napr. pri hladani motivov mozeme sledovat ako casto je ktora pozicia vyskytom motivu
  • generovat nezavisle vzorky z Pr(X|D) moze byt tazke
  • metoda Markov chain Monte Carlo (MCMC) generuje postupnost zavislych vzoriek X^{{(0)}},X^{{(1)}},\dots , konverguje v limite k cielovej distribucii Pr(X|D)
  • Gibbsovo vzorkovanie je specialnym pripadom MCMC

Markovove reťazce

  • Markovov reťazec je postupnosť náhodných premenných X^{{(0)}},X^{{(1)}},\dots , taká, že \Pr(X^{{(t)}}|X^{{(0)}},\dots ,X^{{(t-1)}})=\Pr(X^{{(t)}}|X^{{(t-1)}}), t.j. hodnota v čase t závisí len od hodnoty v čase t-1 a nie ďalších predchádzajúcich hodnôt.
  • Nás budú zaujímať homogénne Markovove reťazce, v ktorých \Pr(X^{{(t)}}|X^{{(t-1)}}) nezávisí od t.
  • Tiez nas zaujimaju len retazce v ktorych nahodne premenne X_{t} nadobudaju hodnoty z konecnej mnoziny (mozne hodnoty X^{{(t)}} nazyvame stavy)
    • Napriklad stavy A,C,G,T
    • V Gibbsovom vzorkovani pre motivy je stav konfiguracia premennych O (t.j. mame (m-L+1)^n stavov)
      • Vzorka v kroku t zavisi od vzorky v kroku t-1 (a lisi sa len v hodnote jedneho o_i)

Matica

  • Pravdepodobnosti prechodu medzi stavmi za jeden krok mozeme vyjadrit maticou pravdepodobnosti P, ktorej prvok p_{{x,y}} oznacuje pravdepodobnost prechodu zo stavu x do stavu y p_{{X,Y}}=\Pr(X_{t}=y|X_{{t-1}}=x)
    • Sucet kazdeho riadku je 1, cisla nezaporne
  • Ako p_{{x,y}}^{t} budeme oznacovat \Pr(X^{{(t)}}=y|X^{{(0)}}=x), tieto hodnoty dostaneme umocnenim matice P na t

Stacionarne rozdelenie

  • Rozdelenie \pi na mnozine stavov sa nazyva stacionarne pre Markovov retazec P, ak pre kazde j plati \sum _{{i}}\pi (i)p_{{i,j}}=\pi (j)\, (alebo v maticovej notacii \pi P=\pi )
  • Ak matica P splna urcite podmienky (je ergodicka), existuje pre nu prave jedno stacionarne rozdelenie \pi . Navyse pre kazde x a y plati \lim _{{t\to \infty }}p_{{x,y}}^{{t}}=\pi (y)\,

Priklady Markovovskych retazcov v bioinformatike

  • V HMM stavy tvoria Markovov retazec
  • Ine varianty: nekonecne stavove priestory (zlozitejsia teoria), spojity cas (videli sme pri evolucnych modeloch), retazce vyssieho radu, kde urcujeme \Pr(X_{t}|X_{{t-r}},\dots ,X_{{t-1}}) a pod.
  • Pouzitie v bioinformatike: charakterizacia nahodnych sekvencii (nulova hypoteza), pre DNA sa pouzivaju rady az do 5, lepsie ako nezavisle premenne

Ergodické Markovove reťazce

  • Vravime ze matica je ergodicka, ak P^{t} pre nejake t>0 ma vsetky polozky nenulove
  • Priklady neergodickych matic
1 0          0.5 0.5          0 1             0.5 0.5
0 1          0   1            1 0             1   0
nesuvisla    slabo suvisla    periodicka      ergodicka
  • V HMM stavy tvoria Markovov retazec; hladanie genov ergodicky stavovy priestor, profilove HMM nie

Markov chain Monte Carlo MCMC

  • Chceme generovať náhodné vzorky z nejakeho cieloveho rozdelenia \pi , ale toto rozdelenie je prilis zlozite.
  • Zostavime ergodicky Markovov retazec, ktoreho stacionarne rozdelenie je rozdelenie \pi , tak aby sme efektivne vedeli vzorkovat X^{{(t)}} ak vieme X^{{(t-1)}}.
  • Ak zacneme z lubovolneho bodu X^{{(0)}}, po urcitom case t rozdelenie X^{{(t)}} priblizne \pi
  • Ale za sebou iduce vzorky nie su nezavisle!
  • Vieme vsak odhadovat ocakavane hodnoty roznych velicin {\frac  {1}{t}}\sum _{{i=1}}^{t}f(X^{{(t)}}) konverguje k E_{\pi }[f(X)]

Gibbsovo vzorkovanie

  • Cielove rozdelenie \pi (X) je cez vektory dlzky n X=(x_{1},...x_{n})
  • V kazdom kroku vzorkujeme jednu zlozku vektora x_{i} z podmienenej pravdepodobnosti \Pr(x_{i}|x_{1},\dots ,x_{{i-1}},x_{{i+1}},\dots x_{n})
  • Ostatne hodnoty nechame rovnake ako v predchadzajucom kroku
  • Hodnotu i zvolime nahodne alebo periodicky striedame i=1,2,\dots ,n

Dôkaz správnosti Gibbsovho vzorkovania

  • Pozor! Gibbsovo vzorkovanie nie je vzdy ergodicke, ak niektore kombinacie hodnot maju nulovu pravdepodobnost!
  • Treba dokazat, ze ak je ergodicky, tak ma ako stacionarnu distribuciu nase zvolene \pi
  • Definicia: Vravime, ze matice P a rozdelenie \pi splnaju detailed balance, ak pre kazde stavy (dva vektory hodnot) x a y mame \pi (x)p_{{x,y}}=\pi (y)p_{{y,x}}
  • Lema: ak pre nejaky retazec P a nejake rozdelenie \pi plati detailed balance, \pi je stacionarna distribucia pre P
    • Dokaz: \sum _{x}\pi (x)p_{{x,y}}=\sum _{x}\pi (y)p_{{y,x}}=\pi (y)\sum _{x}p_{{y,x}}=\pi (y)
  • Lema: pre retazec Gibbsovo vzorkovania plati detailed balance vzhladom na cielove rozdelnie \pi
    • Dokaz: uvazujme dva za sebou iduce vektory hodnot x a y, ktore sa lisia v i-tej suradnici. Nech x_{{-i}} su hodnoty vsetkych ostatnych premennych okrem x_{i}
    • \pi (x)p_{{x,y}}=\pi (x)\Pr(y_{i}|x_{{-i}})=\Pr(x_{{-i}})\Pr(x_{i}|x_{{-i}})\Pr(y_{i}|x_{{-i}})=\pi (y)\Pr(x_{i}|x_{{-i}})=\pi (y)\Pr(x_{i}|y_{{-i}})=\pi (y)p_{{y,x}}

Poriadnejšie Gibbsovo vzorkovanie pre motívy

Uvedene pre zaujimavost - podla clanku Siddharthan R, Siggia ED, van Nimwegen E (December 2005). "PhyloGibbs: a Gibbs sampling motif finder that incorporates phylogeny". PLoS Comput. Biol. 1 (7): e67. doi:10.1371/journal.pcbi.0010067. PMID 16477324.

Pravdepodobnostny model

  • Rozsirime model, aby aj O a W boli nahodne premenne, takze mame rozdelenie Pr(S,W,O)
    • Potom chceme vzorkovat z Pr(O|S) (marginalizujeme cez vsetky hodnoty W)
  • Vygeneruje sa nahodne matica pravdepodobnosti W (napr z roznomernej distribucie cez vsetky matice)
  • V kazdej sekvencii i sa zvoli okno o_{i} dlzky L (rovnomerne z m-L+1 moznosti)
  • V okne sa generuje sekvencia podla profilu W a mimo okna sa generuje sekvencia z nulovej hypotezy (ako predtym)

Gibbsovo vzorkovanie

  • Mame dane S, vzorkujeme O (O^{{(0)}},O^{{(1)}},\dots ) (ak treba, z O^{{(t)}} mozeme zostavit maticu W^{{(t)}})
    • zacni s nahodnymi oknami O^{{(0)}}
    • v kroku t+1 zvol jednu sekvenciu i a pre vsetky pozicie o'_{i} spocitaj \Pr(o'_{i}|O_{{-i}}^{{(t)}},S) (kde O_{{-i}}=o_{1}\dots o_{{i-1}}o_{{i+1}}\dots o_{n}, t.j. všetky pozície výskytov okrem i-tej).
    • nahodne zvol jedno o'_{i} umerne k tymto pravdepodobnostiam
    • O^{{(t+1)}} dostaneme z O^{{(t)}} vymenou pozicie v sekvencii i za prave zvolenu
    • opakuj vela krat
  • Konverguje k cielovemu rozdeleniu \Pr(O|S), ale vzorky nie su nezavisle
  • Dalsie mozne kroky vo vzorkovani: posun vsetky okna o konstantu vlavo alebo vpravo
  • Dalsie moznosti rozsirenia modelu/algoritmu: pridaj rozdelenie cez L a nahodne zvacsuj/zmensuj L, dovol vynechat motiv v niektorych sekvenciach, hladaj viac motivov naraz,...

Ako spocitat \Pr(o_{i}|O_{{-i}},S)?

  • nezaujimaju nas normalizacne konstanty, lahko znormalizujeme scitanim cez vsetky o'_{i}
  • \Pr(o_{i}|O_{{-i}},S)=\Pr(O|S)/\Pr(O_{{-i}}|S), ale menovatel konstanta
  • \Pr(O|S)=\Pr(S|O)\Pr(O)/\Pr(S), kde \Pr(S)=\sum _{{O'}}\Pr(S|O')\Pr(O')
  • Menovatel nas nezaujima (normalizacna konstanta)
  • \Pr(O) je tiez konstanta (rovnomerne rozdelenie pozicii okien)
  • Teda mame \Pr(o_{i}|O_{{-i}},S) je umerne \Pr(S|O)
  • Lahko vieme spocitat \Pr(S|W,O), potrebujeme "zrusit" W, da sa spocitat vzorec...
  • Skusame vsetky mozne hodnoty o'_{i}, pocitame pravdepodobnost \Pr(S|O), vzorkujeme umerne k tomu

Dalsie detaily vypoctu \Pr(S|O):

  • Nech S_{o} su len sekvencie v oknach a S_{{-o}} mimo okien. Mame \Pr(S|O)=\Pr(S_{o}|O)\Pr(S_{{-o}}|O)
  • \Pr(S_{{-o}}|O) lahko spocitame (nezavisi od W)
  • \Pr(S_{o}|O)=\int \Pr(S_{o}|O,W)\Pr(W)dW kde integral ide cez hodnoty, kde w_{{a,i}}\geq 0 a \sum _{a}w_{{a,i}}=1\,
  • \Pr(W) je konstanta (rovnomerne rozdelenie; nejde o pravdepodobnost ale hustotu), \Pr(S_{o}|O,W)=\prod _{{i=1}}^{L}\prod _{a}(w_{{a,i}})^{{n_{{a,i}}}}, kde n_{{a,i}} je pocet vyskytov bazy a na pozicii i v oknach o_{1}\dots o_{n}
  • \Pr(S_{o}|O)=\prod _{{i=1}}^{L}3!/(n+3)!\prod _{a}n_{{a,i}}! (bez dokazu)