1-BIN-301, 2-AIN-501 Methods in Bioinformatics, 2023/24

Introduction · Rules · Tasks and dates · Materials · Moodle
Quizzes can be found in Moodle.
Homework assignments and journal club papers can be found in Tasks and dates.
Exam rules, example questions and syllabus
Groups for journal club have each their own group in Moodle.


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)