MBI 2010/2011

Počítanie fylogenetických stromov, cvičenia pre informatikov

• Ako definujeme strom v teorii grafov? suvisly acyklicky neorientovany graf
• Strom s n vrcholmi ma n-1 hran
• Nezakoreneny binarny fylogeneticky strom: neorientovany suvisl acyklicky graf, v listoch sucasne druhy, vsetky vnutorne vrcholy stupna 3
• Zakoreneny binarny fylogeneticky strom: vsetky hrany orientujeme od korena smerom k listom, kazdy vnutorny vrchol ma dve deti
• Niekedy uvazujeme aj nebinarne stromy, v ktorych mame vnutorne vrcholy vyssieho stupna
• Zakoreneny binarny strom s n listami ma n-1 vnutornch vrcholov, teda 2n-2 hran
• Nezakoreneny binarny strom s n listami ma n-2 vnutornych vrcholov, teda 2n-3 hran
• Pocet nezakorenenych fylogenetickych stromov s n listami:
  • a(3) = 1, a(4) = 3, a(n+1) = a(n) * (2n-3) a teda a(n) = 1 * 3 * 5 * ... * (2n-5) = (2n-5)!!
• Pocet zakorenenych fylogenetickych stromov s n listami:
  • zakoren strom s n listami kazdy 2n-3 sposobmi, teda (2n-3)!!

UCSC genome browser, cvičenia pre biológov

• Používajte čierne počítače, zvoľte v menu Bioinformatika, užívatel Bioinformatika
• V programe Firefox choďte na stránku UCSC genome browser http://genome.ucsc.edu
• Hore v modrom menu zvoľte Genomes
• Na ďalšej stránke zvoľte človeka a v menu Assembly zistite, kedy boli pridané posledné dve verzie ľudského genómu (hg18 a hg19)
• Na tej istej stránke dole nájdete stručný popis zvolenej verzie genómu. Pre ktoré chromozómy máme viacero alternatívnych verzií?
• Zadajte región chr21:31,200,000-31,400,000
• Zapnite si track Mapability na "pack" a track RepeatMasker prepnite na "full"
• Mapability: nakoľko sa daný úsek opakuje v genóme a či teda vieme jednoznačne jeho ready namapovať pri použití Next generation sequencing (detaily keď kliknete na linku "Mapability")
• Približne v strede zobrazeného regiónu je pokles mapovateľnosti. Akému typu opakovania zodpovedá? (pozrite track RepeatMasker)
• Zapnite si tracky "Assembly" a "Gaps" a pozrite si región chr2:110,000,000-110,300,000. Aká dlhá je neosekvenovaná medzera (gap) v strede tohto regiónu? Približnú veľkosť môžete odčítať z obrázku, presnejší údaj zistíte kliknutím na čierny obdĺžnik zodpovedajúci tejto medzere (úplne presnú dĺžku aj tak nepoznáme, nakoľko nie je osekvenovaná).
• Prejdite na genóm Rhesus, región chr7:59,022,000-59,024,000, zapnite si tracky Contigs, Gaps, Quality scores
• Aké typy problémov v kvalite sekvencie v tomto regióne vidíte?

Sekvenovanie genómov, cvičenia pre informatikov

Uvod do pravdepodobnosti, pocitanie pokrytia genomov

• Nas problem: spocitanie pokrytia
  • G = dlzka genomu, napr. 1 000 000
  • N = pocet segmentov (readov), napr. 10 000
  • L = dlzka readu, napr. 1000
  • Celkova dlzka segmentov NL, pokrytie (coverage) NL/G, v nasom pripade 10x
  • V priemere kazda baza pokryta 10x
  • Niektore su ale pokryte viackrat, ine menej.
  • Zaujimaju nas otazky typu: kolko baz ocakavame, ze bude pokrytych menej ako 3x?
  • Dolezite pri planovani experimentov (aku velke pokrytie potrebujem na dosiahnutie urcitej kvality)

• Uvod do pravdepodobnosti
  • Myslienkovy experiment, v ktorom vystupuje nahoda, napr. hod idealnou kockou/korunou
  • Vysledkom experimentu je nejaka hodnota (napr. cislo, alebo aj niekolko cisel, retazec)
  • Tuto neznamu hodnotu budeme volat nahodna premenna
  • Zaujima nas pravdepodobnost, s akou nahodna premenna nadobuda jednotlive mozne hodnoty
  • T.j. ak experiment opakujeme vela krat, ako casto uvidime nejaky vysledok
  • Priklad 1: hodime idealizovanou kockou, premenna X bude hodnota, ktoru dostaneme
  • Mozne hodnoty 1,2,..,6, kazda rovnako pravdepodobna
  • Piseme napr. Pr(X=2)=1/6
  • Priklad 2: hodime 2x kockou, nahodna premenna X bude sucet hodnot, ktore dostaneme
  • Mozne hodnoty: 2,3,...,12
  • Kazda dvojica hodnot na kocke rovnako pravdepodobna, t.j. pr. 1/36
  • Sucet 5 mozeme dostat 1+4,2+3,3+2,4+1 - t.j. P(X=5) = 4/36
  • Sucet 11 mozeme dostat 5+6 alebo 6+5, t.j. P(X=11) = 2/36
  • Rozdelenie pravdepodobnosti: 2: 1/36, 3:2/36, 4: 3/36, ... 7: 6/36, 8: 5/36 ... 12: 1/36
  • Sucet tychot hodnot je 1

• Pokrytie genomu: predpokladame, ze kazdy segment zacina na nahodnej pozicii zo vsetkych moznych G-L+1
• Takze ak premenna Y_i bude zaciatok i-teho segmentu, jej rozdelenie bude rovnomerne
  • P(Y_i=1) = P(Y_i=2) = ... = P(Y_i=G-L+1) = 1/G-L+1
• Uvazujme premennu X_j, ktora udava pocet segmentov pokryvajucich poziciu j
  • mozne hodnoty 0..N
  • i-ty segment pretina poziciu j s pravdepodobnostou L/(G-L+1), oznacme tuto hodnotu p
  • to iste ako keby sme N krat hodili mincou, na ktorej spadne hlava s pravd. p a znak 1-p a oznacili ako X_j pocet hlav
  • taketo rozdelenie pravdepodobnosti sa vola binomicke
  • P(X_j = k) = (N choose k) p^k (1-p)^(N-k)
  • Zle sa pocita pre velke N, preto sa niekedy pouziva aproximacia Poissonovym rozdelenim s parametrom lambda = Np, ktore ma $e^{{-\lambda }}\lambda ^{k}/k!$
  • Spat k sekvenovaniu: vieme spocitat rozdelenie pravdepodobnosti a tiez napr. P(X_i<3) = P(X_i=0)+P(X_i=1)+P(X_i=2) = 0.000045+0.00045+0.0023=0.0028 (v priemere ocakavame 45 baz nepokrytych, 2800 s menej ako 3 segmentami)
  • Takyto graf, odhad, vieme lahko spravit pre rozne pocty segmentov a tak naplanovat, kolko segmentov potrebujeme

• Chceme tiez odhadnut pocet kontigov (podla clanku Lander a Waterman 1988)
  • Ak niekolko baz vobec nie je pokrytych segmentami, prerusi sa kontig
  • Vieme, kolko baz je v priemere nepokrytych, ale niektore mozu byt vedla seba
  • Novy kontig vznikne aj ak sa susedne segmenty malo prekryvaju
  • Predpokladajme, ze na spojenie dvoch segmentov potrebujeme prekryv aspon T=50
  • Aka je pravdepodobnost, ze dany segment i bude posledny v kontigu?
  • Ziaden segment j!=i nesmie zacinat v prvych L-T bazach kontigu i
  • Kazdy segment tam zacina s pravdepodobnostou (L-T)/G, v priemere ich tam zacne N(L-T)/G
  • Pouzijeme Poissonovo rozdelenie pre $\lambda =N(L-T)/G$ a k=0, t.j. pravdepodobnost, ze tam nezacne ziaden je zhruba exp(-N(L-T)/G)
  • Pre N segmentov dostaneme priemerny pocet kontigov N*exp(-N(L-T)/G)
  • Ako keby sme dlzku segmentu skratili o dlzku prekryvu
  • Pre T=50 dostaneme priemerny pocet kontigov 0.75 (v skutocnosti vzdy aspon 1, my sme vsak neuvazovali vplyv konca sekvencie). Ak znizime N na 5000 (5x pokrytie) dostaneme 43 kontigov

• Tento jednoduchy model nepokryva vsetky faktory:
  • Segmenty nemaju rovnaku dlzku
  • Problemy v zostavovani kvoli chybam, opkaovaniam a pod.
  • Segmenty nie su rozlozene rovnomerne (cloning bias a pod.)
  • Vplyv koncov chromozomov
  • Uzitocny ako hruby odhad
  • Na spresnenie mozeme skusat spravit zlozitejsie modely, alebo simulovat data

• Poznamka: pravdepodobnosti z binomickeho rozdelenia mozeme lahko spocitat napr. statistickym softverom R. Tu su prikazy, ktore sa na to hodia, pre pripad, ze by vas to zaujimalo:

dbinom(10,1e4,0.001);  #(12.5% miest ma pokrytie presne 10)
pbinom(10,1e4,0.001,lower.tail=TRUE); #(58% miest ma pokrytie najviac 10)
dbinom(0:30,1e4,0.001); #tabulka pravdepodobnosti
 [1] 4.517335e-05 4.521856e-04 2.262965e-03 7.549258e-03 1.888637e-02
 [6] 3.779542e-02 6.302390e-02 9.007019e-02 1.126216e-01 1.251601e-01
[11] 1.251726e-01 1.137933e-01 9.481826e-02 7.292252e-02 5.207187e-02
[16] 3.470068e-02 2.167707e-02 1.274356e-02 7.074795e-03 3.720595e-03
[21] 1.858621e-03 8.841718e-04 4.014538e-04 1.743354e-04 7.254524e-05
[26] 2.897743e-05 1.112843e-05 4.115040e-06 1.467156e-06 5.050044e-07
[31] 1.680146e-07

Zhrnutie

• Pravdedpobnostny model: myslienkovy experiment, v ktorom vystupuje nahoda, napr. hod idealizovanou kockou
• Vysledok je hodnota, ktoru budeme volat nahodna premenna
• Tabulka, ktora pre kazdu moznu hodnotu nahodnej premennej urci je pravdepodobnost sa vola rozdelenie pravdepodobnosti, sucet hodnot v tabulke je 1
• Znacenie typu P(X=7)=0.1

• Priklad: mame genom dlzky G=1mil., nahodne umiestnime N=10000 segmentov dlzky L=1000
• Nahodna premenna X_i je pocet segmentov pokryvajucich urcitu poziciu i
• Podobne, ako keby sme N krat hodili kocku, ktora ma cca 1 promile sancu padnu ako hlava a 99.9% ako znak a pytame sa, kolko krat padne znak (1 promile sme dostali po zaukruhleni z L/(G-L+1))
• Rozdelenie pravdepobnosti sa v tomto pripade vola binomicke a existuje vzorec, ako ho spocitat
• Takyto model nam moze pomoct urcit, kolko segmentov potrebujeme osekvenovat, aby napr. aspon 95% pozicii bolo pokrytych aspon 4 segmentami

Zostavovanie genomu pomocou eulerovskych tahov

• Opakovanie z teorie grafov:
  • Hamiltonovska kruznica: cyklus, ktory prechadza kazdym vrcholom prave raz
  • Eulerov tah: tah, ktory prechadza po kazdej hrane prave raz
  • Zistit, ci ma graf H.c. je NP-tazke
  • Problem obchodneho cestujuceho: najst najlacnejsiu H.c. je tiez NP-tazky
  • Zistit ci ma graf E.t. a najst ho je lahke
    • neorientovany graf na E.t. <=> je suvisly a vsetky vrcholy okrem najviac dvoch maju parny stupen
    • orientovany graf ma E.t. z u do v <=> po pridani hrany z v do u je silne suvisly a vsetky vrcholy maju rovnako vchadzajucich ako vychadzajucich hran

• Najkratsie spolocne nadslovo: zjednodusena verzia problemu zostavovania genomov zo segmentov
• Mame danu mnozinu retazcov, chceme zostavit najkratsi retazec, ktoremu su vsetky podslova
• Priklad: GCCAAC,CCTGCC,ACCTTC zlozime v poradi 2,1,3: CCTGCCAACCTTC
• Mozeme si predstavit problem ako obmenu problemu obchodneho cestujuceho:
  • retazcom priradim vrcholy
  • dlzka hrany z a do b je $|b|-$ prekryv medzi a a b
  • specialny zaciatocny vrchol, z ktoreho hrany do kazdeho a s cenou |a|
  • specialny koncovy vrchol, do ktoreho hrana z kazdeho a s cenou 0
  • hrana z koncoveho do zaciatocneho vrcholu na uzavretie cyklu
  • Hamiltonovske kruznice zodpovedaju nadslovam
  • Ak najdeme najlacnejsiu Hamiltonovsku kruznicu, mame najkratsie nadslovo
• Ale vieme, ze problem obchodneho cestujuceho je NP-tazky
• Je toto dokaz, ze aj najkratsie spolocne nadslovo je NP-tazke?

• Pevzner, Tang and Waterman 2001 navrhuju namiesto Hamiltonovskej kruznice pouzit Eulerov tah (na inom grafe)
• deBruijnov graf stupna k:
  • vrcholy: podretazce dlzky k vsetkych vstupnych retazcov
  • hrany: nadvazujuce k-tice v ramci kazdeho segmentu (s prekryvom k-1)
• Priklad, k=2
• Chceme prejst po vsetkych hranach, chceme chodit co najmenej, takze po kazdej len raz - Eulerov tah
• V com je finta? Ako sme sa dostali od NP-tazkeho problemu k lahkemu?
• Co ak deBruijnov graf nema Eulerovsky tah? Znasobime niektore hrany, tieto budu zodpovedat opakovaniam
• Čo ak de Bruijnov graf má viacero Eulerovských ťahov?
  • Zoberieme taký ťah, ktorý obsahuje pôvodné segmenty ako podcesty
  • Zase ťažký problém, ale v praxi pomáhajú jednoduché pravidlá
• Opatrné riešenie: Ak z vrcholu 2 cesty, rozdeľ na kontigy
• Pouzitie sparovanych segmentov:
  • Nájdi vrcholy v grafe, ktoré im zodpovedajú
  • Ak je v grafe jediná cesta medzi týmito vrcholmi vhodnej dĺžky, premeň spárované segmenty na jeden veľký segment
• Dalsie problemy, ktore treba riesit
  • sekvenovacie chyby: vytvaraju "bubliny" alebo slepe cesty
  • dve vlakna
• Realne pouzivana technologia, aj pre sekvenovanie novej generacie, napr. Zerbino and Birney 2008 program Velvet

Sekvenovanie genómov, cvičenia pre biológov

• Cvicenia pre biológov zo sekvenovania, pracovna verzia poznamok
• Ospravedlnujem sa za chybajucu diakritiku
• Pozrite tiez grafy k pravdepodobnosti: http://compbio.fmph.uniba.sk/vyuka/mbi/poznamky/cb02.pdf

UCSC genome browser

• http://genome.ucsc.edu/
• Ukazali sme si niekolko trackov, ktore maju suvis so sekvenovanim a skladanim genomov, napr. gaps, quality score a pod. (v genomoch cloveka a macky)

Pravdepodobnost

• Nas problem: spocitanie pokrytia
  • G = dlzka genomu, napr. 1 000 000
  • N = pocet segmentov (readov), napr. 10 000
  • L = dlzka readu, napr. 1000
  • Celkova dlzka segmentov NL, pokrytie (coverage) NL/G, v nasom pripade 10x
  • V priemere kazda baza pokryta 10x
  • Niektore su ale pokryte viackrat, ine menej.
  • Zaujimaju nas otazky typu: kolko baz ocakavame, ze bude pokrytych menej ako 3x?
  • Dolezite pri planovani experimentov (aku velke pokrytie potrebujem na dosiahnutie urcitej kvality)

• Uvod do pravdepodobnosti
  • Myslienkovy experiment, v ktorom vystupuje nahoda, napr. hod idealnou kockou/korunou
  • Vysledkom experimentu je nejaka hodnota (napr. cislo, alebo aj niekolko cisel, retazec)
  • Tuto neznamu hodnotu budeme volat nahodna premenna
  • Zaujima nas pravdepodobnost, s akou nahodna premenna nadobuda jednotlive mozne hodnoty
  • T.j. ak experiment opakujeme vela krat, ako casto uvidime nejaky vysledok
  • Priklad 1: hodime idealizovanou kockou, premenna X bude hodnota, ktoru dostaneme
  • Mozne hodnoty 1,2,..,6, kazda rovnako pravdepodobna
  • Piseme napr. Pr(X=2)=1/6
  • Priklad 2: hodime 2x kockou, nahodna premenna X bude sucet hodnot, ktore dostaneme
  • Mozne hodnoty: 2,3,...,12
  • Kazda dvojica hodnot na kocke rovnako pravdepodobna, t.j. pr. 1/36
  • Sucet 5 mozeme dostat 1+4,2+3,3+2,4+1 - t.j. P(X=5) = 4/36
  • Sucet 11 mozeme dostat 5+6 alebo 6+5, t.j. P(X=11) = 2/36
  • Rozdelenie pravdepodobnosti: 2: 1/36, 3:2/36, 4: 3/36, ... 7: 6/36, 8: 5/36 ... 12: 1/36
  • Overte, ze sucet je 1

• Pokrytie genomu: predpokladame, ze kazdy segment zacina na nahodnej pozicii zo vsetkych moznych G-L+1
• Takze ak premenna Y_i bude zaciatok i-teho segmentu, jej rozdelenie bude rovnomerne
  • P(Y_i=1) = P(Y_i=2) = ... = P(Y_i=G-L+1) = 1/G-L+1
• Uvazujme premennu X_j, ktora udava pocet segmentov pokryvajucich poziciu j
  • mozne hodnoty 0..N
  • i-ty segment pretina poziciu j s pravdepodobnostou L/(G-L+1), oznacme tuto hodnotu p
  • to iste ako keby sme N krat hodili mincou, na ktorej spadne hlava s pravd. p a znak 1-p a oznacili ako X_j pocet hlav
  • Priklad: majme mincu, ktora ma hlavu s pr. 1/4 a hodime je 3x.

HHH 1/64
HHT 3/64
HTH 3/64
HTT 9/64
THH 3/64
THT 9/64
TTH 9/64
TTT 27/64

• P(X_j=3) = 1/64, P(X_j=2)=9/64, P(X_j=1)=27/64, P(X_j=0)=27/64
  • taketo rozdelenie pravdepodobnosti sa vola binomicke
  • P(X_j = k) = (N choose k) p^k (1-p)^(N-k), kde ${N \choose k}={\frac {N!}{k!(N-k)!}}$ a n! = 1*2*...*n
  • napr pre priklad s troma hodmi kockou P(X_j=2) = 3!/(2!*1!) * (1/4)^2 * (3/4)^1 = 9/64
  • Zle sa pocita pre velke N, preto sa niekedy pouziva aproximacia Poissonovym rozdelenim s parametrom lambda = Np, ktore ma $e^{{-\lambda }}\lambda ^{k}/k!$
  • Spat k sekvenovaniu: vieme spocitat rozdelenie pravdepodobnosti a tiez napr. P(X_i<3) = P(X_i=0)+P(X_i=1)+P(X_i=2) = 0.000045+0.00045+0.0023=0.0028 (v priemere ocakavame 45 baz nepokrytych)
  • Takyto graf, odhad, vieme lahko spravit pre rozne pocty segmentov a tak naplanovat, kolko segmentov potrebujeme

• Chceme tiez odhadnut pocet kontigov (nebrali sme na cviceni, uvedene len pre zaujimavost)
  • Ak niekolko baz vobec nie je pokrytych segmentami, prerusi sa kontig
  • Vieme, kolko baz je v priemere nepokrytych, ale niektore mozu byt vedla seba
  • Novy kontig vznikne aj ak sa susedne segmenty malo prekryvaju
  • Predpokladajme, ze na spojenie dvoch segmentov potrebujeme prekryv aspon T
  • Lander a Waterman 1988 odhadli, ze dany segment ma pravdepodobnost zhruba exp(-N(L-T)/G), ze bude posledny v kontigu
  • Pre N segmentov dostaneme priemerny pocet kontigov 1+N*exp(-N(L-T)/G)
  • Ako keby sme dlzku segmentu skratili o dlzku prekryvu
  • Pre T=50 dostaneme priemerny pocet kontigov 1.75, ak znizime N na 5000 (5x pokrytie) dostaneme 44 kontigov

• Tento jednoduchy model nepokryva vsetky faktory:
  • Segmenty nemaju rovnaku dlzku
  • Problemy v zostavovani kvoli chybam, opkaovaniam a pod.
  • Segmenty nie su rozlozene rovnomerne (cloning bias a pod.)
  • Vplyv koncov chromozomov
  • Uzitocny ako hruby odhad
  • Na spresnenie mozeme skusat spravit zlozitejsie modely, alebo simulovat data

• Poznamka: pravdepodobnosti z binomickeho rozdelenia mozeme lahko spocitat napr. statistickym softverom R. Tu su prikazy, ktore sa na to hodia, pre pripad, ze by vas to zaujimalo:

dbinom(10,1e4,0.001);  #(12.5% miest ma pokrytie presne 10)
pbinom(10,1e4,0.001,lower.tail=TRUE); #(58% miest ma pokrytie najviac 10)
dbinom(0:30,1e4,0.001); #tabulka pravdepodobnosti
 [1] 4.517335e-05 4.521856e-04 2.262965e-03 7.549258e-03 1.888637e-02
 [6] 3.779542e-02 6.302390e-02 9.007019e-02 1.126216e-01 1.251601e-01
[11] 1.251726e-01 1.137933e-01 9.481826e-02 7.292252e-02 5.207187e-02
[16] 3.470068e-02 2.167707e-02 1.274356e-02 7.074795e-03 3.720595e-03
[21] 1.858621e-03 8.841718e-04 4.014538e-04 1.743354e-04 7.254524e-05
[26] 2.897743e-05 1.112843e-05 4.115040e-06 1.467156e-06 5.050044e-07
[31] 1.680146e-07

Zhrnutie

• Pravdedpobnostny model: myslienkovy experiment, v ktorom vystupuje nahoda, napr. hod idealizovanou kockou
• Vysledok je hodnota, ktoru budeme volat nahodna premenna
• Tabulka, ktora pre kazdu moznu hodnotu nahodnej premennej urci je pravdepodobnost sa vola rozdelenie pravdepodobnosti, sucet hodnot v tabulke je 1
• Znacenie typu P(X=7)=0.1

• Priklad: mame genom dlzky G=1mil., nahodne umiestnime N=10000 segmentov dlzky L=1000
• Nahodna premenna X_i je pocet segmentov pokryvajucich urcitu poziciu i
• Podobne, ako keby sme N krat hodili kocku, ktora ma cca 1 promile sancu padnu ako hlava a 99.9% ako znak a pytame sa, kolko krat padne znak (1 promile sme dostali po zaukruhleni z L/(G-L+1))
• Rozdelenie pravdepobnosti sa v tomto pripade vola binomicke a existuje vzorec, ako ho spocitat
• Takyto model nam moze pomoct urcit, kolko segmentov potrebujeme osekvenovat, aby napr. aspon 95% pozicii bolo pokrytych aspon 4 segmentami

Eulerovský ťah

• Nestihli sme na cviceniach, uvedene pre zaujimavost
• Orientovany graf
• Silne suvisly orientovany graf
• Eulerovsky tah, uzatvoreny eulerovsky tah
• Pouzitie v de Bruijnovom grafe na zostavovanie genomu
• Silne súvislý orientovany graf ma uzatvoreny eulerovsky tah prave vtedy, ked v kazdom vrchole je rovnako vchadzajucich a vychadzajucich hran
  • ak ma tah, musi byt suvisly
  • ak ma tak, musia sediet pocty hran
  • ak je suvisly a sedia pocty hran:
    • prechadzka po grafe, kym sa da pokracovat nepouzitou hranou
    • musime skoncit, kde sme zacali (z ostatnych vrcholov vieme pokracovat)
    • zvysok grafu stale splna podmienku s poctami hran
    • najdi dalsi vrchol, ktory je na tahu a ma este nepouzite hrany
    • musi taky existovat
    • pokracuj v tom istom, az kym sa nevratime, kde sme zacali
    • cele opakujeme, az kym sa neminu hrany
• Poznamka: v neorientovanych grafoch musia byt stupne vsetkych vrcholov parne
• Ak hladame otvoreny tah, spojime zaciatocny a koncovy vrchol hranou a najdeme uzatvoreny tah

Zarovnávanie sekvencií, cvičenia pre biológov

Opakovanie dynamickeho programovania

• Uvazujme skorovanie zhoda +2, nezhoda -1, medzera -1
• Retazce TAACGG a CACACT

Globalne zarovnanie

• Rekurencia: A[i,j] = max {A[i-1,j]-1, A[i,j-1]-1, A[i-1,j-1]+c(p_i, q_j) }, pricom A[0,i]=-i, A[i,0]=-i

        C   A   C   A   C   T
    0  -1  -2  -3  -4  -5  -6
T  -1  -1  -2  -3  -4  -5  -3
A  -2  -2   1   0  -1  -2  -3
A  -3  -3   0   0   2   1   0
C  -4  -1  -1   2   1   4   3
G  -5  -2  -2   1   1   3   3
G  -6  -3  -3   0   0   2   2

TA-ACGG
CACAC-T

Lokalne zarovnanie

• Rekurencia: A[i,j] = max {0, A[i-1,j]-1, A[i,j-1]-1, A[i-1,j-1]+c(p_i, q_j) }, pricom A[0,i]=0, A[i,0]=0

        C   A   C   A   C   T
    0   0   0   0   0   0   0
T   0   0   0   0   0   0   2
A   0   0   2   1   2   1   1
A   0   0   2   1   3   2   1
C   0   2   1   4   3   5   4
G   0   1   1   3   3   4   4
G   0   0   0   2   2   3   3

ACAC
A-AC

Dotploty

• Dotplot je graf, ktory ma na kazdej osi jednu sekvenciu a ciarky zobrazuju lokalne zarovnania (cesty v matici)
• Na slidoch mame niekolko prikladov dotplotov porovnavajucich rozne mitochondrialne genomy
• Tieto boli vytvorene pomocou nastroja YASS http://bioinfo.lifl.fr/yass/yass.php
• Dalsi priklad je zarovnanie genu Oaz Drosophila zinc finger s genomickym usekom chr2R:10,346,241-10,352,965
• Trochu iny dotplot, ktory funguje pre proteiny a nerobi lokalne zarovnania, iba spocita skore bez medzier v kazdom okne danej vysky a nakresli ciaru ak pre kroci urcenu hodnotu
• http://emboss.bioinformatics.nl/cgi-bin/emboss/dotmatcher
• Vyskusame protein escargot voci sebe s hodnotami http://pfam.sanger.ac.uk/protein/ESCA_DROME window 8 threshold 24
• Pomocou YASSu vyskusame kluster zhlukov PRAME z ludskeho genomu
• Obrázky dotplotov na slidoch

Prakticke cvicenie

• Na stranke UCSC genome browseru http://genome.ucsc.edu/ si zadajte ludsky genom, najnovsiu verziu (hg19)
• V druhom okne/tabe si otvorte YASS server na adrese http://bioinfo.lifl.fr/yass/yass.php
• V genome browseri zadajte region chr21:10,018,860-10,019,178
  • tento región obsahuje Alu repeat, ktorý tvorí cca 10% ľudského genómu, viac ako milión kópií
  • zobrazte si DNA sekvenciu tohto useku takto: na hornej modrej liste zvolte DNA, na dalsej obrazovke tlacidlo get DNA
• DNA sekvenciu Alu opakovania chceme zarovnat samu k sebe programom YASS
  • DNA sekvenciu Alu opakovania skopirujte do okienka "Paste your sequences" v stranke Yass-u a dvakrat stlacte tlacidlo Select vedla okienka
  • Nizsie v casti "Selected DNA sequence(s)" by sa Vam malo v oboch riadkoch objavit "Pasted file 1"
  • Nizsie v casti "Parameters" zvolte "E-value threshold" 0.01 a stlacte "Run YASS"
  • Vo vysledkoch si pozrite Dotplot, co z neho viete usudit o podobnosti jednotlivych casti Alu opakovania?
  • Vo vysledkoch si pozrite Raw: blast, ake su suradnice opakujucej sa casti a kolko zarovnanie obsahuje zhod/nezhod/medzier? (Pozor, prve zarovnanie je cela sekvencia sama k sebe, druhe je asi to, co chcete)
• V genome browseri zadajte poziciu chr21:9,440,000-9,450,000 (10kb sekvencie na chromozome 21, s niekolkymi vyskytmi Alu)
• Chceme teraz porovnat tento usek genomu so sekvenciou Alu pomocou YASSu
  • Ako predtym si stiahnite DNA sekvenciu tohto useku
  • V YASSe chodte sipkou spat na formular
  • Skopirujte DNA sekvenciu do YASSoveho formulara, do okienka vpravo (vyznacit si ju mozete klavesovou kombinaciou Ctrl-A alebo Select All v menu Edit),
  • V casti formulara Selected DNA sequence(s) stlacte Remove pri hornom riadku
  • Pri pravom okienku, kam ste nakopirovali sekvenciu, stlacte Select
  • Zase stlacte Run YASS
  • Pozrite si vysledok ako Dotplot, kolko opakovani Alu ste nasli? Preco je jedno cervene?
  • Pozrite si Raw: blast, na kolko percent sa podoba najpodobnejsia a na kolko druha najpodobnejsia kopia?

Zarovnávanie sekvencií 2, cvičenia pre biológov

Uvod do pravdepodobnosti

Opakovanie E-value:

• Priklady k tejto casti v prezentacii [1]
• Mame dotaz dlzky m, databazu dlzky n, skore najlepsieho lok. zarovnania S
• E-value je ocakavany pocet zarovnani so skore aspon S ak dotaz aj databaza su nahodne
• Hrackarsky priklad: dotaz dlzky m=4, databaza dlzky n=200, S=4, t.j. presny vyskyt dotazu v databaze (pre presne vyskyty sa pravdepodobnosti pocitaju ovela lahsie ako ked dovolime nezhody a medzery)
• Pravdepodobnostny model: Myslienkovy experiment obsahujuci nahodne javy (hod kockou a pod.)
• Priklad: pravdepodobnostny model na generovanie nahodnej sekvencie:
• Mame vrece s gulockami oznacenyni A,C,G,T, pricom gulocok oznacenych A je 30%, C 20%, G 20% a T 30%.
• Vytiahneme gulicku, zapiseme si pismeno, hodime ju naspat, zamiesame a opakujeme s dalsim pismenom atd az kym nevygenerujeme m pismen pre dotaz a n pismen pre databazu
  • V nasom modeli mame obsah GC 40%
• Pre nase vygenerovane sekvencie spocitame, kolkokrat sa dotaz vyskytuje v databaze
• Cely experiment opakujeme vela krat a spocitame priemerny pocet vyskytov, co bude odhad E-value
• Vieme spocitat aj presne, ale na to potrebujeme trochu pravdepodobnosti

Zaklady pravdepodobnosti

• Vytiahnime z mechu 2x gulicku. Prve pismeno ktore nam vyjde oznacme X1 a druhe X2
• Pr(X1=A) = 0.3, Pr(X2=C)=0.2
• Pr(X1=A a X2=C) = Pr(X1=A)*Pr(X2=C) = 0.3*0.2 = 0.06
  • T.j. sanca ze dostaneme sekvenciu AC po dvoch tahaniach je 6%
  • Ak ratame pravdepodobnost, ze sa dve nezavisle udalosti stanu, ich pravdepodobnosti nasobime. V tomto pripade to ci X1=A je nezavisle od toho, ci X2=C
• Pr(X1 je A alebo G) = Pr(X1=A)+Pr(X1=G) = 0.3+0.2 = 0.5
  • Pravdepodobnost ze prve pismeno bude A alebo G je 50%
  • Pravdepodobnosti navzajom sa vylucujucich udalosti (X1=A a X1=G) sa mozu scitat, cim dostaneme pravdepodobnost, ze aspon jedna z nich nastane
• Pr(v sekvencii je aspon jedno A) = Pr(X1=A alebo X2=A) nemozeme pocitat ako Pr(X1=A)+Pr(X2=A), lebo sa navzajom nevylucuju a pripad ze X1=X2=A by sme zapocitali dvakrat
• Spravne je Pr(v sekvencii je aspon jedno A) = Pr(X1=A) + Pr(X1 <> A a X2=A) = 0.3+0.7*0.3 = 0.51
• Pr(X1=X2) = Pr(X1=X2=A) + Pr(X1=X2=C) + Pr(X1=X2=G) + Pr(X1=X2=T) = 0.3*0.3+0.2*0.2+0.2*0.2+0.3*0.3 = 0.26.

Naspat k E-value

• Mame dotaz X1X2...Xm a databazu Y1Y2...Yn. Aka je sanca, ze X sa nachadza na zaciatku Y, t.j. ze X1=Y1, X2=2... Xm=Ym?
• Nezavisle udalosti, t.j. Pr(X1=Y1)*Pr(X2=T2)*...*Pr(Xm=Ym) = 0.26 * 0.26 * ... * 0.26 = 0.26^m (0.26 na m-tu)
• Pravdepodobnost, ze sa X nachadza hocikde v Y je tazsie zratat - moze sa vyskytovat aj viackrat, udalosti za sa nachadza na pozicii i sa navzajom nevylucuju
• Ale priemerny pocet vyskytov na pozicii i je 0.26^m, mame n-m+1 pozicii zaciatku, takze celkovy ocakavany pocet vyskytov je (n-m+1)0.26^m
• Ak je n velke v porovnani s m, clen -m+1 mozeme zanedbat, t.j. mame zhruba n*0.26^m
• Ako sa meni toto cislo s velkostou databazy n? linearne rastie (zdvojnasobime databazu, zdvojnasobi sa e-value)
• Ako sa meni s m? Exponencialne klesa, t.j. ak predlzime dotaz o 1, E-value klesne zhruba na stvrtinu
• Ako sa meni s GC obsahom? Zakodovane v cisle 0.26 - vyskusame si pri pocitaci

BLAT vs BLAST, nets v UCSC browseri

• Program BLAT v UCSC browseri rychlo vyhladava sekvencie v genome, ale nevie najst slabsie podobnosti
  • Vhodne pouzitie: zarovnanie EST ku genomu, presne urcenie suradnic nejakej sekvencie, a pod.
• NCBI BLAST je pomalsi, ale najde viacej podobnosti
• Aj BLAST vsak toho moze vela vynechat
• Net tracky v UCSC genome browseri nam umoznuju prechadzat medzi homologickymi oblastami roznych genomov

Prakticke cvicenie pri pocitaci

• Na stránke http://compbio.fmph.uniba.sk/vyuka/mbi/poznamky/cb04 nájdete súbory k tomuto cvičeniu.

E-values

• Stiahnite si z tejto stranky subor evalue.ods, ulozte si ho na plochu, potom ho otvorte v OpenOffice
• Doplnte chybajuci vzorec v stlpci D, potom rozkopirujte vzorce v stlpcoch B-F na celu dlzku tabulky
• Aky rozsah E-value dostavate v poslednom stlpci? Pre ktore obsahy GC je najvyssi a najnizsi? Ako sa zmenia, ked zvysite hodnotu m o 1? Ak chcete, mozete si nakreslit aj graf.

BLAT/BLAST

• Na stránke so súbormi otvorte subor rtpcr-seq.txt. Obsahuje sekvenciu, ktorá vznikla pomocou RT-PCR na ľudských cDNA knižniciach
• Choďte na UCSC genome browser (http//genome.ucsc.edu/), na modrej lište zvoľte BLAT, zadajte túto sekvenciu a hľadajte ju v ľudskom genóme. Akú podobnosť (IDENTITY) má najsilnejší nájdený výskyt? Aký dlhý úsek genómu zasahuje? (SPAN). Všimnite si, že ostatné výskyty sú oveľa kratšie.
• V stĺpci ACTIONS si pomocou Details môžete pozrieť detaily zarovnania a pomocou Browser si pozrieť príslušný úsek genómu.
• V tomto úseku genómu si zapnite track Vertebrate net na full a kliknutím na farebnú čiaru na obrázku pre tento track zistite, na ktorom chromozóme kuraťa sa vyskytuje homologický úsek.
• Skusme tu istu sekvenciu namapovat do genomu sliepky: stlacte najprv na hornej modrej liste Genomes, zvolte Vertebrates a Chicken a potom na hornej liste BLAT. Do okienka zadajte tu istu sekvenciu. Akú podobnosť a dĺžku má najsilnejší nájdený výskyt teraz? Na ktorom je chromozóme?
• Skúsme to isté v NCBI blaste: Choďte na http://blast.ncbi.nlm.nih.gov/ zvoľte nucleotide blast, database others a z menu reference genomic sequence, organism chicken (taxid:9031), program blastn
• Aka je dlzka, identity a E-value najlepsieho zarovnania? Na ktorom je chromozome?

RT PCR sekvencia z cvičenia vyššie

AACCATGGGTATATACGACTCACTATAGGGGGATATCAGCTGGGATGGCAAATAATGATTTTATTTTGAC
TGATAGTGACCTGTTCGTTGCAACAAATTGATAAGCAATGCTTTCTTATAATGCCAACTTTGTACAAGAA
AGTTGGGCAGGTGTGTTTTTTGTCCTTCAGGTAGCCGAAGAGCATCTCCAGGCCCCCCTCCACCAGCTCC
GGCAGAGGCTTGGATAAAGGGTTGTGGGAAATGTGGAGCCCTTTGTCCATGGGATTCCAGGCGATCCTCA
CCAGTCTACACAGCAGGTGGAGTTCGCTCGGGAGGGTCTGGATGTCATTGTTGTTGAGGTTCAGCAGCTC
CAGGCTGGTGACCAGGCAAAGCGACCTCGGGAAGGAGTGGATGTTGTTGCCCTCTGCGATGAAGATCTGC
AGGCTGGCCAGGTGCTGGATGCTCTCAGCGATGTTTTCCAGGCGATTCGAGCCCACGTGCAAGAAAATCA
GTTCCTTCAGGGAGAACACACACATGGGGATGTGCGCGAAGAAGTTGTTGCTGAGGTTTAGCTTCCTCAG
TCTAGAGAGGTCGGCGAAGCATGCAGGGAGCTGGGACAGGCAGTTGTGCGACAAGCTCAGGACCTCCAGC
TTTCGGCACAAGCTCAGCTCGGCCGGCACCTCTGTCAGGCAGTTCATGTTGACAAACAGGACCTTGAGGC
ACTGTAGGAGGCTCACTTCTCTGGGCAGGCTCTTCAGGCGGTTCCCGCACAAGTTCAGGACCACGATCCG
GGTCAGTTTCCCCACCTCGGGGAGGGAGAACCCCGGAGCTGGTTGTGAGACAAATTGAGTTTCTGGACCC
CCGAAAAGCCCCCACAAAAAGCCG

Zarovnávanie sekvencií, cvičenia pre informatikov

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\}$

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 roznch pribuznych problemov bez toho, aby sme vzdy vymyslali novu rekurenciu.

Lokalne 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 hrnay do t (ako?)

Afinne skore medzier

• Napr. otvorenie medzery -2, pokracovanie medzery -1

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

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)
• DOKONCIT

Linearna pamat: Hirshbergov algoritmus 1975

• DOKONCIT

Vypisanie vsetkych najlepsich rieseni

• 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

Zarovnávanie sekvencií 2, cvičenia pre informatikov

• Cvicenia pre informatikov z jadier s medzerami (spaced seeds)
• Prva cast na slidoch

Vzorec na vypocet senzitivity jadra bez medzier

• Uvazujme jadro bez medzier dlzky w (ako v programe BLAST pre nukleotidy)
• Uvazujme pravdepodobnostny model zarovnania, v ktorom ma kazda pozicia pravdepodobnost p, ze bude zhoda a (1-p), ze bude nezhoda, medzery neuvazujeme.
• Nech f(L) je pravdepodobnost vyskytu jadra v zarovnani dlzky L, t.j. pravdepodobnost w za sebou iducich zhod.
• $f(L)=\left\{{\begin{array}{ll}0&{\mbox{ak }}L<w\\p^{w}+\sum _{{i=0}}^{{w-1}}p^{i}(1-p)f(L-i-1)&{\mbox{ak }}L\geq w\\\end{array}}\right.$

Ine typy jadier

• Vektorove jadra
• zarovnanie ako postupnost skore (napr 0/1 pre zhoda/nezhoda, male cele cisla pre BLOSUM63 maticu skorovania aminokyselin)
• Jadro je vektor vah v_q\dots v_k dlzky k a prahova hodnota T
• Okno zarovnania s_1\dots d_k je vyskyt jadra ak $\sum _{{i=1}}^{k}s_{i}v_{i}\geq T$
• Vyjadrite jadra s medzerami, BLAT-ove a BLASTp jadra ako vektorove jadra
• Ako mozeme rychlo hladat jadra roznych typov v sekvenciach?
  • Jadro bez medzier (BLAST)
    • trivialna hash tabulka velkost 4^w
    • Aho-Corasickovej algoritmus pre vsetky w-tice v prvom retazci
    • sufixovy strom
    • atd
  • Jadro s medzerami
    • Pouzijeme trivialnu hash tabulku, ale ukladame len znaky na neignorovanych poziciach
  • Jadro BLAT-u
    • Hash tabulka pre jednu sekvenciu, pre druhu hladame vsetky retazce v okoli danej k-tice ktore sa lisia najviac na jednom mieste (je ich 3k+1)
    • Ak by sme dovolili viac rozdielov, ich pocet by rychlo narastal

Senzitivita jadra s medzerami

• Da sa rekurentny vzorec vyssie rozsirit na jadro s medzerami?
• $f(L)=\left\{{\begin{array}{ll}0&{\mbox{ak }}L<w\\p^{w}+\sum _{{i=0}}^{{w-1}}p^{i}(1-p)f(?)&{\mbox{ak }}L\geq w\\\end{array}}\right.$
• Co by sme pouzili namiesto otaznika?
• Napr pre jadro 1011 a zaciatok zarovnania 1?10 moze byt vyskyt na pozicii 3, ale jeho pravdepodobnost nie je f(L-2), lebo mame dane prve dva znaky

Vseobecny algoritmus pre jadra s medzerami/vektorove jadra

• Zial vo vseobecnosti exponencialny od dlzky jadra, ale pre male jadra sa da zbehnut v rozumnom case
• Zarovnanie si predstavme ako nejaku postupnost nad abecedou D (napr D={0,1} kde 0 je zhoda, 1 nezhoda). Nech q_d je pravdepodobnost skore d v zarovnani
• Uvazujme jadro dlzky k, co si vieme predstavit ako nejake pravidlo, ktore urci, ci dany retazec dlzky k nad D je alebo nie je vyskytom jadra
• H nech je mnozina vsetkych vyskytov jadra (t.j. mnozina retazcov x dlzky k nad mnozinou D, ktore splnaju podmienku jadra)
• Nech H_p je mnozina prefixov retazcov z H (mozne vyskyty) a H_g je mnozina prefixov retazcov z H, ktorych kazde rozsirenie na dlzku k je vyskyt (garantovane vyskyty)
• $H\subseteq H_{g}\subseteq H_{g}$
• Nech f(L,x) je pravdepodobnost vyskytu jadra v zarovnani dlzky L, ktore zacina na retazec x.
• $f(L,x)=\left\{{\begin{array}{ll}0&{\mbox{ak }}L<k\\1&{\mbox{ak }}L\geq k\wedge x\in H_{g}\\\sum _{{d\in D}}q_{d}f(L,xd)&{\mbox{ak }}L\geq k\wedge x\in H_{p}\setminus H_{g}\\f(L-|y|,z)&{\mbox{ak }}L\geq k\wedge x\notin H_{p}\wedge x=yz,z\in H_{P},z{\mbox{ najdlhsie mozne}}\\\end{array}}\right.$
• f(L,x) budeme pocitat pre rastuce L a pre kazde L od najdlhsich x. Staci uvazovat x z H_p alebo x, ktore dostaneme pridanim jedneho znaku k slovu z H_p.
• vysledok je $f(L,\varepsilon )$
• Pre jadro s g medzerami $|H|=2^{g}$ a $|H_{p}|\leq k2^{g}$.

Hľadanie génov a anotácia sekvencií, cvičenia pre biológov

• Pozri slidy [2]

Príklady stavových automatov pre HMM

Uvazujme HMM so specialnym zaciatocnym stavom b a koncovym stavom e, ktore nic negeneruju.

• Nakreslite HMM (stavovy diagma), ktory generuje sekvencie, ktore zacinaju niekolkymi cervenymi pismenami a potom obsahuju niekolko modrych
• Ako treba zmenit HMM, aby dovoloval ako "niekolko" aj nula?
• Ako treba zmenit HMM, aby pocet cervenych aj modrych bol vzdy parne cislo?
• Ako zmenit HMM, aby sa striedali cervene a modre kusy parnej dlzky?

V dalsich prikladoch uvazujeme aj to, ktore pismena su v ktorom stave povolene (pravdepodobnost emisie > 0) a ktore su zakazane

• cervena sekvencia dlzky dva, ktora zacina na A
• cervena sekvencia dlzky dva, ktora je hocico okrem AA

Dalsi biologicky priklad HMM: topologia transmembranovych proteinov.

Histónové modifikácie

• A. Barski, S. Cuddapah, K. Cui, T. Roh, D. Schones, Z. Wang, G. Wei, I. Chepelev, K. Zhao (2007) High-Resolution Profiling of Histone Methylations in the Human Genome Cell, Volume 129, Issue 4, Pages 823-837 pdf

Gény v ľudskom genóme

• What is a gene, post-ENCODE? History and updated definition. Gerstein MB, Bruce C, Rozowsky JS, Zheng D, Du J, Korbel JO, Emanuelsson O, Zhang ZD, Weissman S, Snyder M.
• Most "dark matter" transcripts are associated with known genes. H Van Bakel, C Nislow, BJ Blencowe, TR Hughes - PLoS Biol, 2010
• Transcribed dark matter: meaning or myth? CP Ponting, TG Belgard - Human molecular genetics, 2010

• Priklady gene finderov v UCSC browseri, ludsky genom
  • chr2:149,242,117-149,247,381 kusok genu
  • chr2:149,118,916-149,120,420 zaciatok genu
  • chr2:149,115,905-149,266,404 cely gen
  • chr21:33,860,819-34,247,408 priklad z encode

Hľadanie génov v prokaryotických genómoch

• ORF: open reading frame, jednoduche hladanie
• ako najst zaciatok, ako rozlisit psedogeny a nahodne ORF-y
• samotrenujuce sa HMM, codon bias, GC%

E. coli http://nar.oxfordjournals.org/content/34/1/1.full

• Prvykrat sekvenovana a anotovana 1997
• Porovnanie s verziou 2005 (oprava sekvenovacích chýb aj chýb v anotácii)
  • 682 zmien v start kodone
  • 31 génov zrušených
  • 48 nových génov
  • Celkovo asi 4464 génov

Programy na anotovanie prokaryotických genómov

• ORFfinder at NCBI [3]
• Glimmer at NCBI [4]
• GeneMark at NCBI [5]
• tRNAscan-SE [6]

Skryté Markovove modely, cvičenia pre informatikov

Opakovanie

Parametre HMM:

• $a_{{u,v}}$: prechodová pravdepodobnosť zo stavu $u$ do stavu $v$
• $e_{{u,x}}$: pravdepodobnosť emisie $x$ v stave $u$
• $\pi _{{u}}$: pravdepodobnosť, že začneme v stave $u$

• Sekvencia $S=S_{1}S_{2}\dots S_{n}$
• Anotácia $A=A_{1}A_{2}\dots A_{n}$

$Pr(S,A)=\pi _{{A_{1}}}e_{{A_{1},S_{1}}}\prod _{{i=2}}^{n}a_{{A_{{i-1,A_{i}}}}}e_{{A_{i},S_{i}}}$

Trénovanie

Proces, pri ktorom sa snažíme čo najlepšie odhadnúť pravdepodobnosti $a_{{u,v}}$ a $e_{{u,x}}$ v modeli podľa trénovacích dát

Usudzovanie (inferencia)

Proces, pri ktorom sa snažíme pre sekvenciu $S$ nájsť anotáciu $A$, ktorá sekvenciu $S$ emituje s veľkou pravdepodobnosťou.

Inferencia pomocou najpravdepodobnejšej cesty, Viterbiho algoritmus

Hľadáme najpravdepodobnejšiu postupnosť stavov $A$, teda $\arg \max _{A}\Pr(A,S)$. Úlohu budeme riešiť dynamickým programovaním.

• Podproblém $V[i,u]$: Pravdepodobnosť najpravdepodobnejšej cesty končiacej po $i$ krokoch v stave $u$, pričom vygeneruje $S_{1}S_{2}\dots S_{i}$.

• Rekurencia:
  • $V[1,u]=\pi _{u}e_{{u,s_{1}}}$ (*)
  • $V[i,u]=\max _{w}V[i-1,w]a_{{w,u}}e_{{u,x_{i}}}$ (**)

Algoritmus:

1. Nainicializuj $V[1,*]$ podľa (*)
2. for i=2 to n=dĺžka reťazca

for u=1 to m=počet stavov

vypočítaj $V[i,u]$ pomocou (**)

1. Maximálne $V[n,j]$ je pravdepodobnosť najpravdepodobnejšej cesty

Aby sme vypísali anotáciu, pamätáme si pre každé $V[i,u]$ stav $w$, ktorý viedol k maximálnej hodnote vo vzorci (**).

Zložitosť: $O(nm^{2})$

Poznámka: pre dlhé sekvencie budú čísla $V[i,u]$ veľmi malé a môže dôjsť k podtečeniu. V praxi teda používame zlogarimované hodnoty, namiesto násobenia súčet.

Inferencia - dopredný algoritmus

Aká je celková pravdepodobnosť, že vygenerujeme sekvenciu $S$, t.j. $\sum _{A}Pr(A,S).$ Podobný algoritmus ako Viterbiho.

Podproblém $F_{{i,u}}$: pravdepodobnosť, že po $i$ krokoch vygenerujeme $s_{1},s_{2},\dots s_{i}$ a dostaneme sa do stavu $u$.

$F_{{i,u}}=\Pr(A_{i}=u\wedge S_{1},S_{2},\dots ,S_{i})=\sum _{{A_{1},A_{2},\dots ,A_{i}=j}}\Pr(A_{1},A_{2},...,A_{i}\wedge S_{1},S_{2},...,S_{i})$

$F_{{1,u}}=\pi _{u}e_{{u,x_{1}}}$

$F_{{i,u}}=\sum _{v}F_{{i-1,v}}a_{{v,u}}e_{{u,x_{i}}}$

Celková pravdepodobnosť $\sum _{u}F_{{n,u}}$

Inferencia - posterior decoding

Aposteriórna pravdepodobnosť stavu u na pozícii i: $Pr(A_{i}=u|S_{1}\dots S_{n})$

Pre každý index i chceme nájsť stav u s najväčšiou aposteriórnou pravdepodobnosťou, dostaneme tak inú možnú anotáciu.

Spustíme dopredný algoritmus a jeho symetrickú verziu, spätný algoritmus, ktorý počíta hodnoty $B[i,u]=\Pr(A_{i}=u\wedge S_{{i+1}}\dots S_{n})$

Aposteriórna pravdepodobnosť stavu u na pozícii i: $Pr(A_{i}=u|S_{1}\dots S_{n})=F[i,u]B[i,u]/\sum _{u}F_{{n,u}}.$

Posterior decoding uvažuje všetky anotácie, nielen jednu s najvyššou pravdepodobnosťou. Môže však vypísať anotáciu, ktorá má sama o sebe nulovú pravdepodobnosť (napr. počet kódujúcich báz v géne nie je deliteľný 3).

Trénovanie HMM

• Stavový priestor + povolené prechody väčšinou ručne
• Parametre (pravdepodobnosti prechodu, emisie a počiatočné) automaticky z trénovacích sekvencií
  • Ak máme anotované trénovacie sekvencie, jednoducho počítame frekvencie
  • Ak máme iba neanotované sekvencie, snažíme sa maximalizovať vierohodnosť trénovacích dát v modeli. Používajú sa heuristické iteratívne algoritmy, napr. Baum-Welchov, ktorý je verziou všeobecnejšieho algoritmu EM (expectation maximization).
• Čím zložitejší model a viac parametrov máme, tým potrebujeme viac trénovacích dát, aby nedošlo k preučeniu, t.j. k situácii, keď model dobre zodpvedá nejakým zvláštnostiam trénovacích dát, nie však ďalším dátam.
• Presnosť modelu testujeme na zvláštnych testovacích dátach, ktoré sme nepoužili na trénovanie.

Tvorba stavového priestoru modelu

• Promótor + niekoľko eukaryotických génov
• Repeaty v intrónoch: multiple path problem
• Intrón má dĺžku aspoň 10

Zovšeobecnené HMM

• Problém s modelovaním rozdelenia dĺžok - v základnom modeli je geometrické
• Zovseobecnene HMM v jednom stave vygenerujú viac znakov
• Viterbi alg. a spol budú pomalšie (kvadratické od dĺžky sekvencie)

Evolúcia, cvičenia pre biológov

• V tomto cviceni potrebujete pristup k niekolkym suborom a linkam, takze si otvorte tento text na stranke predmetu http://compbio.fmph.uniba.sk/vyuka/mbi/ (Prednasky a poznamky, v tabulke tento tyzden, CB: Poznamky)

Geny, evolucia a komparativna genomika v UCSC genome browseri

• Chodte na stranku http://genome.ucsc.edu/
• Zvolte starsiu verziu ludskeho genomu hg18, ktora ma viac informacii
• Do okienka position zadajte gen MAGEA2B a potom zvolte vyskyt na pozicii cca chrX:151636040-151637735 (ma dva vyskyty)
  • Dostanete sa tam aj touto linkou: [7]
• Ak date 3x zoom out, mozete si vsimnut, ze tento gen ma viacero foriem zostrihu, ktore sa ale lisia iba v 5' UTR
• Vsimnite si ENCODE Regulation Super-track, ktory zobrazuje ChIP data pre niektore histonove modifikacie
• Vela veci sa mozete dozvediet klikanim na rozne casti broswera: napr, kliknutim na gen si mozete precitat o jeho funkcii, kliknutim na listu ku tracku (lavy okraj obazku) sa dozviete viac o tracku a mozete nastavovat parametre zobrazenia

• V casti Genes and Gene Prediction Tracks zapnite track Pos Sel Genes, ktory obsahuje geny s pozitivnym vyberom (cervenou, pripadne slabsie fialovou a modrou)
• Ked kliknete na cerveny obdlznik pre tento gen, uvidite, v ktorych castiach fylogenetickeho stromu bol detekovany pozitivny vyber

• V casti multiz alignments vidite zarovnania k roznym inym genomom (da sa zapinat, ze ku ktorym). Mozete si pozriet, ako sa uroven zarovnania zmeni ked sa priblizujeme a vzdalujeme (zoom in/zoom out).
• Ked sa priblizite spat na gen MAGEA2B a potom tak, aby ste boli na urovni "base", t.j. zobrazenych cca 100bp, v obdlzniku multiz alignment uvidite zarovnanie s homologickym usekom v inych genomoch. Konkretne v MAGEA2B vidime pomerne dost rozdielov v proteine medzi clovekom a makakom rezus, vdaka ktorym bol zrejme klasifikovany ako pod pozitivnym vyberom.
• V casti conservation by PhyloP vidime graf toho, ako silne su zachovane jednotlive stlpce zarovnania
• Da sa zapnut aj track 28-Way Most Cons, ktory zobrazuje konkretne useky, ktore su najvac konzervovane
  • Ak chceme zistit, kolko percent genomu tieto useky pokryvaju, ideme na modrej liste do casti Tables, zvolime group Comparative genomics, track 28-way most cons, table Mammal alebo Vertebrate, region zvolime genome (v celom genome) a stlacime tlacidlo Summary/statistics

Objavenie genu HAR1 pomocou komparativnej genomiky

• Pollard KS, Salama SR, Lambert N, et al. (September 2006). "An RNA gene expressed during cortical development evolved rapidly in humans". Nature 443 (7108): 167–72. doi:10.1038/nature05113. PMID 16915236. pdf
• Zobrali vsetky regiony dlzky aspon 100 s > 96% podobnostou medzi simpanzom a mysou/potkanom (35,000)
• Porovnali s ostatnymi cicavcami, zistili, ktore maju vela mutaci v cloveku ale malo inde (pravdepodobnostny model)
• 49 statisticky signifikantych regionov, 96% nekoduje oblasti
• Najvyznamnejsi HAR1: 118nt, 18 substitucii u cloveka, ocakavali by sme 0.27. Iba 2 zmeny medzi simpanzom a sliepkou (310 milionov rokov), ale nebol najdeny v rybach a zabe.
• Nezda sa byt polymorfny u cloveka
• Prekryvajuce sa RNA geny HAR1R a HAR1F
• HAR1F je exprimovany v neokortexe u 7 a 9 tyzdennych embrii, neskor aj v inych castiach mozgu (u cloveka aj inych primatov)
• Vsetky substitucie v cloveku A/T->C/G, stabilnejsia RNA struktura (ale tiez su blizko k telomere, kde je viacej takychto mutacii kvoli rekombinacii a biased gene conversion)
• Mozete si pozriet tento region v browseri: chr20:61,203,911-61,204,071, pricom ak sa este priblizite, uvidite zarovnanie aj s bazami a mozete vidiet, ze vela zmien je specifickych pre cloveka
• Vynimkou je slon, niektore zmeny v slonovi su sposobene nizkou kvalitou sekvencie. Ked pomocou nastroja Convert na hornej liste premapujete do najnovsej verzie ludskeho genomu (hg19), uvidite, ze aj v najnovsej verzii genomu slona su mnohe zmeny, nechyba tam uz vsak cast sekvencie, ako vo verzii pouzitej v hg18.

Prakticka ukazka tvorby stromov

• V UCSC browseri mozeme ziskavat viacnasobne zarovnania jednotlivych genov (nukleotidy alebo proteiny). Nasledujuci postup nemusite robit, subor si stiahnite tu: http://compbio.fmph.uniba.sk/vyuka/mbi/poznamky/cb06/cb06-aln.fa
  • UCSC browseri si pozrieme usek ludskeho genomu chr6:136,214,527-136,558,402 s genom PDE7B (phosphodiesterase 7B)
  • Na modrej liste zvolime Tables, v nej RefSeq genes, zaklikneme Region: position, a Output fomat: CDS FASTA alignment a stlacime Get output
  • Na dalsej obrazovke zaklikneme show nucleotides. Z primatov zvolime chimp, rhesus, tarsier, z inych cicavcov mouse, rat, dog, elephant a z dalsich organizmov opposum, platypus, chicken, lizard, stlacime Get output.
  • Vystup ulozime do suboru, z mien sekvencii zmazeme spolocny prefix NM_018945_, pripadne celkovo prepiseme mena na anglicke nazvy

• Skusme zostavit strom na stranke http://mobyle.pasteur.fr/cgi-bin/portal.py
• Pouzijeme program quicktree, metodu neighbor joining, bootstrap 100
• Na zobrazenie stromu vysledok dalej prezenieme cez zobrazovacie programy drawtree alebo newicktops (zvolit v menu pri tlacidle further analysis)
  • Vysledok z drawtree, nezakoreneny, nezobrazuje bootstrap hodnoty
  • Vysledok z newicktops, zakoreneny na nahodnom mieste (nie spravne) zobrazuje bootstrap hodnoty
  • v drawtree sme nastavili sme formát výstupu MS-Windows Bitmap a X,Y resolution aspoň 1000, v newicktops sme nastavili show bootstrap values
• "Spravny strom" [8] v nastaveniach Conservation track-u v UCSC browseri (podla clanku Murphy WJ, Eizirik E, O'Brien SJ, Madsen O, Scally M, Douady CJ, Teeling E, Ryder OA, Stanhope MJ, de Jong WW, Springer MS. Resolution of the early placental mammal radiation using Bayesian phylogenetics. Science. 2001 Dec 14;294(5550):2348-51.)
• Nas strom ma long branch attraction (zle postavenie hlodavcov, ktori maju dlhu vetvu aj slona, co moze byt zapricene sekvenovacimi chybami).
• Ine programy, ktore mozete skusit na mobyle
  • phyml: metoda maximalnej vierohodnosti (daju sa nastavit detaily modelu, bootstraps, ktory ale moze dost dlho trvat, typy operacii na strome pri heuristickom hladani najlepsieho stromu)
  • dnapars alebo protpars na parsimony
  • viacnasobne zarovnanie pomocou clustalw alebo modernejsou alternativou muscle
  • Ak chcete skusat zarovnania, zacnite z nezarovnanych sekvencii: http://compbio.fmph.uniba.sk/vyuka/mbi/poznamky/cb06/cb06-seq.fa

Bootstrap

• Nahodne vyberieme niektore stlpce zarovnania, zostrojime strom vybranou metodou
• Cele to opakujeme vela krat
• Znacime si, kolkokrat sa ktora vetva opakuje v stromoch, ktore dostavame
  • Pri nezakoreneneych stromoch je vetva rozdelenie listov na dve skupiny
• Nakoniec zostavime strom s casto sa vyskytujucich hran (napr tych, ktore su vo viac ako 50% stromov)
• Bootstrap hodnoty nam daju urcity odhad spolahlivosti, hlavne ak mame celkovo malo dat (kratke zarovnanie)
• Ak vsak data velmi dobre nezodpovedaju vybranej metode tak aj pre zly strom mozeme dostat vysoky bootstrap

Evolúcia 2, cvičenia pre biológov

Fylogeneticke stromy

• Terminologia: zakoreneny/nezakoreneny strom, hrana, vrchol, list, vnutorny vrchol, koren
• Ak zmazeme hranu, strom sa rozpadne na dve suvisle casti. Strom mozeme teda povazovat za sadu biparticii.
• Ak mame zakoreneny strom s n listami a kazdy vnutorny vrchol ma 2 deti, bude mat n-1 vnutornych vrcholov (preco?)
• V nezakorenenom strome s n listami, ak ma kazdy vrchol 3 susedov, bude mat n-2 vnutornych vrcholov
• Kolkymi sposobmi mozeme zakorenit nezakoreneny strom?
• Ak nakreslime zakoreneny strom obvyklym sposobom, listy su usporiadane zhora nadol. Kolko roznych poradi listov vieme dostat roznym zakreslovanim toho isteho stromu? Kolko je vobec roznych poradi listov, ak neberime do uvahy strom?

• Ine pouzitie stromov v informatike: uvidime hierarchicke clustrovanie, bayesovske siete, ale tiez efektivne datove struktury

Fitchov algoritmus

• Parsimony/uspornost
  • Vstup: fylogeneticky strom, 1 stlpec zarovnania (jedna baza v kazdom liste stromu)
  • Vystup: priradenie baz predkom minimalizujuce pocet substitucii
• Priklad - obr 1

• Uvazujme, co vieme povedat o strome s dvoma susednymi listami vo vacsom strome (oznacenie: obr. 2, listy v1 a v2, hrany do listov e1, e2, ich predok v3, hrana z v3 vyssie e3).
• Ak oba listy maju bazu rovanku bazu, napr. A, predok v3 v optimalnom rieseni bude urcite mat bazu A
  • Dokaz sporom: nech to tak nie je, nech optimalne riesenie ma nejaku inu bazu, napr. C. Vymenme v tomto rieseni toto C za A. Moze nam pribudnut jedna mutacia na hrane e3, ale ubudnu dve na hranach e1 a e2. Tym celkova cena riesenia klesne o 1, takze nebolo optimalne.
• Ak tieto dva listy maju rozne bazy, napr. A a C, tak existuje optimalne riesenie, ktore ma v predkovi v3 bazu A alebo C.
  • Dokaz: vezmime optimalne riesenie. Ak ma v3 bazu A alebo C, tvrdenie plati. Ak ma v3 nejaku inu bazu, napr T, mozeme ju vymenit napr. za A, ci mozno pribudne jedna muracia na e3 ale urcite ubudne mutaci na e1. Teda celkovy pocet mutacii sa nezvysi a nase nove riesenie je stale optimalne. Pozor, vo vseobecnosti nevieme povedat, ci ma v3 mat bazu A alebo C. V niektorych pripadoch su optimalne obe, v niektorych len jedna z nich.

• Fitchov algoritmus 1971
• Kazdemu vrcholu v priradime mnozinu baz M(v)
• M(v) pocitame od listov smerom ku korenu
• Pre list v bude M(v) obsahovat bazu v tomto liste
• Uvazujme vnutorny vrchol v s detmi x a y. Mame uz spocitane M(x) a M(y), chceme M(v)
• Ak M(x) a M(y) maju nejake spolocne bazy, vsetky tieto spolocne bazy dame do M(v), t.j. $M(v)=M(x)\cap M(y)$
• Ak M(x) a M(y) nemaju spolocne bazy, do M(v) dame vsetky bazy z M(x) aj M(v), t.j. $M(v)=M(x)\cup M(y)$

V tomto pripade pocet mutacii vzrastie o jedna

• Ked mame M(v) spocitane pre vsetky vrcholy, ideme od korena smerom k listom a vyberieme vzdy jednu bazu z M(v).
• Ak sme vybrali pre rodica bazu x a x je v M(v), zvolime x aj pre v, inak zvolime lubovolnu bazu z M(v).
• Priklad algoritmu na obr 3

Opakovanie pravdepodobnostnych modelov

Ake sme doteraz videli modely

• E-value v BLASTe: nahodne vygenerujeme databazu a dotaz (query), kolko bude v priemere medzi nimi lokalnych zarovnani so skore aspon T?
• Hladanie genov: model generujuci sekvenciu+anotaciu naraz (parametre nastavene na znamych genoch). Pre danu sekvenciu, ktora pravdedpodobnost je najpravdepodobnejsia?
• Evolucia, Jukes-Cantorov model: model generujuci stlpec zarovnania. Nezname parametre: strom, dlzky hran. Pre danu sadu stlpcov zarovnania, ktore parametre povedu k najvacsej pravdepodobnosti? $\max _{{param}}\Pr(data|param)$
• Trochu detailov: pravdepodobnost zmeny/nezmeny na hrane dlzky t: $P(A|A,t)=(1+3e^{{-{\frac {4}{3}}t}})/4$, $P(C|A,t)=(1-e^{{-{\frac {4}{3}}t}})/4$
• Ak pozname ancestralne sekvencie, vieme spocitat pravdepodobnost dat
• Ancestralne sekvencie su nahodne premenne, ktore nas nezaujimaju: marginalizujeme ich (uvazujeme vsetky ich mozne hodnoty)

Zlozitejsie evolucne modely

• Jukes-Cantorov model uvazuje len dlzku hrany udanu v priemernom pocte substitucii (vratane tych, ktore nevidime, kvoli tomu, ze su dve na tom istom mieste)
• Nie vsetky substitucie sa deju rovnako casto: tranzicie (v ramci pyrimidinov T<->C, v ramci purinov A<->G) su pravdepodobnjesie ako transverzie (A,G)<->(C,T)
• Nie vsekty nukleotidy sa v danom genome objavuju rovnako casto (napr. mitochondrialne genomy velmi male C,G)
• Tieto javy zachytava HKY model
• Matica rychlosti

$\left({\begin{array}{cccc}-\mu _{A}&\beta \pi _{C}&\alpha \pi _{G}&\beta \pi _{T}\\\beta \pi _{A}&-\mu _{C}&\beta \pi _{G}&\alpha \pi _{T}\\\alpha \pi _{A}&\beta \pi _{C}&-\mu _{G}&\beta \pi _{T}\\\beta \pi _{A}&\alpha \pi _{C}&\beta \pi _{G}&-\mu _{T}\\\end{array}}\right)$

• $\kappa =\alpha /\beta$ je pomer rychlosti, ktorymi sa deju tranzicie vs. transverzie
• $\pi _{j}$ je frekvencia bazy $j$ v sekvencii
• Rychlost, ako sa deje substitucia z X do Y je sucin pravdepodobnosti Y a faktoru, ktory zavisi od toho, ci ide o tranziciu alebo transverziu
• Sucet kazdeho stlpca matice ma byt 0, t.j. $\mu _{A}=\beta \pi _{A}+\alpha \pi _{A}+\beta \pi _{A}$
• Matica sa znormalizuje tak, aby priemerny pocet substitucii za jednotku casu bol 1
• Matica ma styri parametre: $\kappa$ a tri frekvencie (stvrta musi doplnit do 1) plus dlzka hrany
• Zlozitejsi model lepsie zodpoveda skutocnym procesom, ale na odhad viac parametrov potrebujeme viacej dat.
• Existuju metody, ktore pre dany cas t z matice spocitaju pravdepodobnost, ze baza X zmutuje na bazu Y Pr(Y|X,t)
• Napr. pre velmi velmi maly cas $\epsilon$ mame $\Pr(C|A,\epsilon )$ je zhruba $\epsilon \beta \pi _{C}$
• Pre rozumne dlhe casy toto neplati, preto sa pouzivaju algebraicke metody, ktore beru do uvahy moznost viacerych substitucii na tom istom mieste
• Je aj vela inych modelov s mensim alebo vacsim poctom parametrov

Evolúcia, cvičenia pre informatikov

Felsensteinov algoritmus 1981

• Mame dany strom T s dlzkami hran a bazy v listoch (jeden stlpec zarovnania) a maticu rychlosti R. Spocitajme pravdepodobnost, ze z modelu dostaneme prave tuto kombinaciu baz v listoch.
• Nech X_v je premenna reprezentujuca bazu vo vrchole v a nech x_v je konkretna baza v liste v. Nech listy su 1..n a vnut. vrcholy n+1..2n-1. Nech dlzka hrany z v do rodica je t_v. Nech P(a|b,t) je pravdepodobnost, ze b sa zmeni na a za cas t (spocitame z matice R, ako nabuduce). Nech q_a je pravdepodobnost bazy a v koreni (ekvilibrium matice R)
• Napr. v Jukes-Cantorovom modeli $P(A|A,t)=(1+3e^{{-{\frac {4}{3}}t}})/4$, $P(C|A,t)=(1-e^{{-{\frac {4}{3}}t}})/4$

• Chceme pravdepodobnost $P(X_{1}=x_{1},X_{2}=x_{2},\dots X_{n}=x_{n}|T,R)=\sum _{{x_{{n+1}}\dots x_{{2n-1}}}}P(X_{1}=x_{1}\dots X_{{2n-1}}=x_{{2n-1}}|T,R)$
• Spocitame rychlejsie dynamickym programovanim.
• Nech A[v,a] je pravdedpodobnost dat v podstrome s vrcholom v ak X_v=a
• A[v,a] pocitame od listov ku korenu
• v liste A[v,a] = [a=x_v]
• Vo vnut. vrchole s detmi x a y mame $A[v,a]=\sum _{{b,c}}A[x,b]A[y,c]P(b|a,t_{x})P(c|a,t_{y})$
• Celkova pravdepodobnost je $P(X_{1}=x_{1},X_{2}=x_{2},\dots X_{n}=x_{n}|T,R)=\sum _{a}A[r,a]q_{a}$ pre koren r.

Zlozitost, zlepsenie

• Zlozitost $O(n|\Sigma |^{3})$
• Pre nebinarne stromy exponencialne
• Zlepsenie $A[v,a]=(\sum _{{b}}A[x,b]P(b|a,t_{x}))(\sum _{c}A[y,c](c|a,t_{y}))$
• Zlozitost $O(n|\Sigma |^{2})$ aj pre nebinarne stromy

Chybajuce data

• Ak v niektorom liste mame neznamu bazu N, nastavime A[v,a]=1
• Podobne sa spracovavaju medzery v zarovnani, aj ked mohli by sme mat aj model explicitne ich modelujuci

Aposteriorna pravdepodobnost

• Co ak chceme spocitat pravdepodobnost P(X_v=a|X_1=x_1, X_2=x_2,\dots X_n=x_n,T,R)?
• Potrebujeme B[v,a]=pravdpodobnost dat ak podstrom v nahradim listom s bazou a.
• B[v,a] pocitame od korena k listom
• V koreni B[v,a] = q_a
• Vo vrchole v s rodicom u a surodencom x mame $B[v,a]=\sum _{{b,c}}B[u,b]A[x,c]P(a|b,t_{v})P(c|b,t_{v})$
• Ziadana pravdepodobnost je $B[v,a]A[v,a]/P(X_{1}=x_{1},X_{2}=x_{2},\dots X_{n}=x_{n}|T,R)$

Fylogeneticke HMM

• Spojenie HMM a fylogenetickeho stromu
• Konecne mnozina stavov
• $a_{{u,v}}$: prechodová pravdepodobnosť zo stavu $u$ do stavu $v$
• $\pi _{{u}}$: pravdepodobnosť, že začneme v stave $u$
• Namiesto emisnej pravdepodobnosti mame pre kazdy stav dlzky hran, maticu rychlosti R

Priklad:

• PhastCons: 2 stavy: zachovana a nezachovana sekvencia
• Komparativne hladanie genov

Opakovanie: Viterbiho algoritmus pre HMM

• Podproblém $V[i,u]$: Pravdepodobnosť najpravdepodobnejšej cesty končiacej po $i$ krokoch v stave $u$, pričom vygeneruje $S_{1}S_{2}\dots S_{i}$.
• $V[i,u]=\max _{v}V[i-1,v]a_{{v,u}}e_{{u,x_{i}}}$

Modifikacia pre fylogeneticke HMM:

• Pre kazdy stlpec zarovnania i a stav u spustime Felsensteinov algoritmus
• Vyslednu pravdepodobnost pouzijeme namiesto $e_{{u,x_{i}}}$
• Zlozitost pre k stavov, dlzku zarovnania n, pocet organizmov m: O(nmk + nk^2)

Parove HMM

• Emituje dve sekvencie
• V jednom kroku moze emitovat:
  • pismenka v oboch sekvenciach naraz
  • pismenko v jednej skevencii
  • pismenko v druhej sekvencii

Priklad: HMM s jednym stavom v, takym, ze

• $e_{{v,x,x}}=p_{1}$
• $e_{{v,x,y}}=p_{2}(x\neq y)$,
• $e_{{v,x,-}}=p_{3}$,
• $e_{{v,-,x}}=p_{3}$
• tak, aby sucet emisnych pravdepodobnosti bol 1
• Co reprezentuje najpravdepodobnejsia cesta v tomto HMM?

Zlozitejsi HMM: tri stavy M, X, Y, uplne navzajom poprepajane

• $e_{{M,x,x}}=p_{1}$
• $e_{{M,x,y}}=p_{2}(x\neq y)$,
• $e_{{X,x,-}}=1/4$,
• $e_{{Y,-,y}}=1/4$,
• Co reprezentuje najpravdepodobnejsia cesta v tomto HMM?

Viterbiho algoritmus pre parove HMM

• V[i,j,u] = pravdepodobnost najpravdepodobnejsej postupnosti stavov, ktora vygeneruje x1..xi a y1..yj a skonci v stave u
• $V[i,j,u]=\max _{w}\left\{{\begin{array}{l}V[i-1,j-1,w]\cdot a_{{w,u}}\cdot e_{{u,x_{i},y_{j}}}\\V[i-1,j,w]\cdot a_{{w,u}}\cdot e_{{u,x_{i},-}}\\V[i,j-1,w]\cdot a_{{w,u}}\cdot e_{{u,-,y_{j}}}\\\end{array}}\right.$
• Casova zlozitost O(mnk^2) kde m a n su dlzky vstupnych sekvencii, k je pocet stavov

Ako by sme spravili parove HMM na hladanie genov v dvoch sekvenciach naraz?

• Predpokladajme rovnaky pocet exonov
• V exonoch medzery len cele kodony (oboje zjednodusuje)
• Inde hocijake medzery

Evolúcia 2, cvičenia pre informatikov

Substitucne modely - odvodenie

• Nech $P(b|a,t)$ je pravdepodobnost, ze ak sme zacneme s bazou a, tak po case t budeme mat bazu b.
• Pre dane t mozeme take pravdepodobnosti usporiadat do matice 4x4 (ak studujeme DNA), kde $S(t)_{{a,b}}=P(b|a,t)$
• Intuitivne cim vacsie t, tym vacsia pravdepodobnost zmeny; $S(0)=I$ (jednotkova matica), $S(t),t\rightarrow \infty$ ma vsetky riadky rovnake, napr. 1/4, 1/4, 1/4, 1/4
• Ak mame matice pre casy $t_{1}$ a $t_{2}$, vieme spocitat maticu pre cas $t_{1}+t_{2}$: $P(b|a,t_{1}+t_{2})=\sum _{x}P(x|a,t_{1})\cdot P(b|x,t_{2})$ a teda v maticovej notacii $S(t_{1}+t_{2})=S(t_{1})S(t_{2})$. Takyto model nazyvame multiplikativny a predpoklada, ze pravdepodobnost mutacie zavisi len od aktualnej bazy, nie od minulych stavov.
• Ak by sme uvazovali iba diskretne (celocislene) casy, stacilo by nam urcit iba $S(1)$ a vsetky ostatne casy dostaneme umocnenim tejto matice. Je vsak elegantnejsie mat $S(t)$ definovane aj pre realne t.

• Jukes-Cantor-ov model evolucie predpoklada, ze vsetky substitucie su rovnako pravdepodobne.

• Pre velmi maly cas mame maticu

$S(\epsilon )=\left({\begin{array}{cccc}1-3p&p&p&p\\p&1-3p&p&p\\p&p&1-3p&p\\p&p&p&1-3p\\\end{array}}\right)$

kde $p$ je tiez velmi male cislo.

• Pre cas $2\epsilon$ dostavame maticu

$S(2\epsilon )=S(\epsilon ^{2})=\left({\begin{array}{cccc}1-6p+12p^{2}&2p-4p^{2}&2p-4p^{2}&2p-4p^{2}\\\dots \end{array}}\right)$

• Ale cleny s $p^{2}$ su ovela mensie ako cleny s $p$, takze tato matica je priblizne

$S(2\epsilon )=S(\epsilon ^{2})\approx \left({\begin{array}{cccc}1-6p&2p&2p&2p\\\dots \end{array}}\right)$

• Vytvorme si teraz maticu rychlosti (rate matrix)

$R=\left({\begin{array}{cccc}-3\alpha &\alpha &\alpha &\alpha \\\alpha &-3\alpha &\alpha &\alpha \\\alpha &\alpha &-3\alpha &\alpha \\\alpha &\alpha &\alpha &-3\alpha \\\end{array}}\right)$

• Dostavame, ze pre velmi male casy plati $S(\epsilon )\approx I+R\epsilon$ ($p$ pouzivane vyssie by malo hodnotu $\alpha /\epsilon$).
• $S(t+\epsilon )=S(t)S(\epsilon )\approx S(t)(I+R\epsilon )$ a teda $(S(t+\epsilon )-S(t))/\epsilon \approx S(t)R$ a v limite pre $\epsilon \rightarrow o$ dostavame $S'(t)=S(t)R$ (diferencialne rovnice, pociatocny stav $S(0)=I$).
• Ak diagonalne prvky $S(t)$ oznacime $r(t)$ a nediagnoalne $s(t)$, dostavame, ze diagonalny prvok $S(t)R$ je $-3\alpha r(t)+3\alpha s(t)$ a nediagonalny $-\alpha s(t)+\alpha r(t)$. Takze dostavame diferencialne rovnice $r'(t)=-3\alpha (s(t)-r(t))$ a $s'(t)=\alpha (r(t)-s(t))$.
• Overme, ze riesenim tejto rovnice je $r(t)=(1+3e^{{-4\alpha t}})/4$ a $s(t)=(1-e^{{-4\alpha t}})/4$:

Zderivujeme $r'(t)=-3\alpha e^{{-4\alpha t}}$ a $s'(t)=\alpha e^{{-4\alpha t}}$, dosadime do rovnic.

• Takze mame maticu:

$S(t)=\left({\begin{array}{cccc}(1+3e^{{-4\alpha t}})/4&(1-e^{{-4\alpha t}})/4&(1-e^{{-4\alpha t}})/4&(1-e^{{-4\alpha t}})/4\\\dots \end{array}}\right)$

• V case $t$ je pravdepodobnost, ze uvidime zmenenu bazu $D_{t}={\frac {3}{4}}(1-e^{{-4\alpha t}})$ a teda ak v skutocnosti vidime $d\%$ zmenenych baz, vieme spatne zratat t, ktore by hodnote $D_{t}=d$ prinalezalo.
• Aby sme nemali naraz aj $\alpha$ aj $t$, zvykneme maticu R normalizovat tak, aby priemerny pocet substitucii za jednotku casu bol 1. V pripade Jukes-Cantorovho modelu je to ked $\alpha =1/3$.
• Dostavame teda vzorec pre vzdialenost, ktory sme videli na prednaske $t=-{\frac {3}{4}}\log \left(1-{\frac {4}{3}}d\right)$
• Preco sme ten vzorec odvodili takto? V skutocnosti chceme najst najvierohodnejsiu hodnotu t, t.j. taku, pre ktore hodnota P(data|t) bude najvacsia. Zhodou okolnosti vyjde takto.

• V praxi sa pouzivaju komplikovanejsie substitucne modely, ktore maju vseobecnejsiu maticu rychlosti R
• Kimurov model napr. zachytava, ze puriny sa castejsie menia na ine puriny (A a G) a pyrimidiny na ine pyrimidiny (C a T) ma dva parametre:

$R=\left({\begin{array}{cccc}-2\beta -\alpha &\beta &\alpha &\beta \\\beta &-2\beta -\alpha &\beta &\alpha \\\alpha &\beta &-2\beta -\alpha &\beta \\\beta &\alpha &\beta &-2\beta -\alpha \\\end{array}}\right)$

• HKY model (Hasegawa, Kishino & Yano) tiez umoznuje rozne pravdepodobnosti A, C, G a T v ekvilibriu.
• Vo vseobecnosti pre rate matrix $R$ dostavame $S(t)=e^{{Rt}}$. Ak R diagonalizujeme (urcite sa da pre symetricke R) $R=UDU^{{-1}}$, kde D je diagonalna matica (na jej diagonale budu vlastne hodnoty R), tak $e^{{Rt}}=Ue^{{Dt}}U^{{-1}}$, t.j. exponencialnu funkciu uplatnime iba na prvky na uhlopriecke matice D.

Proteíny, cvičenia pre biológov

• Toto cvičenie je z časti inšpirované stránkou [9]
• Pozrieme sa na enzým Bis(5'-adenosyl)-triphosphatase
• Nájdime ho na stránke http://www.uniprot.org/ pod názvom FHIT_HUMAN
  • Pozrime si podrobne jeho stránku, ktoré časti boli predpovedané bioinformatickými metódami z prednášky?
  • Všimnime si Pfam doménu a pozrime si jej stránku, do akej super-rodiny (klanu) patrí?

• Tento enzým je vzdialene podobný na enzým galactose-1-phosphate uridylyltransferase (GALT/GAL7)
• Skúsme nájsť túto podobnosť v BLASTe: http://blast.ncbi.nlm.nih.gov/ v časti proteíny, zvoľme databázu Swissport, ako Query zadajme Accesion nášho proteínu P49789
• GAL gén (konkrétne GAL7_HAEIN) sa nachádza medzi výsledkami, ale s vysokou E-value
• Odložme si tieto výsledky a v novom okne spustíme PSI-BLAST, dve iterácie
• Aká je E-value teraz? V čom sa líšia zarovnania?

• Nájdime v Uniprote proteín GAL7_HAEIN.
• Ak sa pozrieme na jeho Pfam domény, v akom sú klane?

Evolúcia 3 a proteíny, cvičenia pre informatikov

Substitucne modely - odvodenie

• Nech $P(b|a,t)$ je pravdepodobnost, ze ak sme zacneme s bazou a, tak po case t budeme mat bazu b.
• Pre dane t mozeme take pravdepodobnosti usporiadat do matice 4x4 (ak studujeme DNA), kde $S(t)_{{a,b}}=P(b|a,t)$
• Intuitivne cim vacsie t, tym vacsia pravdepodobnost zmeny; $S(0)=I$ (jednotkova matica), $S(t),t\rightarrow \infty$ ma vsetky riadky rovnake, napr. 1/4, 1/4, 1/4, 1/4
• Ak mame matice pre casy $t_{1}$ a $t_{2}$, vieme spocitat maticu pre cas $t_{1}+t_{2}$: $P(b|a,t_{1}+t_{2})=\sum _{x}P(x|a,t_{1})\cdot P(b|x,t_{2})$ a teda v maticovej notacii $S(t_{1}+t_{2})=S(t_{1})S(t_{2})$. Takyto model nazyvame multiplikativny a predpoklada, ze pravdepodobnost mutacie zavisi len od aktualne bazy, nie od minulych stavov.
• Ak by sme uvazovali iba diskretne (celocislene) casy, stacilo by nam urcit iba $S(1)$ a vsetky ostatne casy dostaneme umocnenim tejto matice. Je vsak elegantnejsie mat $S(t)$ definovane aj pre realne t.

• Jukes-Cantor-ov model evolucie predpoklada, ze vsetky substitucie su rovnako pravdepodobne.

• Pre velmi maly cas mame maticu

$S(\epsilon )=\left({\begin{array}{cccc}1-3p&p&p&p\\p&1-3p&p&p\\p&p&1-3p&p\\p&p&p&1-3p\\\end{array}}\right)$

kde $p$ je tiez velmi male cislo.

• Pre cas $2\epsilon$ dostavame maticu

$S(2\epsilon )=S(\epsilon )^{2}=\left({\begin{array}{cccc}1-6p+12p^{2}&2p-4p^{2}&2p-4p^{2}&2p-4p^{2}\\\dots \end{array}}\right)$

• Ale cleny s $p^{2}$ su ovela mensie ako cleny s $p$, takze tato matica je priblizne

$S(2\epsilon )=S(\epsilon ^{2})\approx \left({\begin{array}{cccc}1-6p&2p&2p&2p\\\dots \end{array}}\right)$

• Vytvorme si teraz maticu rychlosti (rate matrix)

$R=\left({\begin{array}{cccc}-3\alpha &\alpha &\alpha &\alpha \\\alpha &-3\alpha &\alpha &\alpha \\\alpha &\alpha &-3\alpha &\alpha \\\alpha &\alpha &\alpha &-3\alpha \\\end{array}}\right)$

• Dostavame, ze pre velmi male casy plati $S(\epsilon )\approx I+R\epsilon$ ($p$ pouzivane vyssie by malo hodnotu $\alpha /\epsilon$).
• $S(t+\epsilon )=S(t)S(\epsilon )\approx S(t)(I+R\epsilon )$ a teda $(S(t+\epsilon )-S(t))/\epsilon \approx S(t)R$ a v limite pre $\epsilon \rightarrow o$ dostavame $S'(t)=S(t)R$ (diferencialne rovnice, pociatocny stav $S(0)=I$).
• Ak diagonalne prvky $S(t)$ oznacime $r(t)$ a nediagnoalne $s(t)$, dostavame, ze diagonalny prvok $S(t)R$ je $-3\alpha r(t)+3\alpha s(t)$ a nediagonalny $-\alpha s(t)+\alpha r(t)$. Takze dostavame diferencialne rovnice $r'(t)=-3\alpha (s(t)-r(t))$ a $s'(t)=\alpha (r(t)-s(t))$.
• Overme, ze riesenim tejto rovnice je $r(t)=(1+3e^{{-4\alpha t}})/4$ a $s(t)=(1-e^{{-4\alpha t}})/4$:

Zderivujeme $r'(t)=-3\alpha e^{{-4\alpha t}}$ a $s'(t)=\alpha e^{{-4\alpha t}}$, dosadime do rovnic.

• Takze mame maticu:

$S(t)=\left({\begin{array}{cccc}(1+3e^{{-4\alpha t}})/4&(1-e^{{-4\alpha t}})/4&(1-e^{{-4\alpha t}})/4&(1-e^{{-4\alpha t}})/4\\\dots \end{array}}\right)$

• V case $t$ je pravdepodobnost, ze uvidime zmenenu bazu $D_{t}={\frac {3}{4}}(1-e^{{-4\alpha t}})$ a teda ak v skutocnosti vidime $d\%$ zmenenych baz, vieme spatne zratat t, ktore by hodnote $D_{t}=d$ prinalezalo.
• Aby sme nemali naraz aj $\alpha$ aj $t$, zvykneme maticu R normalizovat tak, aby priemerny pocet substitucii za jednotku casu bol 1. V pripade Jukes-Cantorovho modelu je to ked $\alpha =1/3$.
• Dostavame teda vzorec pre vzdialenost, ktory sme videli na prednaske $t=-{\frac {3}{4}}\log \left(1-{\frac {4}{3}}d\right)$
• Preco sme ten vzorec odvodili takto? V skutocnosti chceme najst najvierohodnejsiu hodnotu t, t.j. taku, pre ktore hodnota P(data|t) bude najvacsia. Zhodou okolnosti vyjde takto.

• V praxi sa pouzivaju komplikovanejsie substitucne modely, ktore maju vseobecnejsiu maticu rychlosti R
• Kimurov model napr. zachytava, ze puriny sa castejsie menia na ine puriny (A a G) a pyrimidiny na ine pyrimidiny (C a T) ma dva parametre:

$R=\left({\begin{array}{cccc}-2\beta -\alpha &\beta &\alpha &\beta \\\beta &-2\beta -\alpha &\beta &\alpha \\\alpha &\beta &-2\beta -\alpha &\beta \\\beta &\alpha &\beta &-2\beta -\alpha \\\end{array}}\right)$

• HKY model (Hasegawa, Kishino & Yano) tiez umoznuje rozne pravdepodobnosti A, C, G a T v ekvilibriu.
• Vo vseobecnosti pre rate matrix $R$ dostavame $S(t)=e^{{Rt}}$. Ak R diagonalizujeme (urcite sa da pre symetricke R) $R=UDU^{{-1}}$, kde D je diagonalna matica (na jej diagonale budu vlastne hodnoty R), tak $e^{{Rt}}=Ue^{{Dt}}U^{{-1}}$, t.j. exponencialnu funkciu uplatnime iba na prvky na uhlopriecke matice D.

Protein threading a integer linear programming

• Obcas chceme najst optimalne riesenie nejakeho NP-tazkeho problemu
• Jedna moznost je previest na iny NP tazky problem, pre ktory existuju pomerne dobre prakticke programy, napriklad integer linear programming (ILP)

Prakticke programy na NP tazke problemy

• najdu optimalne riesenie, mnohe instancie zrataju v rozumnom case, ale mozu bezat aj velmi dlho
• CPLEX [10] a Gurobi [11] komercne baliky na ILP, akademicka licencia zadarmo
• SYMPHONY v projekte COIN-OR [12] open source
• Minisat [13] open source SAT solver
• Concorde TSP solver [14] - riesi problem obchodneho cestujuceho so symetrickymi vzdialenostami, zadarmo na akademicke ucely
  • Pre zaujimavost: TSP art [15]

ILP

Linearny program:

• Mame realne premenne x_1...x_n, minimalizujeme nejaku ich linearnu kombinaciu $\sum _{i}a_{i}x_{i}\,$ kde a_i su dane vahy.
• Mame tiez niekolko podmienok v tvare linearnych rovnosti alebo nerovnosti, napr. $\sum _{i}b_{i}x_{i}\leq c$
• Hladame teda hodnoty premennych, ktore minimalizuju cielovu sumu, ale pre ktore platia vsetky podmienky

Integer linear program

• Program, v ktorom vsetky premenne musia mat celociselne hodnoty, alebo dokonca povolime oba hodnoty 0 a 1.

Ako zapisat (NP-tazke) problemy ako ILP

Knapsack

• Problem: mame dane predmety s hmotnostami w_1..w_n a cenami c_1..c_n, ktore z nich vybrat, aby celkova hmotnost bola najviac T a cena bola co najvyssia?
• Pouzijeme binarne premenne x_1..x_n, kde x_i = 1 prave vtedy ked sme zobrali i-ty predmet.
• Chceme maximalizovat $\sum _{i}c_{i}x_{i}\,$
• za podmienky ze $\sum _{i}w_{i}x_{i}\leq T$

Set cover:

• Mame n mnozin S_1...S_n nad mnozinou {1...m}. Chceme vybrat co najmensiu pomnozinu A tak, aby kazde S_i obsahovalo aspon jeden prvok z A
• Binarne premenne x_i=1 ak vyberieme i-ty prvok
• Chceme maximalizovat $\sum _{i}x_{i}\,$
• za podmienky ze pre kazde i plati $\sum _{{j\in S_{i}}}x_{j}\geq 1$

Protein threading

• Ciel: protein A ma znamu sekvenciu aj strukturu, protein B iba sekvenciu. Chceme zarovnat proteiny A a B, pricom budeme brat do uvahy znamu strukturu, t.j. ak su dve amino kyseliny blizko v A tak ich ekvivalenty v B by mali byt "kompatibilne".
• Tento problem chceme riesit tak, ze v strukture A urcime nejake jadra, ktore by v evolucii mali zostat zachovane bez inzercii a delecii a v rovnakom poradi. Tieto jadra su oddelene sluckami, ktorych dlzka sa moze lubovolne menit a ktorych zarovnania nebudeme skorovat.
• Formulacia problemu: Mame danu sekvenciu B=b1..bn, dlzky m jadier c_1...c_m a skorovacie tabulky E_ij, ktora vyjadruje, ako dobre bj..b_{j+c_i-1} sedi do sekvencie jadra i a E_ijkl ktora vyjadruje, ako dobre by jadra i a k intergaovali, keby mali sekvencie zacinajuce v B na poziciach j a l. Uloha je zvolit polohy jadier x_1<x_2<...<x_m tak, aby sa ziadne dve jadra neprekryvali a aby sme dosiahli najvyssie skore.
• Poznamka: nevraveli sme, ako konkretne zvoli jadra a skorovacie tabulky, co je modelovaci, nie algoritmicky problem (mozeme skusit napr. nejake pravdepodonostne modely)

Protein threading ako ILP

• Premenne v programe:
  • x_ij=1 ak je zaciatok i-teho jadra zarovnane s b_j
  • y_ijkl=1 ak je zaciatok i-teho jadra na b_j a zaciatok k-teho na b_l (i<k, j<l)
• Chceme maximalizovat $\sum E_{{ij}}x_{{ij}}+\sum E_{{ijkl}}y_{{ijkl}}$
• Podmienky:
  • $\sum _{j}x_{{ij}}=1\,$ pre kazde i
  • $x_{{il}}+x_{{i+1,k}}=\leq 1$ pre vsetky i,k,l, kde k<l+c_i
  • $y_{{ijkl}}\leq x_{{ij}}$ pre vsetky i,j,k,l, kde i<k, j<l
  • $y_{{ijkl}}\leq x_{{kl}}$ pre vsetky i,j,k,l, kde i<k, j<l
  • $y_{{ijkl}}\geq x_{{ij}}+x_{{kl}}$ pre vsetky i,j,k,l, kde i<k, j<l

Na zamyslenie:

• Aka bude velkost programu ako funkcia n a m?
• Co ak nie vsetky jadra navzajom interaguju? Mozeme na velkosti programu usetrit?
• Preco asi vobec autori zaviedli jadra a ako by sme zmenili program, ak by sme chceli uvazovat kazdu aminokyselinu zvlast?

Proteíny, cvičenia pre informatikov

Úvod do bioinformatických databáz a on-line nástrojov

NCBI, Genbank, Pubmed, blast

• National Center for Biotechnology Information http://www.ncbi.nlm.nih.gov/
• Zhromazduje verejne pristupne data z molekularnej biologie
• Mozeme hladat klucove slova v rzonych databazach
  • Pubmed: databaza clankov, napr. najdime phastcons
  • Books: najdime Alu
  • Gene: najdime DNA polymerazu
  • BLAST: najdime nasledujucu sekvenciu v genome kurata (zvoľme nucleotide blast, database others a z menu reference genomic sequence, organism chicken (taxid:9031), program blastn)
  • Ide o osekvenovany kusok ludskej mRNA, kde v kuracom genome sme nasli homolog, ake ma dlzku, skore, E-value, % zhdonych baz?

AACCATGGGTATATACGACTCACTATAGGGGGATATCAGCTGGGATGGCAAATAATGATTTTATTTTGAC
TGATAGTGACCTGTTCGTTGCAACAAATTGATAAGCAATGCTTTCTTATAATGCCAACTTTGTACAAGAA
AGTTGGGCAGGTGTGTTTTTTGTCCTTCAGGTAGCCGAAGAGCATCTCCAGGCCCCCCTCCACCAGCTCC
GGCAGAGGCTTGGATAAAGGGTTGTGGGAAATGTGGAGCCCTTTGTCCATGGGATTCCAGGCGATCCTCA
CCAGTCTACACAGCAGGTGGAGTTCGCTCGGGAGGGTCTGGATGTCATTGTTGTTGAGGTTCAGCAGCTC
CAGGCTGGTGACCAGGCAAAGCGACCTCGGGAAGGAGTGGATGTTGTTGCCCTCTGCGATGAAGATCTGC
AGGCTGGCCAGGTGCTGGATGCTCTCAGCGATGTTTTCCAGGCGATTCGAGCCCACGTGCAAGAAAATCA
GTTCCTTCAGGGAGAACACACACATGGGGATGTGCGCGAAGAAGTTGTTGCTGAGGTTTAGCTTCCTCAG
TCTAGAGAGGTCGGCGAAGCATGCAGGGAGCTGGGACAGGCAGTTGTGCGACAAGCTCAGGACCTCCAGC
TTTCGGCACAAGCTCAGCTCGGCCGGCACCTCTGTCAGGCAGTTCATGTTGACAAACAGGACCTTGAGGC
ACTGTAGGAGGCTCACTTCTCTGGGCAGGCTCTTCAGGCGGTTCCCGCACAAGTTCAGGACCACGATCCG
GGTCAGTTTCCCCACCTCGGGGAGGGAGAACCCCGGAGCTGGTTGTGAGACAAATTGAGTTTCTGGACCC
CCGAAAAGCCCCCACAAAAAGCCG

Uniprot

• Prehladnejsi pohlad na proteiny, vela linkov na ine databazy, cast vytvarana rucne
  • Pozrieme sa na enzým Bis(5'-adenosyl)-triphosphatase
  • Nájdime ho na stránke http://www.uniprot.org/ pod názvom FHIT_HUMAN
  • Pozrime si podrobne jeho stránku, ktoré časti boli predpovedané bioinformatickými metódami z prednášky?
  • Všimnime si Pfam doménu a pozrime si jej stránku, do akej super-rodiny (klanu) patrí?

UCSC genome browser

• http://genome.ucsc.edu/

Sekvenovanie

• Hore v modrom menu zvoľte Genomes
• Na ďalšej stránke zvoľte človeka a v menu Assembly zistite, kedy boli pridané posledné dve verzie ľudského genómu (hg18 a hg19)
• Na tej istej stránke dole nájdete stručný popis zvolenej verzie genómu. Pre ktoré chromozómy máme viacero alternatívnych verzií?
• Zadajte región chr21:31,200,000-31,400,000
• Zapnite si track Mapability na "pack" a track RepeatMasker prepnite na "full"
• Mapability: nakoľko sa daný úsek opakuje v genóme a či teda vieme jednoznačne jeho ready namapovať pri použití Next generation sequencing (detaily keď kliknete na linku "Mapability")
• Približne v strede zobrazeného regiónu je pokles mapovateľnosti. Akému typu opakovania zodpovedá? (pozrite track RepeatMasker)
• Zapnite si tracky "Assembly" a "Gaps" a pozrite si región chr2:110,000,000-110,300,000. Aká dlhá je neosekvenovaná medzera (gap) v strede tohto regiónu? Približnú veľkosť môžete odčítať z obrázku, presnejší údaj zistíte kliknutím na čierny obdĺžnik zodpovedajúci tejto medzere (úplne presnú dĺžku aj tak nepoznáme, nakoľko nie je osekvenovaná).

Geny

• Zvolte starsiu verziu ludskeho genomu hg18, ktora ma viac informacii
• Do okienka position zadajte gen MAGEA2B a potom zvolte vyskyt na pozicii cca chrX:151636040-151637735 (ma dva vyskyty)
  • Dostanete sa tam aj touto linkou: [16]
• Ak date 3x zoom out, mozete si vsimnut, ze tento gen ma viacero foriem zostrihu, ktore sa ale lisia iba v 5' UTR
• Vela veci sa mozete dozvediet klikanim na rozne casti broswera: napr, kliknutim na gen si mozete precitat o jeho funkcii, kliknutim na listu ku tracku (lavy okraj obazku) sa dozviete viac o tracku a mozete nastavovat parametre zobrazenia

Komparativna genomika

• V casti Genes and Gene Prediction Tracks zapnite track Pos Sel Genes, ktory obsahuje geny s pozitivnym vyberom (cervenou, pripadne slabsie fialovou a modrou)
• Ked kliknete na cerveny obdlznik pre tento gen, uvidite, v ktorych castiach fylogenetickeho stromu bol detekovany pozitivny vyber

• V casti multiz alignments vidite zarovnania k roznym inym genomom (da sa zapinat, ze ku ktorym). Mozete si pozriet, ako sa uroven zarovnania zmeni ked sa priblizujeme a vzdalujeme (zoom in/zoom out).
• Ked sa priblizite spat na gen MAGEA2B a potom tak, aby ste boli na urovni "base", t.j. zobrazenych cca 100bp, v obdlzniku multiz alignment uvidite zarovnanie s homologickym usekom v inych genomoch. Konkretne v MAGEA2B vidime pomerne dost rozdielov v proteine medzi clovekom a makakom rezus, vdaka ktorym bol zrejme klasifikovany ako pod pozitivnym vyberom.
• V casti conservation by PhyloP vidime graf toho, ako silne su zachovane jednotlive stlpce zarovnania
• Da sa zapnut aj track 28-Way Most Cons, ktory zobrazuje konkretne useky, ktore su najvac konzervovane
  • Ak chceme zistit, kolko percent genomu tieto useky pokryvaju, ideme na modrej liste do casti Tables, zvolime group Comparative genomics, track 28-way most cons, table Mammal alebo Vertebrate, region zvolime genome (v celom genome) a stlacime tlacidlo Summary/statistics

Blat

• Choďte na UCSC genome browser (http//genome.ucsc.edu/), na modrej lište zvoľte BLAT, zadajte túto sekvenciu a hľadajte ju v ľudskom genóme. Akú podobnosť (IDENTITY) má najsilnejší nájdený výskyt? Aký dlhý úsek genómu zasahuje? (SPAN). Všimnite si, že ostatné výskyty sú oveľa kratšie.
• V stĺpci ACTIONS si pomocou Details môžete pozrieť detaily zarovnania a pomocou Browser si pozrieť príslušný úsek genómu.
• V tomto úseku genómu si zapnite track Vertebrate net na full a kliknutím na farebnú čiaru na obrázku pre tento track zistite, na ktorom chromozóme kuraťa sa vyskytuje homologický úsek.
• Skusme tu istu sekvenciu namapovat do genomu sliepky: stlacte najprv na hornej modrej liste Genomes, zvolte Vertebrates a Chicken a potom na hornej liste BLAT. Do okienka zadajte tu istu sekvenciu. Akú podobnosť a dĺžku má najsilnejší nájdený výskyt teraz? Na ktorom je chromozóme?
• Ako sa to porovna s hodnotami, ktore sme dostali pomocou BLASTu na NCBI?

Práca s tabuľkami, sťahovanie anotácií

• Položka Tables na hornej lište umožnuje robiť rafinované veci s tabuľkami, ktoré obsahujú súradnice génov a pod.
• Základná vec: vyexportovať napr. všetky gény v zobrazenom výseku v niektorom formáte:
  • sequence: fasta súbor proteínov, génov alebo mRNA s rôznymi nastaveniami
  • GTF: súradnice
  • Hyperlinks to genome browser: klikacia stránka
• Namiesto exportu si môžeme pozrieť rôzne štatistiky

• Zložitejšie: prienik dvoch tabuliek, napr. gény, ktoré sú viac než 50% pokryté simple repeats
  • V intersection zvolíme group: Variation and repeats, track: RepeatMasker, nastavíme records that have at least 50% overlap with RepeatMasker
  • V summary/statistics zistíme, kolko ich je v genóme, môžeme si ich preklikať cez Hyperlinks to genome browser

• Filter na tabuľku, napr. gény, ktoré majú v názve ribosomal (postup pre drozofilu):
  • V casti hg19.kgXref based filters políčko description dáme *ribosomal*

Dotploty

• Na stranke UCSC genome browseru http://genome.ucsc.edu/ si zadajte ludsky genom, najnovsiu verziu (hg19)
• V druhom okne/tabe si otvorte YASS server na adrese http://bioinfo.lifl.fr/yass/yass.php
• V genome browseri zadajte region chr21:10,018,860-10,019,178
  • tento región obsahuje Alu repeat, ktorý tvorí cca 10% ľudského genómu, viac ako milión kópií
  • zobrazte si DNA sekvenciu tohto useku takto: na hornej modrej liste zvolte DNA, na dalsej obrazovke tlacidlo get DNA
• DNA sekvenciu Alu opakovania chceme zarovnat samu k sebe programom YASS
  • DNA sekvenciu Alu opakovania skopirujte do okienka "Paste your sequences" v stranke Yass-u a dvakrat stlacte tlacidlo Select vedla okienka
  • Nizsie v casti "Selected DNA sequence(s)" by sa Vam malo v oboch riadkoch objavit "Pasted file 1"
  • Nizsie v casti "Parameters" zvolte "E-value threshold" 0.01 a stlacte "Run YASS"
  • Vo vysledkoch si pozrite Dotplot, co z neho viete usudit o podobnosti jednotlivych casti Alu opakovania?
  • Vo vysledkoch si pozrite Raw: blast, ake su suradnice opakujucej sa casti a kolko zarovnanie obsahuje zhod/nezhod/medzier? (Pozor, prve zarovnanie je cela sekvencia sama k sebe, druhe je asi to, co chcete)
• V genome browseri zadajte poziciu chr21:9,440,000-9,450,000 (10kb sekvencie na chromozome 21, s niekolkymi vyskytmi Alu)
• Chceme teraz porovnat tento usek genomu so sekvenciou Alu pomocou YASSu
  • Ako predtym si stiahnite DNA sekvenciu tohto useku
  • V YASSe chodte sipkou spat na formular
  • Skopirujte DNA sekvenciu do YASSoveho formulara, do okienka vpravo (vyznacit si ju mozete klavesovou kombinaciou Ctrl-A alebo Select All v menu Edit),
  • V casti formulara Selected DNA sequence(s) stlacte Remove pri hornom riadku
  • Pri pravom okienku, kam ste nakopirovali sekvenciu, stlacte Select
  • Zase stlacte Run YASS
  • Pozrite si vysledok ako Dotplot, kolko opakovani Alu ste nasli? Preco je jedno cervene?
  • Pozrite si Raw: blast, na kolko percent sa podoba najpodobnejsia a na kolko druha najpodobnejsia kopia?

Fylogeneticke stromy, mobyle portal

• V UCSC browseri mozeme ziskavat viacnasobne zarovnania jednotlivych genov (nukleotidy alebo proteiny). Nasledujuci postup nemusite robit, subor si stiahnite tu: http://compbio.fmph.uniba.sk/vyuka/mbi/poznamky/cb06/cb06-aln.fa
  • UCSC browseri si pozrieme usek ludskeho genomu chr6:136,214,527-136,558,402 s genom PDE7B (phosphodiesterase 7B)
  • Na modrej liste zvolime Tables, v nej RefSeq genes, zaklikneme Region: position, a Output fomat: CDS FASTA alignment a stlacime Get output
  • Na dalsej obrazovke zaklikneme show nucleotides. Z primatov zvolime chimp, rhesus, tarsier, z inych cicavcov mouse, rat, dog, elephant a z dalsich organizmov opposum, platypus, chicken, lizard, stlacime Get output.
  • Vystup ulozime do suboru, z mien sekvencii zmazeme spolocny prefix NM_018945_, pripadne celkovo prepiseme mena na anglicke nazvy

• Skusme zostavit strom na stranke http://mobyle.pasteur.fr/cgi-bin/portal.py
• Pouzijeme program quicktree, metodu neighbor joining, bootstrap 100
• Na zobrazenie stromu vysledok dalej prezenieme cez zobrazovacie programy drawtree alebo newicktops (zvolit v menu pri tlacidle further analysis)
  • Vysledok z drawtree, nezakoreneny, nezobrazuje bootstrap hodnoty
  • Vysledok z newicktops, zakoreneny na nahodnom mieste (nie spravne) zobrazuje bootstrap hodnoty
  • v drawtree sme nastavili sme formát výstupu MS-Windows Bitmap a X,Y resolution aspoň 1000, v newicktops sme nastavili show bootstrap values
• "Spravny strom" [17] v nastaveniach Conservation track-u v UCSC browseri (podla clanku Murphy WJ, Eizirik E, O'Brien SJ, Madsen O, Scally M, Douady CJ, Teeling E, Ryder OA, Stanhope MJ, de Jong WW, Springer MS. Resolution of the early placental mammal radiation using Bayesian phylogenetics. Science. 2001 Dec 14;294(5550):2348-51.)
• Nas strom ma long branch attraction (zle postavenie hlodavcov, ktori maju dlhu vetvu aj slona, co moze byt zapricene sekvenovacimi chybami).
• Ine programy, ktore mozete skusit na mobyle
  • phyml: metoda maximalnej vierohodnosti (daju sa nastavit detaily modelu, bootstraps, ktory ale moze dost dlho trvat, typy operacii na strome pri heuristickom hladani najlepsieho stromu)
  • dnapars alebo protpars na parsimony
  • viacnasobne zarovnanie pomocou clustalw alebo modernejsou alternativou muscle
  • Ak chcete skusat zarovnania, zacnite z nezarovnanych sekvencii: http://compbio.fmph.uniba.sk/vyuka/mbi/poznamky/cb06/cb06-seq.fa

Expresia génov, cvičenia pre biológov

Normalizacia dat

• (nepreberané na cvičeniach, uvádzané pre zaujímavosť)
• Pozri prednasku http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic/rigat/st416/

Nadreprezentacia, obohatenie (enrichment)

• Mnohe celogenomove analyzy nam daju zoznam genov, ktore sa v nejakom ukazovateli vyrazne lisia od priemeru.
• Napriklad geny s pozitivnym vyberom v komparativnej genomike, geny vyrazne nadexprimovane alebo podexprimovane v microrarray experimentoch, geny regulovane urcitym transkripcnym faktorom a pod.
• Niektore z nich budu preskumanejsie (znama funkcia a pod.), niektore mozu mat nejake udaje o funkcii prenesene z homologov a dalsie mozu byt uplne nezname
• Co s takym zoznamom "zaujimavych genov"?
• moznost 1: vybrat si z neho niekolko malo zaujimavych kandidatov a preskumat ich podrobnejsie (experimentalne alebo informaticky)
• moznost 2: zistit, ci tato cela skupina je obohatena o geny urcitych skupin
  • napr. v pripade pozitivneho vyberu nam casto vychadzaju geny suvisiace s imunitou, lebo su pod velkym evolucnym tlakom od patogenov
  • takato analyza nam teda da informaciu o suvislostiach medzi roznymi procesmi
• Priklad (Kosiol et al)
  • 16529 genov celkovo, 70 genov v GO kategorii innate immune response (0.4% zo vsetkych genov)
  • 400 genov s pozivnym vyberom, mame 8 genov s innate immune response (2% zo vsetky genov s poz. vyb.)
• Celkovy pocet genov n, imunitnych ni, pozitivny vyber np, imunitnych s poz. vyb. nip.
• Nulova hypoteza: geny v nasom zozname boli nahodne vybrane z celeho genomu, t.j. ak v celom genome je frekvencia imunitnych genov ni/n, vo vzorke velkosti np (geny s pozitivnym vyberom) ocakavame cca np * ni / n imunitnych genov.
  • v nasom priklade by sme ocakavali 1.7 genu s innate immune response, ale mame 8 (4.7xviac)
• Aj v nulovej hypoteze vsak vzorka velkosti ni cisto nahodou moze obsahovat viac alebo menej takych genov.
• Rozdelenie pravdepodbnosti nip je hypergeometricke, t.j. pravdepodobnost ze nip=x je ${n_{i} \choose x}{n-n_{i} \choose n_{i}-x}/{n \choose n_{p}}$
• Aka je pravdepodobnost, ze v nulovej hypoteze bude nip tolko, kolko sme namerali alebo viac? (Chvost rozdelenia). V nasom pripade p-value 2.8e-4.
• Hypergeometric or Fisher's exact test, pripadne ich aproximacie pre velke hodnoty v tabulke (chi^2 test) zisti, ci sa nasa tabulka velmi lisi od toho, co by sme ocakavali v nulovej hypoteze
• Ak skusame vela (vsetky) GO kategorie, musime uplatnit "multiple testing correction": napr. skusame 1000 kategorii, ocakavali by sme cisto nahodou 50 kategorii, ktore budu vyznamne pre p<=0.05. Ak najdeme takych 100, tak cca polovica z nasich vysledkov su falosne, cisto nahodne asociacie. Takze chceme znizit p-value tak, aby ocakavane nahodne asociacie netvorili vyznamnu cast vysledkov.

• Suvisiace clanky
  • Rivals I, Personnaz L, Taing L, Potier MC (February 2007). "Enrichment or depletion of a GO category within a class of genes: which test?". Bioinformatics (Oxford, England) 23 (4): 401–7. doi:10.1093/bioinformatics/btl633. PMID 17182697.
  • Huang da W, Sherman BT, Lempicki RA (January 2009). "Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists". Nucleic Acids Research 37 (1): 1–13. doi:10.1093/nar/gkn923. PMID 19033363.

• Existuju web servery, napr. GOrilla pre ludske geny: http://cbl-gorilla.cs.technion.ac.il/, DAVID (http://david.niaid.nih.gov), g:Profiler http://biit.cs.ut.ee/gprofiler/
• Treba dat pozor, ci pocitaju to co chceme

• Kod v statistickom systeme R na pocitanie hypergeometrickeho rozdelenia

> dhyper(0:70, 70, 16529-70, 400);
 [1]  1.793421e-01  3.126761e-01  2.679872e-01  1.505169e-01  6.231088e-02
 [6]  2.027586e-02  5.400796e-03  1.210955e-03  2.332580e-04  3.920215e-05
[11]  5.818723e-06  7.702558e-07  9.166688e-08  9.873221e-09  9.678760e-10
[16]  8.677204e-11  7.143849e-12  5.420388e-13  3.802134e-14  2.472342e-15
[21]  1.493876e-16  8.405488e-18  4.412274e-19  2.164351e-20  9.935473e-22
[26]  4.273662e-23  1.724446e-24  6.533742e-26  2.326517e-27  7.791092e-29
[31]  2.455307e-30  7.285339e-32  2.036140e-33  5.361856e-35  1.330660e-36
[36]  3.112566e-38  6.862558e-40  1.426089e-41  2.792792e-43  5.153006e-45
[41]  8.955105e-47  1.465159e-48  2.255667e-50  3.265636e-52  4.442631e-54
[46]  5.674366e-56  6.797781e-58  7.629501e-60  8.012033e-62  7.860866e-64
[51]  7.193798e-66  6.129013e-68  4.851139e-70  3.558526e-72  2.412561e-74
[56]  1.506983e-76  8.641725e-79  4.530590e-81  2.161126e-83  9.326620e-86
[61]  3.617279e-88  1.250737e-90  3.817900e-93  1.016417e-95  2.323667e-98
[66] 4.469699e-101 7.034762e-104 8.698702e-107 7.924236e-110 4.728201e-113
[71] 1.386176e-116
phyper(7, 70, 16529-70, 400, lower.tail=FALSE);

Prakticke cvicenie

Data o expresii ludskych genov v roznych tkanivach a podobne v UCSC genome browseri

• Chodte na stranku http://genome.ucsc.edu/, najdite PTPRZ1 gen v ludskom genome
• Zvolte Gene Sort, sort by nechajme Expression (GNF Atlas 2), v configure nechajme iba name, search PTPRZ1, output text
• http://biit.cs.ut.ee/gprofiler/ mena genov skopirujme do policka Query, stlacte g:Profile!
• Co by sme na zaklade nadreprezentovanych kategorii usudzovali o tomto gene?
• Najdite tento gen v Uniprote (http://www.uniprot.org/), potvrdzuje nase domnienky?

• Vratme sa do do genome browsera, najdime si PTPRZ1 gen v genome
• V browseri su rozne tracky tykajuce sa expresie, napr. GNF Atlas 2. Precitajte si, co je v tomto tracku zobrazene, zapnite si ho a pozrite si expresiu okolitych genov okolo PTPRZ1
• Kliknite na gen v tracku UCSC known genes. V tabulke uvodite zase prehlad expresie v roznych tkanivach (podla GNF Atlasu), linku na Visigene.

NCBI Gene Expression Omnibus http://www.ncbi.nlm.nih.gov/geo/

• Databaza gene expression dat na NCBI
• Do okienka Data sets zadajme GDS2925
• Mali by sme dostat (Various weak organic acids effect on anaerobic yeast chemostat cultures)
• Mozeme si pozriet zakladne udaje, napr. citation, platform
• Link expression profiles nam zobrazi grafy pre rozne geny
• Pri kazdom profile mozeme kliknut na profile neighbors, aby sme videli geny s podobnym profilom
• Data analysis tools, cast Cluster heatmaps, K-means, skuste rozne pocty clustrov

Proteíny 2, cvičenia pre informatikov

• Viac informacii: [18], Bafna, Reinert 2004

Uvod do proteomiky

Gélová elektroforéza (gel electrophoresis)

• Izolovanie jednotlivých proteínov, porovnávanie ich množstva.
• Negatívne nabité proteíny migrujú v géli v elektrickom poli. Väčšie proteíny migrujú pomalšie, dochádza v oddeleniu do pruhov. Táto metóda sa používa aj na DNA a RNA. Pre proteíny možno tiež robiť 2D gél (podľa hmotnosti a náboja).
• Bioinformatický problém: zisti, ktoré fliačiky na dvoch 2D géloch zodpovedajú tým istým proteínom.
• Automatizovanejšia technológia: kvapalinová chromatografia (liquid chromatography) - separácia proteínov v tenkom stĺpci

Hmotnostná spektrometria (mass spectrometry)

• Hmotnostná spektrometria meria pomer hmostnosť/náboj molekúl vo vzorke.
• Používa sa na identifikáciu proteínov, napr. z 2D gélu.
• Proteín nasekáme enzýmom trypsín (seká na [KR]{P}) na peptidy
• Meriame hmostnosť kúskov, porovnáme s databázou proteínov.
• Tandemová hmotnostná spektrometria (MS/MS) ďalej fragmentuje každý kúsok a dosiahne podrobnejšie spektrum, ktoré obsahuje viac informácie

Sekvenovanie proteinov pomocou MS/MS

Berme do uvahy len y-iony, ktore zodpovedaju hmotnosti prefixu

• Vsetky hmotnosti budu cele cisla, pozname celkovu hmotnost peptidu M a hmotnosti aminokyselin
• Mame skorovaciu funkciu f(m), ktora pre prefix peptidu s hmotnostou m urci skore, ako dobre sedi so spektrom
• Chceme najst peptid, ktory maximalizuje sucet skore svojich prefixov
• Dynamicke programovanie s podproblemom A[m] je skore najlepsieho prefixu s hmotnostou m
  • Rekurencia? Zlozitost? Je to polynomialny algoritmus? (Co je vlastne vstup?)

Co ak berieme do uvahy aj b-iony, ktore meraju hmotnost sufixu?

• Idea 1: upravme skorovaciu funkciou f(m): skontroluje signal na pozicii m aj M-m
• Problem: jeden signal sa moze ratat dvakrat, raz ako b-ion, raz ako y-ion, algoritmus ma tendenciu pridavat taketo artefakty
• Nove skorovanie: ak hmotnost nejakeho prefixu a nejakeho sufixu su rovnake, zarataj ich skore iba raz (skore peptidu je skore mnoziny hmotnosti jeho prefixov a sufixov)
• Definuj novy podproblem: A[p,s] je najlepsie skore, ktore moze dosiahnut prefix s hmotnostou p a sufix s hmotnostou s
• Pre p+s<=M a p<=s mame A[p,s] = max_a A[p, s-a] + f(s) * (1-delta(p,s)) a naopak pre p>s mame A[p,s] = max_a A[p-a, s] + f(p)
• Ako ukoncime dynamicke programovanie? Zlozitost?
• Zrychlenie: staci uvazovat s od p-w po p+w kde w je maximalna hmotnost aminokyseliny

Detekcia znamych proteinov pomocou MS

• Predikcia spektra pre dany peptid, porovnanie s realnym spektrom, zlozite skorovacie schemy
• Filtrovanie kandidatov na proteiny, ktore obsahuju peptidy s pozorovanou hmotnostou
• Problem: mame danu databazu proteinov a cielovu hmotnost peptidu M, pozname hmotnost kazdej aminokyseliny. Najdite vsetky podretazce s hmotnostou M.
• Databazu proteinov si vieme predstavit aj ako postupnost cisel - hmotnosti aminokyselin, hladame intervaly su suctom M.
• Trivialny algoritmus: zacni na kazdej pozicii, pricitavaj kym nedosiahnes hmotnost>=M. Zlozitost? Vieme zlepsit?
• Predspracovanie: pocitajme hmotnosti vsetkych podretazcov, potom vyhladajme binarne. Zlozitost?
  • Dat sa aj cas na predspracovanie pomocou FFT (Fast Fourier Transform) Bansal, Cieliebak, Liptak 2004

Hľadanie motívov, cvičenia pre informatikov

Hladanie motivov zadefinovanych pravdepodobnostnou maticou

• Mame danych n sekvencii S=(S1..Sn), kazda dlzky m, dlzku motivu L, nulova hypoteza q (rozdelenie nukleotidov)
• Hladame motiv vo forme pravdepodobnostneho profilu dlzky L a jeho vyskyt v kazdej sekvencii
• Nech $w_{{i,a}}$ je pravdepodobnost, ze na pozicii i motivu bude baza a, W cela matica
• o_i je pozicia vyskytu v sekvencii S_i, O su vsetky vyskyty dokopy
• Pr(S|W.O) je jednoduchy sucin, kde pre pozicie v oknach pouzijeme pravdepodobnosti z W, pre pozicie mimo okna pouzijeme q
• Hladame W,O, ktore maximalizuju tuto pravdepodobnost (najvierhododnejsie)
  • Nepozname efektivny algoritmus, ktory by to 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
  • EM algoritmus: iterativne zlepsuje W, pricom berie Vsetky O vahovane podla ich pravdepodohnosti vzhladom na W z minuleho kola
• Casto sa v zlozitych pravdepdobnostnych modeloch pouziva namiesto optimalizacie parametrov vzorkovanie (sampling)
  • Rozsirime model, aby aj O a W boli nahodne premenne
  • Vzorkujeme z Pr(W,O|S)
  • Mozeme pre kazdu vzorku spocitat Pr(S|W,O) a zobrat najlepsiu
  • Alebo mozeme robit priemery roznych velicin, skumat tvar rozdelenia
  • My pouzijeme metodu Markov chain Monte Carlo (MCMC) a konkretne jej variant Gibbsovo vzorkovanie (Gibbs sampling)

Hladanie motivov pomocou Gibbsovho vzorkovania

Zadefinujeme cele ako velky pravdepodobnostny model:

• 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 nageneruje sekvencia podla profilu w
• mimo okna sa nageneruje sekvencia z nulovej hypotezy (nezavisle bazy)

Gibbsovo vzorkovanie:

• Mame dane S, vzorkujeme O ($O^{{(0)}},O^{{(1)}},\dots$) (z nich uz lahko zostavime maticu W)
  • 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^{{(t}})_{{-i}},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 poyicie v sekvencii i za prave zvolenu
  • opakuj vela krat
• Konverguje k cielovemu rozdeleniu $\Pr(O|S)$, ale vzorky nie su nezavisle
  • mozeme si pamatat, ktore okno ako casto vnutri motivu alebo konfiguraciu s najvacsou pravdepodobnostou, aku sme stretli pocas vzorkovania
• 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)
• 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_{{i,a}}\geq 0$ a $\sum _{a}w_{{i,a}}=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)

• Celkovy trik: vieme lahko spocitat pravdepodobnost dat ak vieme vsetky parametre $\Pr(S|W,O)$
• Pomocou integralu vieme tiez odmarginalizovat W, t.j. spocitat $\Pr(S|O)$
• Skusame vsetky mozne hodnoty $o'_{i}$, pocitame pravdepodobnost dat, vzorkujeme umerne k tomu

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.

Gibbs sampling, Gibbsovo vzorkovanie vseobecne

• Cielove rozdelenie ma n premennych $\pi (x_{1},...x_{n})$
• V kazdom kroku vzorkujeme jednu premennu 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
• Premennu $x_{i}$ zvolime nahodne alebo periodicky $i=1,2,\dots ,n$
• Vzorky nie su nezavisle, no vieme dokazat nieco o konvergencii k $\pi$ (pozri nizsie)

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 je $X_{t}$ z konecnej abecedy (mozne hodnoty $X_{t}$ nazyvame stavy)
• Napriklad stavy A,C,G,T
• V nasom Gibbsovom vzorkovani bude 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)
• Pravdepodobnosti prechodu medzi stavmi za jeden krok mozeme vyjadrit maticou pravdepodobnosti P, ktorej prvok $p_{{i,j}}$ oznacuje pravdepodobnost prechodu zo stavu i do stavu j $p_{{i,j}}=\Pr(X_{t}=j|X_{{t-1}}=i)$
• Sucet kazdeho riadku je 1, cisla nezaporne
• Ako $p_{{ij}}^{{(t)}}$ budeme oznacovat $\Pr(X_{t}=i|X_{{0}}=j)$, tieto hodnoty dostaneme umocnenim matice P na t
• Vravime ze matica je ergodicka (niekedy volane regularna), 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

• 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$)
• Pre ergodicke matice existuje prave jedno stacionarne rozdelenie, pre ktore navyse plati, ze pre kazde i a j plati $\lim _{{t\to \infty }}p_{{i,j}}^{{(t)}}=\pi (j)\,$
• V HMM stavy tvoria Markovov retazec; hladanie genov ergodicky stavovy priestor, profilove HMM nie
• 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

Markov chain Monte Carlo MCMC

• Chceme generovať náhodné vzorky z nejakeho cieloveho rozdelenia $\pi$, ale toto rozdelenie 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, 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)]$

Dokaz spravnosti Gibbsovo 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$
  • Dokazeme tzv. detailed balance, ze pre kazde dva vektory hodnot x a x' mame $\pi (x)p_{{x,x'}}=\pi (x')p_{{x',x}}$
  • ak pre nejaky retazec a nejaku distribuciu $\pi$ plati detailed balance, $\pi$ je stacionarna distribucia: $\sum _{i}\pi (i)p_{{i,j}}=\sum _{i}\pi (j)p_{{j,i}}=\pi (j)\sum _{i}p_{{j,i}}=\pi (j)\,$
  • Pre nas retazec plati detailed balance: uvazujme dva za sebou iduce vektory hodnot x a x', ktore sa lisia v i-tej suradnici. Nech $x_{{-i}}$ su hodnoty vsetkych ostatnych premennych okrem $x_{i}$
  • $\pi (x)p_{{x,x'}}=\pi (x)\Pr(x'_{i}|x_{{-i}})=\Pr(x_{{-i}})\Pr(x_{i}|x_{{-i}})\Pr(x'_{i}|x_{{-i}})=\pi (x')\Pr(x_{i}|x_{{-i}})=\pi (x')\Pr(x_{i}|x'_{{-i}})=\pi (x')p_{{x',x}}$

Hľadanie motívov a populačná genetika, cvičenia pre biológov

Zhrnutie semestra

Tvorba bioinformatickeho nastroja:

• sformulujeme biologicke ciele (ake mame data, ake typy otazok sa chceme pytat)
• sformulujeme informaticky/matematicky, napr. ako pravdepodobnostny model, dostaneme informaticke zadanie problemu, v ktorom je presne dany vztah medzi vstupom a zelanym vystupom (napr. najst zarovnaie s maximalnym skore v urcitej skorovacej scheme)
• hladame efektivne (rychle) algoritmy na riesenie informatickeho problemu
• ak sa nam nepodari najst dost rychly algoritmus, pouzijeme heuristiky, ktore davaju priblizne riesenia

• Konkretne nastroje a webstranky sa rychlo menia, celkove principy sa menia pomalsie
• Pre spravne interpretovanie vysledkov je dolezite poznat model, predpoklady, ktore autori nastroja pouzili, resp. zdroj dat v prislusnej databaze

Prehlad preberanych tem:

• Zarovnania (skorovanie ako pravdepodobnostny model, dynamicke programovanie, heuristicke zarovnanvanie, E-value, lokalne/globalne, parove vs. viacnasobne, celogenomove)
• Hladanie genov (skryte Markovove modely)
• Evolucia (pravdepodobnostne modely substitucii, metoda maximalnej vierohodnosti, metoda maximalnej uspornosti, metoda spajania susedov)
• Komparativna genomika (hladanie konzervovanych oblasti, komparativne hladanie genov, pozitivny vyber, fylogeneticke HMM, kodonove matice)
• Proteiny (predikcia struktury, profily a profilove HMM rodin/domen, protein threading)
• Expresia génov (zhlukovanie, klasifikácia, regulačné siete)
• Transkripčné faktory (hladanie vyskytov znamych motivov, hladanie novych motivov, regularne vyrazy, profily)
• Populačná genetika (geneticky drift, mapovanie asociacii, linkage disquilibrium, struktura populacie)
• RNA struktura (dynamicke programovanie, stochasticke bezkontextove gramatiky)

Hľadanie motívov

Opakovanie z prednasky

• Hladame vazobne miesta TF
• Reprezentujeme ako sekvencne motivy (vo forme konsenzus retazcov, regularnych vyrazov alebo PSSM)
• Dva bioinformaticke problemy:
• Hladanie novych motivov
  • hladame motiv, ktory sa opakuje vo vstupnych sekvenciach
  • sekvencie ziskame pomocou ChIP, alebo ako promotery koregulovanych genov
  • vypoctovo narocny problem
  • nasli sme to, co sme chceli?
• Hladanie novych vyskytov znameho motivu
  • vypoctovo jednoduchy problem
  • zvycajne najdeme vela vyskytov, nie vsetky biologicky relevantne

Programy na pracu s motivmi: Mobyle

• Nejake programy na pracu s motivmi najdete na stranke http://mobyle.pasteur.fr/
• Hladanie vyskytov znamych motivov
  • iANTPatScan
  • pftools, sig
  • scan_for_matches
  • prophet, profit (matice)
  • jaspscan (jaspar db), tfscan (old version of transfac)
• Hladanie novych motivov
  • SMILE
  • CONSENSUS
  • prophecy (matrix from alignment)
• Programy roznej kvality, treba skontrolovat, ci robia to, co chcete

MEME

• Znamy program na hladanie novych motivov, webserver MEME http://meme.nbcr.net/ (na cviceni pro pocitaci)

Komparativne hladanie motivov

• Z vela vyskytov motivu chceme najst tie skutocne funkcne
• Purifikacny vyber: funkcne sekvencie sa menia pomalsie ako zvysok genomu
• Toto by malo platit aj pre vazobne miesta TF
• Chceme teda najst vyskyty motivu zachovane aj v pribuznych organizmoch
• Mala fylogeneticka vzdialenost (napr. clovek-simpanz): mala sanca, ze bude mutacia v kratkom okne, nie je mozne statisticky rozpoznat "zachovane" sekvencie pre kratke motivy
• Velka fylogeneticka vzdialenost: nekodujuce oblasti tazko zarovnatelne, takze aj homologicke vyskyty motivu nebudu zarovnane (lebo su kratke)
• Ako zostavujeme celogenomove zarovnania:
  • najdeme signifikantne lokalne zarovnania (mala e-value)
  • tie spajame dokopy do retazi a sieti
  • ak je jeden kratky motiv uprostred neutralne sa vyvijajucej DNA, nebude mat signifikantne lokalne zarovnanie
• Mnohe programy dovoluju posun motivu v urcitom okne zarovnania
• Priklad: Kheradpour P, Stark A, Roy S, Kellis M (December 2007). "Reliable prediction of regulator targets using 12 Drosophila genomes". Genome Res. 17 (12): 1919–31. doi:10.1101/gr.7090407. PMID 17989251.
  • Drosophila melanogaster ako referencny genom
  • pre kazdy vyskyt motivu hladaj v ostatnych genomoch v urcitom okne
  • spocitaj v akom percente stromu najdene vyskyty
  • spocitaj p-value pomocou motivov s nahodne poprehadzovanymi stlpcami

RNA struktura

• Buduca prednaska: ako predpovedat strukturu RNA molekul, resp. hladat RNA geny

Prakticke cvicenie

MEME

• Chodte na stranku http://meme.nbcr.net/
• Zvolte nastroj MEME a do okienka "actual sequences" zadajte tieto sekvencie
• Pozrite si ostatne nastavenia. Co asi robia?
• Ak server pocita dlho, mozete si pozriet vysledky tu

RNA

• Znama databaza rodin RNA genov je Rfam
• Najdite si v nej rodinu RF00015 (U4 spliceosomal RNA)
• V casti secondary structure si mozete pozriet obrazky farebne kodovane podla roznych kriterii
• Jedna z mnohych ludskych kopii je tato:

AGCTTTGCGCAGTGGCAGTATCGTAGCCAATGAGGTTTATCCGAGGCGCG
ATTATTGCTAATTGAAAACTTTTCCCAATACCCCGCCATGACGACTTGAA
ATATAGTCGGCATTGGCAATTTTTGACAGTCTCTACGGAGA

• Skuste ju najst v ludskom genome nastrojom BLAT v UCSC genome browseri
• Pozrite si tracky Mapability, RNA genes, Mammalian conservation, Repeats v jej okoli
• Pozrite si track "CSHL Sm RNA-seq" ktory obsahuje RNASeq kratky RNA z roznych casti buniek, zapnite si v jeho nastaveniach aj zobrazenie RNA z jadra (nucleus)
• Zadajte sekvenciu do na RNAfold serveri [19]
• Ak vypocet dlho trva, pozrite si vysledok tu
• Podoba sa na strukturu zobrazenu v Rfame? v com sa lisi?

Populacna genomika v UCSC genome browseri

• Verzia hg18 na UCSC genome browseri ma viacero trackov tykajucich sa populacnej genomiky a polymorfizmov
• V casti Phenotype and Disease Associations si zapnite GAD view
• V casti Variation and Repeats si zapnite
  • HGDP Allele Freq na Pack (po kliknuti na SNP zobrazi mapu sveta s distribuciou alel)
  • "HapMap LD Phased" na Full
  • "DGV Struct Var" na Pack
• Pozrite si napriklad region chr2:164,862-426,468
• Track Genome Variants obsahuje genomy niekolkych ludi, napr Jim Watsona

• Browser diverzity u S.cerevisae: [20]

Populačná genetika, cvičenia pre informatikov

Určovanie štruktúry populácie

• Pritchard JK, Stephens M, Donnelly P (June 2000). "Inference of population structure using multilocus genotype data". Genetics 155 (2): 945–59. PMID 10835412.

• Majme N diploidnych jedincov, L genotypovanych SNPov, SNPy navzajom nezavisle (v stave LE), pocet populacii K
• X[i,l,j] - gentotyp jedinca i v SNPe l pre j=1,2 (kazdy jedinec dve alely lubovolne ocislovane 1,2) (zvacsa binarna premenna)
• Z[i,l,j] - z ktorej subpopulacie pochadza alela X[i,l,j] (cislo z {1...k})
• Q[i,k] - aka cast genomu jedinca i pochadza z populacie k (realne cislo)
• P[k,l,a] - frekvencia alely a v SNPe l v populacii k (realne cislo)

• P pochadza z nejakeho apriorneho rozdelenia, napr. rovnomerne rozdelenie, Dirichletovo rozdelenie, nezavisle pre kazde k,l
• Podobne Q (nezavisle pre kazde i)
• Pr(Z[i,j,l] = k|P,Q) = Q[i,k] a Pr(Z|P,Q) je sucin takychto clenov
• Pr(X[i,l,j]=a|Z,P,Q) = P[Z[i,l,j],l,a] a Pr(X|Z,P,Q) je sucin takychto clenov (SNPy nezavisle)
• To nam urcuje Pr(P,Q,Z,X) = Pr(P)Pr(Q)Pr(Z|P,Q)Pr(X|Z,P,Q)
• My chceme Pr(Q|X)

Algoritmus Gibbsovho vzorkovania:

• Zvol pociatocne $Z^{{(0)}}$
• Opakuj:
  • Zvol nahodne $P^{{(m)}},Q^{{(m)}}$ z $Pr(P,Q|X,Z^{{(m-1)}})$
  • Zvol nahodne $Z^{{(m)}}$ z $Pr(Z|X,P^{{(m)}},Q^{{(m)}})$
• Vzorce v clanku

• Mierna komplikacia: ak aproximujeme $E[Q[i,k]|X]$ pomocou priemeru $Q^{{(m)}}$, mali by sme dostat 1/K kvoli symetrii (K! symetrickych rieseni)
• Nastastie sa Gibbsovo vzorkovanie malokedy presuva medzi roznymi oznackovaniami tych istych populacii
• Inak musime pouzit nejako preznacit populacie vo vysledku aby boli ekvivalentne v roznych vzorkach

Zhrnutie semestra

Precvičenie dynamického programovania

• zarovnavanie sekvencii, viacero variantov (globalne, lokalne, afinne medzery)
• skryte Markovove modely (Viterbiho algoritmus, dopredny a spatny algoritmus)
• vypocty na stromoch (uspornost, vierhodnost - Felsensteinov algoritmus)
• hmotnostna spektrometria (MCMC)
• sekundarna struktura RNA

Pravdepodobnostne modely

• skryte Markovove modely (hladanie genov, konzervovanych oblasti, fylogeneticke a parove HMM, profilove HMM, transmembranove proteiny)
• fylogeneticke stromy a substitucne modely
• Markovove retazce
• Gibbsovo vzorkovanie/MCMC
• metoda maximalnej vierohodnosti

Ine

• integer linear programming
• kombinatorika na stromoch
• deBruijnove grafy
• statisticke testy (LD, pozitivny vyber)

Ako modelovat problemy realneho sveta

• rozmysliet si, ake mame data, co by sme chceli ako vysledok
• sformulovat ako informaticky problem (napr. optimalizacia nejakeho skore)
• pravdepodobnostne modely nam casto dovolia zvolit skorovaciu schemu systematickym sposobom
• vysledny problem casto NP tazky
  • heuristiky, aproximacne algoritmy
  • ILP a ine techniky na presne riesenie
  • neda sa problem trochu preformulovat?
• testovanie: su vypoctove vysledky relevantne v danej domene? (bola formulacia dostatocne realisticka?)

Články na journal club

Svoje preferencie ohľadom článkov na journal club odovzdajte pomocou tejto linky najneskôr do 20.10.2010. Vo forme uveďte čísla troch článkov z tohto zoznamu v poradí podľa preferencie.

1. Lunter G (2007). "Probabilistic whole-genome alignments reveal high indel rates in the human and mouse genomes.". Bioinformatics 23 (13): i289-96. doi:10.1093/bioinformatics/btm185. PMID 17646308. Čítajú Švec, Ďudák, Jurík, Slavkovský, Vlčák, Porubský
2. Wapinski I, Pfeffer A, Friedman N, Regev A (2007). "Natural history and evolutionary principles of gene duplication in fungi". Nature 449 (7158): 54-61. doi:10.1038/nature06107. PMID 17805289.
3. Clamp M, Fry B, Kamal M, Xie X, Cuff J, Lin MF et al. (2007). "Distinguishing protein-coding and noncoding genes in the human genome.". Proc Natl Acad Sci U S A 104 (49): 19428-33. doi:10.1073/pnas.0709013104. PMID 18040051.
4. Harbison CT, Gordon DB, Lee TI, Rinaldi NJ, Macisaac KD, Danford TW et al. (2004). "Transcriptional regulatory code of a eukaryotic genome.". Nature 431 (7004): 99-104. doi:10.1038/nature02800. PMID 15343339. Čítajú Goffa, Balážfyová, Trojáková, Bača, Budiš
5. Bystroff C, Thorsson V, Baker D (August 2000). "HMMSTR: a hidden Markov model for local sequence-structure correlations in proteins". J. Mol. Biol. 301 (1): 173–90. doi:10.1006/jmbi.2000.3837. PMID 10926500.
6. Andronescu M, Fejes AP, Hutter F, Hoos HH, Condon A (2004). "A new algorithm for RNA secondary structure design.". J Mol Biol 336 (3): 607-24. doi:10.1016/j.jmb.2003.12.041. PMID 15095976.

Ukážkové príklady na skúšku

Aspoň polovicu bodov na skúške bude možné pre biológov aj informatikov možné získať z príkladov typu uvedeného nižšie. Samozrejme, na skúške nepoužijeme všetky tieto príklady a konkrétne reťazce, čísla, stromy a pod. budú iné. Na skúške máte dovolené používať kalkulačku s bežnými matematickými funkciami (nie však zložitejšie výpočtové zariadenia) a ťahák na dvoch listoch formátu A4, ktoré môžu byť z oboch strán popísané alebo potlačené ľubovoľným obsahom. Ostatné príklady na skúške budú prekvapením.

• Nájdite najkratšie spoločné nadslovo reťazcov GACAATAA, ATAACAC, GTATA, TAATTGTA.
• Zostavte deBruijnov graf stupňa k=3 pre reťazce GACAATAA, ATAACAC, GTATA, TAATTGTA a zistite, či má Eulerov ťah. Ak áno, akému reťazcu zodpovedá? Ak nie, prečo?
• Vyplňte maticu dynamického programovania pre lokálne (resp. globálne) zarovnanie reťazcov TACGT a CAGGATT, pričom zhodu skórujeme ako +3, nezhodu -1, medzeru -2. Napíšte aj optimálne zarovnanie, ktoré ste takto našli.
• Spočítaje skóre nižšieuvedeného zarovnania, pričom použijete skórovaciu maticu uvedenú nižšie, začatie medzery -5, rozšírenie medzery o jednu ďalšiu bázu -2. Nájdite globálne zarovnanie s vyšším skóre pre tieto dve sekvencie (netreba nájsť optimálne zarovnanie; pri hľadaní môžete použiť ľubovoľný postup alebo úvahu) a spočítajte aj skóre vášho nového zarovnania.

Zarovnanie:                             Matica:
ATAGTTTAA                                 A   C   G   T
A-GGG--AA                             A   2  -2  -1  -2
                                      C  -2   1  -2  -1    
                                      G  -1  -2   1  -2
                                      T  -2  -1  -2   2

• Uvažujme BLASTn, ktorý začína z jadier veľkosti w=3. Koľko jadier nájde pri porovnávaní sekvencií GATTACGGAT a CAGGATT? Ktoré to budú?
• Pre model na strane 16 prednášky o hľadaní génov (bol by v zadaní) spočítajte pravdepodobnosť vygenerovania báz AGT a stavov modrý,červený,modrý.
• Na strome na strane 6 v prednáške o evolúcii (bol by v zadaní) nájdite najúspornejšie ancestrálne znaky pre stĺpec zarovnania TTAAA (v poradí glum, hobit, človek, elf, ork). Nemusíte použiť algoritmy z prednášky resp. cvičení.
• Nájdite najúspornejší strom pre zarovnanie uvedené nižšie. Aká je jeho cena (koľko mutácií je nutných na vysvetlenie týchto sekvencií)? Odpoveď môžete spočítať ľubovoľným spôsobom.

vtáčik biely       ACAACGTCT
vtáčik čierny      TCTGAATCA
vtáčik sivý        TGTGAAAGA
vtáčik modrý       ACTACGTCT
vtáčik zelený      TGTGAAAGA

• Uvažujme maticu vzdialeností uvedenú nižšie. Ktorú dvojicu vrcholov spojí metóda spájania susedov ako prvú a aká bude nová matica po spojení?

                biely   čierny  sivý   modrý
vtáčik biely      0       5       7       4
vtáčik čierny     5       0       8       5
vtáčik sivý       7       8       0       5
vtáčik modrý      4       5       5       0

• Uvažujme strom ako na strane 6 v prednáške o evolúcii (bol by v zadaní), pričom každá hrana má rovnakú dĺžku a pravdepodobnosť každej mutácie na jednej hrane je 0.1 (t.j. napr. Pr(C|A,t)=0.1) a teda pravdepodobnosť zachovania tej istej bázy je 0.7, pravdepodobnosť každej bázy v koreni je 0.25. Aká je pravdepodobnosť, že v listoch dostaneme TTAAA a vo vnútorných vrcholoch samé Áčka? Nájdite priradenie ancestrálnych báz vo vnútorných vrcholoch, ktoré má väčšiu pravdepodobnosť a spočítajte, aká tá pravdepodobnosť je (nemusíte nájsť najlepšie možné priradenie).
• Zostavte profil (PSSM) pre zarovnanie sekvencií uvedené nižšie, pričom predpokladáme, že v celej databáze A tvorí 60% a T 40% všetkých sekvencií (iné bázy neuvažujeme). Použite prirodzený logaritmus (ln) a nepoužívajte pseudocounty.

AATA
TATA
TAAA
TTAT
TTAA

• Uvažujme microarray experimenty pre 5 génov. Medzi každými dvomi profilmi sme spočítali vzdialenosť (napr. pomocou Pearsonovho korelačného koeficientu) a dostali sme tabuľku vzdialeností uvedenú nižšie. Nájdite hierarchické zhlukovanie týchto génov, pričom vzdialenosť medzi dvoma zhlukmi (clustrami) bude vzdialenosť najbližších génov v nich.

         A    B    C    D    E
gén A    0   0.6  0.1  0.3  0.7    
gén B   0.6   0   0.5  0.5  0.4
gén C   0.1  0.5   0   0.6  0.6
gén D   0.3  0.5  0.6   0   0.8
gén E   0.7  0.4  0.6  0.8   0

• Uvažujte motív reprezentovaný profilom (skórovaciou maticou, PSSM) uvedenou nižšie. Spočítajte skóre reťazca GGAG. Ktorá sekvencia dĺžky 4 bude mať najmenšie a ktorá najväčšie skóre?

A   -3    3   -2   -2
C   -2   -2    1   -2
G    0   -2   -1    3
T    1   -1    1   -2

• Nájdite všetky výskyty regulárneho výrazu TA[CG][AT]AT v sekvencii GACGATATAGTATGTACAATATGC.
• Doplňte chýbajúce hodnoty za otázniky v matici dynamického programovania (Nussinovov algoritmus) pre nájdenie najväčšieho počtu dobre uzátvorkovaných spárovaných báz v RNA sekvencii GAACUAUCUGA (dovoľujeme len komplementárne páry A-U, C-G) a nakreslite sekundárnu štruktúru, ktorú algoritmus našiel.

 0 0 0 1 1 2 2 3 3 ? ?
   0 0 0 1 1 2 2 3 3 ?
     0 0 1 1 2 2 2 3 3
       0 0 1 1 1 1 2 3
         0 1 1 ? 1 2 3
           0 1 1 1 2 2
             0 0 0 1 2
               0 0 1 1
                 0 0 1
                   0 0
                     0

• Pre dvojice SNPov, ktorých tabuľky sú uvedené nižšie, určite, či môžeme štatisticky vylúčiť hypotézu, že sú v stave LE (linkage equilibrium) pri hladine významnosti p=0.05, resp. $\chi ^{2}>3.841$. Pre každú dvojicu spočítajte veličinu $\chi ^{2}$.

    Q   q              Q  q             Q  q
P  100 200          P 10  20         P  1  2
p  300 200          p 30  20         p  3  2