Aktuálna revízia z 20:38, 18. september 2023

Literature:

BV: Brejová, Vinař: Metódy v bioinformatike. (preliminary version of lecture notes in Slovak, only several lectures)
DEKM: Durbin, Eddy, Krogh, Mitchison: Biological sequence analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press 1998. Can be studied in the FMFI library under code I-INF-D-21
ZB: Zvelebil, Baum: Understanding Bioinformatics. Taylor & Francis 2008. Can be studied in the FMFI library under code I-INF-Z-2

For each lecture, we list book chapters best corresponding to the covered material. However, the lecture may differ substantially from the listed chapters which serve as the source of additional information.

Recordings of lectures in Slovak from 2018/19

playlist on youtube

Presentations

Notes

L: lecture (everybody), TI: tutorial for computer science/informatics students, TB: tutorial for biology/chemistry/physics students

Sept. 22
L:	Introduction, course rules, sequencing and genome assembly
	BV chapter 1, video 1 video 2
TI:	Introduction to biology	notes
	ZB chapter 1, video
TB:	Introduction to computer science, UCSC genome browser	notes
Sept. 29
L:	Genome assembly 2
	video
TI:	Introduction to probability, genome coverage by sequencing reads	notes
TB:	Introduction to dynamic programming, introduction to probability	notes
Oct. 6
L:	Sequence alignment: Smith-Waterman, Needleman-Wunsch, scoring
	BV chapter 2, DEKM chapter 2.1-2.4, 2.8, ZB chapter 4.1-4.4, 5.1-5.2, video
TI:	Introduction to dynamic programming, proteomics	notes
TB:	Dynamic programming for sequence alignment, dotplots	notes
Oct. 13
L:	Sequence alignment: heuristic alignment (BLAST), statistical significance of alignments, whole genome alignments, multiple alignments
	BV chapter 2, DEKM chapter 2.5, 2.7, 6.1-6.3; ZB chapter 4.5-4.7, 5.3-5.5, video
TI:	Advanced algorithms for sequence alignment	notes
TB:	Programs for sequence alignment, scoring schemes	notes
Oct. 20
L:	Gene finding, hidden Markov models
	BV chapter 4, DEKM chapter 3; ZB chapter 9.3, 10.4-10.7, video
TI:	Fast similarity search, BLAST, MinHash	notes
TB:	Hidden Markov models, E-value	notes
Oct. 27
L:	Phylogenetic tree reconstruction (parsimony, neighbor joining, models of evolution)
	BV chapter 3, DEKM chapter 7,8; ZB chapter 7, 8.1-8.2, video
TI:	Algorithms for HMM	notes
TB:	Substitution models, bootstrap, tree rooting	notes
Nov. 3
L:	Comparative genomics, detection of positive and purification selection, comparative gene finding, phylogenetic HMMs
	BV chapter 5, ZB chapter 9.8, 10.8, video
TI:	Substitution models	notes
TB:	Phylogenetic trees	notes
Nov. 10
L:	Gene expression, clustering, classification, regulatory networks, transcription factors, sequence motifs
	DEKM chapter 5.1, 11.5, ZB chapter 6.6,15.1,16.1-16.5,17.1, SK video
TI:	Felsenstein algorithm, algorithms for HMM and phyloHMM	notes
TB:	K-means clustering, enrichment, multiple testing correction	notes
Nov. 17 No lecture, no tutorials
Nov. 24
L:	Protein structure and function
	DEKM chapter 5; ZB chapter 4.8-4.10, 6.1-6.2, 13.1-13.2, SK video
TI:	Motif finding by EM and Gibbs sampling	notes
TB:	Introduction to context-free grammars	notes
Dec. 1
L:	RNA, secondary structure, Nussinov algorithm, stochastic context-free grammars, RNA family profiles
	DEKM chapter 10, ZB chapter 11.9, SK video
TI:	Examples of biological databases, HW1 model solutions, introduction to context-free grammars	notes
TB:	Genes, comparative genomics, Pfam	notes
Dec. 8
L:	Population genetics
	SK video
TI:	RNA structure	notes
TB:	Example of command-line tools	notes
Dec. 15
L:	Optional journal club presentations
TI:	Course summary, protein threading via integer linear programming	notes
TB:	PSI-BLAST, microarray data, RNA structure, MEME, transcription factors in SGD, population genetics, course summary, graphs	notes

Obsah

1 CI01
- 1.1 Úvod do biológie pre informatikov
2 CB01
- 2.1 Používanie počítačov v M 217
- 2.2 UCSC genome browser
- 2.3 Sekvenovanie v UCSC genome browseri
- 2.4 QUAST: program na štatistiky o kvalite poskladania genómu
- 2.5 Prehľad systému Galaxy
3 CI02
- 3.1 Úvod do pravdepodobnosti
- 3.2 Počítanie pokrytia genómov
- 3.3 Zhrnutie
4 CB02
- 4.1 Úvod do pravdepodobnosti
  - 4.1.1 Pravdepodobnostný model náhodnej sekvencie
- 4.2 Použitie pravdepodobnosti na analýzu potrebného pokrytia pri sekvenovaní
- 4.3 Dynamické programovanie
5 CI03
- 5.1 Dynamické programovanie
- 5.2 Uvod do proteomiky
- 5.3 Sekvenovanie proteinov pomocou MS/MS
- 5.4 Detekcia znamych proteinov pomocou MS (nerobili sme)
6 CB03
- 6.1 Zarovnávanie sekvencií, opakovanie
- 6.2 Dotploty
- 6.3 Dynamické programovanie v Exceli
7 CI04
- 7.1 Opakovanie dynamického programovania pre globálne zarovnanie
- 7.2 Reprezentácia pomocou grafu
- 7.3 Krátka vsuvka o acyklických orientovaných grafoch
- 7.4 Lokálne zarovnanie
- 7.5 Afínne skóre medzier
  - 7.5.1 Nesprávne riešenie pomocou dynamického programovania
  - 7.5.2 Správne riešenie pomocou dynamického programovanania
- 7.6 Lineárna pamäť: Hirshbergov algoritmus 1975
- 7.7 Vypísanie všetkých najlepších riešení
8 CB04
- 8.1 Skórovacie matice
  - 8.1.1 Cvičenie pri počítači
- 8.2 Praktické cvičenie pri počítači: dotploty
  - 8.2.1 Yass a dotploty
- 8.3 Príklady praktických programov
- 8.4 Praktické cvičenie pri počítači: BLAT vs BLAST
  - 8.4.1 BLAT/BLAST
  - 8.4.2 RT PCR sekvencia z cvičenia vyššie
- 8.5 Uniprot pre projekty
9 CI05
- 9.1 Vzorec na výpočet senzitivity jadra
- 9.2 Minimizery: ako ušetriť pamäť a čas
- 9.3 MinHash
10 CB05
- 10.1 Príklady stavových automatov pre HMM
- 10.2 E-hodnota (E-value) zarovnania
- 10.3 Praktické príklady z predchádzajúcich cvičení
11 CI06
- 11.1 HMM opakovanie
- 11.2 Inferencia pomocou najpravdepodobnejšej cesty, Viterbiho algoritmus
- 11.3 Inferencia - dopredný algoritmus
- 11.4 Inferencia - posterior decoding
- 11.5 Trénovanie HMM
- 11.6 Tvorba stavového priestoru modelu
- 11.7 Zovšeobecnené HMM
- 11.8 Párové HMM (pair HMM)
12 CB06
- 12.1 Fylogenetické stromy
- 12.2 Bootstrap
  - 12.2.1 Príklad
- 12.3 Opakovanie pravdepodobnostných modelov
- 12.4 Zložitejšie evolučné modely
- 12.5 Fitchov algoritmus
13 CI07
- 13.1 Substitučné modely - odvodenie
- 13.2 Substitučné matice, zhrnutie
- 13.3 Použitie na odhad evolučnej vzdialenosti
- 13.4 Zložitejšie modely
14 CB07
- 14.1 Praktická ukážka tvorby stromov
- 14.2 Gény, evolúcia a komparatívna genomika v UCSC genome browseri (cvičenie pri počítači)
- 14.3 Objavenie génu HAR1 pomocou komparatívnej genomiky
  - 14.3.1 Cvičenie pri počítači
- 14.4 Praktické príklady z predchádzajúcich cvičení
- 14.5 Hľadanie génov
15 CI08
- 15.1 Felsensteinov algoritmus 1981
16 CB08
- 16.1 Zhlukovanie
  - 16.1.1 K-Means
    - 16.1.1.1 Algoritmus
- 16.2 Nadreprezentacia, obohatenie (enrichment)
- 16.3 Multiple testing correction
- 16.4 Nadreprezentácia, Uniprot (cvičenie pri počítači)
17 CI09
- 17.1 Hľadanie motívov zadefinovaných pravdepodobnostnou maticou
  - 17.1.1 EM algoritmus
  - 17.1.2 Gibbsovo vzorkovanie (Gibbs sampling)
- 17.2 Vzorkovanie z pravdepodobnostného modelu vo všeobecnosti
- 17.3 Poriadnejšie Gibbsovo vzorkovanie pre motívy
18 CB09
- 18.1 Bezkontextové gramatiky
19 CI10
- 19.1 Úvod do bioinformatických databáz a on-line nástrojov
- 19.2 Úvod do bezkontextových gramatík
20 CB10
- 20.1 RNA štruktúra
- 20.2 PSI BLAST
- 20.3 Nadreprezentácia, Uniprot (cvičenie pri počítači)
- 20.4 Objavenie génu HAR1 pomocou komparatívnej genomiky
  - 20.4.1 Cvičenie pri počítači
- 20.5 Expresia génov
- 20.6 Sekvenčné motívy, program MEME
- 20.7 Kvasinkové transkripčné faktory v SGD
- 20.8 Uniprot
- 20.9 Nussinovovej algoritmus (nerobili sme)
21 CI11
- 21.1 RNA struktura
22 CB11
- 22.1 Ukážka práce v Linuxe
23 CI12
- 23.1 Protein threading
24 CB12
- 24.1 Zhrnutie semestra a úvod do teórie grafov
- 24.2 Populacna genomika v UCSC genome browseri
  - 24.2.1 Zopar zaujimavych polymorfizmov v ludskom genome
- 24.3 Sekvenčné motívy, program MEME
- 24.4 Gény, evolúcia a komparatívna genomika v UCSC genome browseri (cvičenie pri počítači)

CI01

Úvod do biológie pre informatikov

Užitočné linky

Pesničky z YouTube:

CB01

Používanie počítačov v M 217

V textovom menu pri štarte zvoľte Linux, v prihlasovacom menu zadajte užívatela bioinf, heslo dostanete
Na dolnom okraji obrazovky je lišta s často používanými nástrojmi, napr. internetový prehliadač Firefox
Vo Firefoxe si otvorte stránku predmetu http://compbio.fmph.uniba.sk/vyuka/mbi/ čast Prednášky a poznámky, nalistujte materiály k dnešnému cvičeniu

UCSC genome browser

On-line grafický nástroj na prezeranie genómov
Konfigurovateľný, veľa možností, ale pomerne málo organizmov
V programe Firefox choďte na stránku UCSC genome browser http://genome-euro.ucsc.edu/ (európsky mirror stránky http://genome.ucsc.edu/ )
Hore v modrom menu zvoľte Genomes, potom zvoľte ľudský genóm verzia hg38. Do okienka search term zadajte HOXA2. Vo výsledkoch hľadania (Known genes) zvoľte gén homeobox A2 na chromozóme 7.
- Pozrime si spolu túto stránku
- V hornej časti sú ovládacie prvky na pohyb vľavo, vpravo, približovanie, vzďaľovanie
- Pod tým schéma chromozómu, červeným vyznačená zobrazená oblasť
- Pod tým obrázok vybranej oblasti, rôzne tracky
- Pod tým zoznam všetkých trackov, dajú sa zapínať, vypínať a konfigurovať
- Po kliknutí na obrázok sa často zobrazí ďalšia informácia o danom géne alebo inom zdroji dát
- V génoch exony hrubé, UTR tenšie, intróny vodorovné čiary

Koľko má HOXA2 exónov? Na ktorom chromozóme a pozícii je? Pozor, je na opačnom vlákne. Ako je táto skutočnosť naznačená na obrázku?
V tracku GENCODE kliknite na gén, mali by ste sa dostať na stránku popisujúcu jeho rôzne vlastnosti. Čo ste sa dozvedeli o jeho funkcii?
- Na tejto stránke nájdite linku na stiahnutie proteínovej sekvencie. Aké sú prvé štyri aminokyseliny?

Sekvenovanie v UCSC genome browseri

Vráťte sa na UCSC genome browser http://genome-euro.ucsc.edu/
Pozrieme si niekoľko vecí týkajúcich sa sekvenovania a skladania genómov
Hore v modrom menu zvoľte Genomes, časť Other
Na ďalšej stránke zvoľte človeka a pomocou menu Human Assembly zistite, kedy boli pridané posledné dve verzie ľudského genómu (hg19 a hg38)
Na tej istej stránke dole nájdete stručný popis zvolenej verzie genómu. Pre ktoré oblasti genómu máme v hg38 najviac alternatívnych verzií? (haplotypes)
Prejdite na región chr21:31,250,000-31,300,000 v hg19 touto linkou: [1]
Zapnite si tracky Mapability a RepeatMasker na "full"
Mapability: nakoľko sa daný úsek opakuje v genóme a či teda vieme jednoznačne jeho čítania namapovať pri použití Next generation sequencing
Ako a prečo sa pri rôznych dĺžkach čítaní líšia? (Keď kliknete na linku "Mapability", môžete si prečítať bližšie detaily.)
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 v hg19: [2] 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ĺžka aj tak nebola známa, nakoľko nebola osekvenovaná).
Cez menu položku View, In other genomes si pozrite, ako zobrazený úsek vyzerá vo verzii hg38. Ako sa zmenila dĺžka z pôvodných 300kb?
Prejdite na genóm Rhesus, verzia rheMac2, región chr7:59,022,000-59,024,000 [3], zapnite si tracky Contigs, Gaps, Quality scores
Aké typy problémov v kvalite sekvencie v tomto regióne vidíte?
- Opäť si môžete pozrieť, či sa problémy odstránili a ako sa zmenila dĺžka sekvencie v najnovšej verzii rheMac8

QUAST: program na štatistiky o kvalite poskladania genómu

Čítania technológie Illumina MiSeq z 500kbp oblasti genómu E.coli
Boli poskladané programom SPAdes [4]
Vzniknuté kontigy sú tu: [5]
Pozrime si štatistiky tohto poskladania v nástroji Quast,
- Predpočítané výsledky a report v pdf formáte
Teraz si tento výpočet spustíme v systéme Galaxy

Prehľad systému Galaxy

https://usegalaxy.eu/
Obsahuje veľa bioinformatických nástrojov, ktoré môžete spúšťať
Ale na výsledky treba niekedy dlho čakať
V ľavom stĺpci hľadanie nástroja alebo nahrávanie dát
V pravom stĺpci zoznam nahratých dát, bežiacich programov a hotových výsledkov (výsledky si pozriete ikonou oka alebo stiahnete ikonou diskety)
V strede nastavenia nástroja alebo prezeranie výsledkov
Pri serióznom používaní odporúčam vytvoriť si konto a prihlásiť sa

Stiahnite si kontigy [6], uložte ako súbor
V ľavom menu zvolíme Upload Data a nahráme stiahnutý súbor
V časti Tools v ľavom menu zadáme do vyhľadávania Quast, zvolíme Quast
- Ako Contigs/scaffolds file zadáme nahratý súbor, ostatné položky necháme predvolené, stlačíme Execute
- Predpočítané výsledky [7]
Druhá analýza: porovnanie poskladaných kontigov so správnou odpoveďou (ak je známa)
- Skutočná E.coli sekvencia, ktorú sme chceli dostať
- Dá sa zadať do nástroja Quast, ak zvolíte Yes v Use a reference genome a tento súbor nahráte ako Reference genome

Ďalšie dáta pre záujemcov:

Použité čítania: prvé čítania z páru, druhé čítania z páru
Galaxy obsahuje aj program SPAdes na skladanie

CI02

Úvod do pravdepodobnosti

Myšlienkový experiment, v ktorom vystupuje náhoda, napr. hod ideálnou kockou/mincou
Výsledkom experimentu je nejaká hodnota (napr. číslo, alebo aj niekoľko čísel, reťazec)
Túto neznámu hodnotu budeme volať náhodná premenná
Zaujíma nás pravdepodobnosť, s akou náhodná premenná nadobúda jednotlivé možné hodnoty
T.j. ak experiment opakujeme veľa krát, ako často uvidíme nejaký výsledok

Príklad 1: hodíme idealizovanou kockou, premenná X bude hodnota, ktorú dostaneme

Možné hodnoty 1,2,..,6, každá rovnako pravdepodobná
Pišeme napr. Pr(X=2)=1/6

Príklad 2: hodíme 2x kockou, náhodná premenná X bude súčet hodnôt, ktoré dostaneme

Možné hodnoty: 2,3,...,12
Každá dvojica hodnôt (1,1), (1,2),...,(6,6) na kocke rovnako pravdepodobná, t.j. pravdepodobnosť 1/36
Súčet 5 môžeme dostať 1+4,2+3,3+2,4+1 - t.j. P(X=5) = 4/36
Súčet 11 môžeme dostať 5+6 alebo 6+5, t.j. P(X=11) = 2/36
Rozdelenie pravdepodobnosti: (tabuľka udávajúca pravdepodobnosť pre každú možnú hodnotu)

hodnota i:   2     3     4     5     6     7     8     9     10    11    12
Pr(X=i):    1/36  2/36  3/36  4/36  5/36  6/36  5/36  4/36  3/36  2/36  1/36

Overte, ze súčet pravdepodobností je 1

Stredná hodnota E(X):

priemer z možných hodnôt váhovaných ich pravdepodobnosťami
v našom príklade $E(X)=2\cdot {\frac {1}{36}}+3\cdot {\frac {2}{36}}+4\cdot {\frac {3}{36}}+5\cdot {\frac {4}{36}}+6\cdot {\frac {5}{36}}+7\cdot {\frac {6}{36}}+8\cdot {\frac {5}{36}}+9\cdot {\frac {4}{36}}+10\cdot {\frac {3}{36}}+11\cdot {\frac {2}{36}}+12\cdot {\frac {1}{36}}=7$
Ak by sme experiment opakovali veľa krát a zrátali priemer hodnôt X, ktoré nám vyšli, dostali by sme číslo blízke E(X)
Iný výpočet strednej hodnoty:
- X=X1+X2, kde X1 je hodnota na prvej kocke a X2 je hodnota na druhej kocke
- $E(X_{1})=1\cdot {\frac {1}{6}}+...+6\cdot {\frac {1}{6}}=3.5$ , podobne aj E(X2) = 3.5
- Platí, že E(X1+X2)=E(X1) + E(X2) a teda E(X) = 3.5 + 3.5 = 7
- Pozor, pre súčin a iné funkcie takéto vzťahy platiť nemusia, napr. $E(X_{1}\cdot X_{2})$ nie je vždy $E(X_{1})\cdot E(X_{2})$

Počítanie pokrytia genómov

Pozrite tiež grafy k pravdepodobnosti: pdf
Náš problém: spočítanie pokrytia
- G = dĺžka genómu, napr. 1 000 000 (predpokladajme, že je cirkulárny)
- N = počet čítaní (readov), napr. 10 000
- L = dlzka čítania, napr. 1000
- Celková dĺžka čítaní NL, pokrytie (coverage) NL/G, v nasom pripade 10x
- V priemere každá báza pokrytá 10x
- Niektoré sú ale pokryté viackrát, iné menej.
- Zaujímajú nás otázky typu: koľko báz očakávame, že bude pokrytých menej ako 3x?
- Dôležité pri plánovaní experimentov (aké veľké pokrytie potrebujem na dosiahnutie určitej kvality)

Pokrytie genómu: predpokladáme, že každé čítanie začína na náhodnej pozícii zo všetkých možných G
Takže ak premenná Y_i bude začiatok i-tého čítania, jej rozdelenie bude rovnomerné
- P(Y_i=1) = P(Y_i=2) = ... = P(Y_i=G) = 1/G

Aká je pravdepodobnosť že nejaké konkrétne i-te čítanie pokrýva konkrétnu pozíciu j?
- P(Y_i>=j-L+1 and Y_i<=j) = P(Y_i=j-L+1)+...+P(Y_i=j) = L/G, označme túto hodnotu p, v našom príklade p=0.001 (1 promile)

Uvazujme premennu X_j, ktora udava pocet čítaní pokryvajucich poziciu j
- mozne hodnoty 0..N
- i-te čítanie pretina poziciu j s pravdepodobnostou p=L/G
- 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 ju 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 $P(X_{j}=k)=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
Stredna hodnota poctu baz v celom genome s pokrytim k je G*P(X_i=k)
- V priemere teda ocakavame 45 baz nepokrytych, 2800 pokrytých menej ako 3 krát a pod.
- Takyto graf, odhad, vieme lahko spravit pre rozne pocty čítaní a tak naplanovat, kolko čítaní potrebujeme

Chceme tiež odhadnúť počet kontigov (podľa článku E.S. Lander and M.S. Waterman. "Genomic mapping by fingerprinting random clones: a mathematical analysis." Genomics 2.3 (1988): 231-239 [8])

Ak niekoľko báz vôbec nie je pokrytých čítaniami, preruší sa kontig
Vieme, koľko báz je v priemere nepokrytých, ale niektoré môžu byť vedľa seba
Nový kontig vznikne aj ak sa susedné čítania málo prekrývajú
Predpokladajme, že na spojenie dvoch čítaní potrebujeme prekryv aspoň T=50
Nech p je pravdepodobnosť, ze dané čítanie i bude posledné v kontigu
Aby sa to stalo, žiadne čítanie j!=i nesmie začínať v prvých L-T bázach kontigu i
Každé čítanie tam začína s pravdepodobnosťou q=(L-T)/G
Ak X je počet čítaní, ktoré zacinaju v tomto useku, tak p = Pr(X=0) = (1-q)^(N-1) podla binomickeho rozdelenia
v priemere ich tam zacne E(X) = (N-1)(L-T)/G co je zhruba N(L-T)/G
Jednoduchší vzorec pre p dostaneme ak binomické rozdelenie premennej X aproximujeme Poissonovým s parametrom $\lambda =N(L-T)/G$ (t.j. aby mali rovnakú strednú hodnotu)
V Poissonovom rozdelení p = Pr(X=0) = exp(-lambda) = exp(-N(L-T)/G)
Presnosť aproximácie: pre parametre N,L,G,T uvedené vyššie dostaneme z binomického rozdelenia p=7.459e-5, z Poissonovho 7.485e-5
Pre N čítaní dostaneme priemerný počet kontigov N*p = N*exp(-N(L-T)/G)
NL/G je pokrytie, N(L-T)/G je pokrytie, ak by sme dĺžku každého čítania skrátili o dĺžku prekryvu
Pre T=50 dostaneme priemerný počet koncov kontigov 0.75 (ak pokryjeme celý kruh, máme nula koncov, preto je hodnota menšia ako 1). Ak znížime N na 5000 (5x pokrytie) dostaneme 43 kontigov

Môže sa zdať zvláštne, ze pri priemernom pocte nepokrytych baz 45 mame pocet koncov v priemere menej ako jedna. Situacia je vsak taka, ze pri opakovaniach tohto experimentu casto dostavame jeden suvisly kontig, ale ak je uz aspon jeden koniec kontigu, byva tam pomerne velka medzera. Tu je napriklad 50 opakovani expertimentu s T=0, priemerny pocet koncov je 0.55, priemerny pocet nepokrytych baz je 49.

nepokr: 0 koncov: 0     nepokr: 0 koncov: 0     nepokr: 0 koncov: 0      
nepokr: 274 koncov: 2   nepokr: 282 koncov: 1   nepokr: 0 koncov: 0      
nepokr: 0 koncov: 0     nepokr: 0 koncov: 0     nepokr: 8 koncov: 1      
nepokr: 0 koncov: 0     nepokr: 12 koncov: 1    nepokr: 0 koncov: 0      
nepokr: 122 koncov: 1   nepokr: 135 koncov: 1   nepokr: 111 koncov: 1    
nepokr: 13 koncov: 1    nepokr: 1 koncov: 1     nepokr: 56 koncov: 1     
nepokr: 265 koncov: 1   nepokr: 0 koncov: 0     nepokr: 10 koncov: 1     
nepokr: 0 koncov: 0     nepokr: 0 koncov: 0     nepokr: 130 koncov: 1    
nepokr: 217 koncov: 1   nepokr: 3 koncov: 1     nepokr: 0 koncov: 0      
nepokr: 0 koncov: 0     nepokr: 0 koncov: 0     nepokr: 86 koncov: 1     
nepokr: 139 koncov: 2   nepokr: 0 koncov: 0     nepokr: 0 koncov: 0      
nepokr: 76 koncov: 1    nepokr: 221 koncov: 1   nepokr: 26 koncov: 1     
nepokr: 0 koncov: 0     nepokr: 1 koncov: 1     nepokr: 0 koncov: 0      
nepokr: 0 koncov: 0     nepokr: 0 koncov: 0     nepokr: 0 koncov: 0      
nepokr: 0 koncov: 0     nepokr: 0 koncov: 0     nepokr: 12 koncov: 1     
nepokr: 103 koncov: 2   nepokr: 0 koncov: 0     nepokr: 71 koncov: 1     
nepokr: 69 koncov: 1    nepokr: 0 koncov: 0

Tento jednoduchy model nepokryva vsetky faktory:
- čítania nemaju rovnaku dlzku
- Problemy v zostavovani kvoli chybam, opakovaniam a pod.
- čítania nie su rozlozene rovnomerne (cloning bias a pod.)
- Vplyv koncov chromozomov pri linearnych chromozomoch
- 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

Pravdepobnostny 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 jej 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 čítaní dlzky L=1000
Nahodna premenna X_i je pocet čítaní 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 čítaní potrebujeme osekvenovat, aby napr. aspon 95% pozicii bolo pokrytych aspon 4 čítaniami

CB02

Úvod do pravdepodobnosti

Myšlienkový experiment, v ktorom vystupuje náhoda, napr. hod ideálnou kockou/korunou
Výsledkom experimentu je nejaká hodnota (napr. číslo, alebo aj niekoľko čísel, reťazec)
Túto neznámu hodnotu budeme volať náhodná premenná
Zaujíma nás pravdepodobnosť, s akou náhodná premenná nadobúda jednotlivé možné hodnoty
T.j. ak experiment opakujeme veľa krát, ako často uvidíme nejaký výsledok

Príklad 1: hodíme idealizovanou kockou, premenná X bude hodnota, ktorú dostaneme

Možné hodnoty 1,2,..,6, každá rovnako pravdepodobná
Píšeme napr. Pr(X=2)=1/6

Príklad 2: hodíme 2x kockou, náhodná premenná X bude súčet hodnôt, ktoré dostaneme

Možné hodnoty: 2,3,...,12
Každá dvojica hodnôt (1,1), (1,2),...,(6,6) na kocke rovnako pravdepodobná, t.j. pravdepodobnosť 1/36
Súčet 5 môžeme dostať 1+4,2+3,3+2,4+1 - t.j. P(X=5) = 4/36
Súčet 11 môžeme dostať 5+6 alebo 6+5, t.j. P(X=11) = 2/36
Rozdelenie pravdepodobnosti: (tabuľka udávajúca pravdepodobnosť pre každú možnú hodnotu)

hodnota i:   2     3     4     5     6     7     8     9     10    11    12
Pr(X=i):    1/36  2/36  3/36  4/36  5/36  6/36  5/36  4/36  3/36  2/36  1/36

Overte, ze súčet pravdepodobností je 1

Stredná hodnota E(X):

priemer z možných hodnôt váhovaných ich pravdepodobnosťami
v našom príklade $E(X)=2\cdot {\frac {1}{36}}+3\cdot {\frac {2}{36}}+4\cdot {\frac {3}{36}}+5\cdot {\frac {4}{36}}+6\cdot {\frac {5}{36}}+7\cdot {\frac {6}{36}}+8\cdot {\frac {5}{36}}+9\cdot {\frac {4}{36}}+10\cdot {\frac {3}{36}}+11\cdot {\frac {2}{36}}+12\cdot {\frac {1}{36}}=7$
Ak by sme experiment opakovali veľa krát a zrátali priemer hodnôt X, ktoré nám vyšli, dostali by sme číslo blízke E(X)
Iný výpočet strednej hodnoty:
- X=X1+X2, kde X1 je hodnota na prvej kocke a X2 je hodnota na druhej kocke
- $E(X_{1})=1\cdot {\frac {1}{6}}+...+6\cdot {\frac {1}{6}}=3.5$ , podobne aj E(X2) = 3.5
- Platí, že E(X1+X2)=E(X1) + E(X2) a teda E(X) = 3.5 + 3.5 = 7
- Pozor, pre súčin a iné funkcie takéto vzťahy platiť nemusia, napr. $E(X_{1}\cdot X_{2})$ nie je vždy $E(X_{1})\cdot E(X_{2})$

Pravdepodobnostný model náhodnej sekvencie

Napríklad chceme modelovať náhodnú DNA sekvenciu dĺžky n s obsahom GC 40%
Máme vrece s guľôčkami označenými A,C,G,T, pričom guľôčok označených A je 30%, C 20%, G 20% a T 30%.
Vytiahneme guľôčku, zapíšeme si písmeno, hodíme ju naspäť, zamiešame a opakujeme s ďalším písmenom atď, až kým nevygenerujeme n písmen

Vytiahnime z mechu 2x guľôčku. Prvé písmeno, ktoré nám vyjde, označme X1 a druhé 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. šanca, že dostaneme sekvenciu AC po dvoch ťahoch je 6%
- Ak rátame pravdepodobnosť, že sa dve nezávislé udalosti stanú súčasne, ich pravdepodobnosti násobíme. V tomto prípade to, či X1=A je nezávislé od toho, či X2=C
Pr(X1 je A alebo C) = Pr(X1=A)+Pr(X1=C) = 0.3+0.2 = 0.5
- Pravdepodobnosť, že prvé písmeno bude A alebo C je 50%
- Pravdepodobnosti navzájom sa vylučujúcich udalostí (X1=A a X1=C) sa môžu sčítať, čím dostaneme pravdepodobnosť, že aspoň jedna z nich nastane
Pr(v sekvencii je aspoň jedno A) = Pr(X1=A alebo X2=A) nemôžeme počítať ako Pr(X1=A)+Pr(X2=A), lebo sa navzájom nevylučujú a prípad, že X1=A a X2=A by sme započítali dvakrát
Správne je Pr(X1 je A alebo X2 je A) = Pr(X1=A) + Pr(X1 <> A a X2=A) = Pr(X1=A) + Pr(X1 <> A) * Pr(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.
Ak u označíme pravdepodobnosť u = Pr(X1=A)=Pr(X1=T)=Pr(X2=A)=Pr(X2=T) a v=Pr(X1=C)=Pr(X1=G)=Pr(X2=C)=Pr(X2=G), aký bude vzorec pre Pr(X1=X2)?

Príklad použitia modelu: Máme krátky primer AACAT. Koľko bude mať v priemere výskytov v sekvencii dĺžky 1000 v našom modeli?

Pravdepodobnosť, ze AACAT je v náhodnej sekvencii hneď na začiatku je Pr(X1=A a X2=A a X3=C a X4=A a X5=A) = 0.3*0.3*0.2*0.3*0.3 = 0.00162
Rovnaká pravdepodobnosť aj na pozícii 2,3,...996
Nech V je počet výskytov v celej sekvencii (náhodná premenná s možnými hodnotami 0,1,...,996, aj keď napr. 996 to určite nemôže byť)
Ideálne by sme chceli spočítať celú tabuľku pravdepodobností pre V, ale uspokojíme sa aj so strednou hodnotou E(V)
Nech Vi je počet výskytov na pozícii i (co je vzdy 0 alebo 1)
$V=V_{1}+V_{2}+\dots +V_{{996}}=\sum _{{i=1}}^{{996}}V_{i}$
$E(V)=E(V_{1})+E(V_{2})+\dots +E(V_{{996}})=996E(V_{1})$
$E(V_{1})=0\cdot \Pr(V_{1}=0)+1\cdot \Pr(V_{1}=1)=\Pr(V_{1}=1)=0.00162$
E(V) = 996*0.00162 = 1.61352
Takze primer AACAT sa v priemere bude v náhodnej sekvencii dĺžky 1000 s 40% obsahom GC vyskytovať v priemere cca 1,6 krát
Primery byvaju dlhsie, takze sanca nahodnych vyskytov je ovela mensia, co je to co vacsinou chceme (chceme primer cielit na konkretnu poziciu, nie na vela nahodnych zhod)

Použitie pravdepodobnosti na analýzu potrebného pokrytia pri sekvenovaní

Pozri cvičenia pre informatikov

Dynamické programovanie

Túto techniku uvidíme na ďalšej prednáške na hľadanie zarovnaní (alignments)
Uvažujme problém platenia pomocou najmenšieho počtu mincí
Napr. máme mince hodnoty 1,2,5 centov, z každej dostatok kusov
Ako môžeme zaplatiť určitú sumu, napr. 13 centov, s čo najmenším počtom mincí?
Aké je riešenie? 5+5+2+1 (4 mince)
Všeobecná formulácia:
- Vstup: hodnoty k mincí m_1,m_2,...,m_k a cieľová suma X (všetko kladné celé čísla)
- Výstup: najmenší počet mincí, ktoré potrebujeme na zaplatenie X
V našom príklade k=3, m_1 = 1, m_2 = 2, m_3 = 5, X=13
Jednoduchý spôsob riešenia: použi najväčšiu mincu, ktorá je najviac X, odčítaj od X, opakuj
Príklad: najprv použijeme mincu 5, zostane X=8, použijeme opäť mincu 5, zostane X=3, použijeme mincu 2, zostane X=1, použijeme mincu 1.
Nefunguje vždy: zoberme mince hodnôt 1,3,4. Pre X=6 najlepšie riešenie je 2 mince: 3+3, ale náš postup (algoritmus) nájde 3 mince 4+1+1
Ukážeme si algoritmus na základe dyn. programovania, ktorý pre každý vstup nájde najlepšie riešenie
Zrátame najlepší počet mincí nielen pre X, ale pre všetky možné cieľové sumy 1,2,3,...,X-1,X
To sa zdá byť ťažšia úloha, ale ukáže sa, že z riešenia pre menšie sumy vieme zostaviť riešenie pre väčšie sumy, takže nám to vlastne pomôže
Spravíme si tabuľku, kde si pre každú sumu i=0,1,2,...X pamätáme A[i]=najmenší počet mincí, ktoré treba na vyplatenie sumy i
Ukážme si to na príklade s mincami 1,3,4

i      0    1    2    3    4    5    6    7    8    9  
A[i]   0    1    2    1    1    2    2    2    2    3

Nevypĺňali sme ju žiadnym konkrétnym postupom, nejde o algoritmus
Ale predstavme si, ze teraz chceme vyplniť A[10].
V najlepšom riešení je prvá minca, ktorú použijeme 1,3, alebo 4
ak je prvá minca 1, máme ešte zaplatiť sumu 10-1=9, tú podľa tabuľky vieme najlepšie zaplatiť na 3 mince, takže potrebujeme 4 mince na zaplatenie 10
ak je prvá minca 3, máme ešte zaplatiť 10-3 = 7, na čo potrebujeme podla tabuľky 2 mince, takže spolu 3 mince na zaplatenie 10
ak je prvá minca 4, máme ešte zaplatiť 10-4 = 6, na čo treba 2 mince, t.j. 3 mince na 10
Nevieme, ktorá z týchto možností je naozaj v najlepšom riešení, ale pre druhé dva prípady dostávame menej mincí, takže výsledok budu 3 mince (napr. 3+3+4)
Zovšeobecníme: A[i] = min { A[i-1]+1, A[i-3]+1, A[i-4]+1 }
A[11] = min { 3+1, 2+1, 2+1} = min {4, 3, 3 } = 3
Pre sústavu mincí 1,2,5, máme A[i] = 1+ min { A[i-1], A[i-2], A[i-5] }
Vo všeobecnosti A[i] = 1+ min { A[i-m_1], A[i-m_2], ..., A[i-m_k] }
Vzorec treba modifikovať pre malé hodnoty i, ktoré sú menšie ako najväčšia minca, lebo A[-1] a pod. nie je definované
Pre zaujímavosť, program v Pythone, stačí meniť hodnoty m a X:

m = [1,3,4]
X = 11
k = len(m)
nekonecno = math.inf
A = [0]
for i in range(1, X + 1):
  min = nekonecno
  for j in range(k):
     if i >= m[j] and A[i - m[j]] < min:
       min = A[i - m[j]]
  A.append(1 + min)
print(A)

Ako nájsť, ktoré mince použiť?
Pridáme druhú tabuľku B, kde v B[i] si pamätáme, ktorá bola najlepšia prvá minca, keď sme počítali A[i] (ak je viac možností, zoberieme ľubovoľnú, napr. najväčšiu)

i      0    1    2    3    4    5    6    7    8    9   10   
A[i]   0    1    2    1    1    2    2    2    2    3    3
B[i]   -    1    1    3    4    4    3    4    4    4    4

Potom ak chceme nájsť napr. mince pre 10, vidíme, že prvá bola B[10]=4. Zvyšok je 6 a prvá minca na vyplatenie 6 je B[6]=3. Zostáva nám 3 a B[3]=3. Potom nám už zostáva 0, takže sme hotoví. Takže najlepšie vyplatenie je 4+3+3
Celý program v Pythone:

m = [1,3,4]
X = 11
k = len(m)
nekonecno = math.inf
A = [0]
B = [-1]
for i in range(1, X + 1):
  min = nekonecno
  min_minca = -1
  for j in range(k):
     if i >= m[j] and A[i - m[j]] < min:
       min = A[i - m[j]]
       min_minca = m[j]
  A.append(1 + min)
  B.append(min_minca)

while X > 0:
    print(B[X])
    X = X - B[X]

Dynamické programovanie vo všeobecnosti

Okrem riešenia celého problému vyriešime aj veľa menších podproblémov
Riešenia podproblémov ukladáme do tabuľky
Pri riešení väčšieho podproblému používame už vypočítané hodnoty pre menšie podproblémy

Aká je časová zložitosť tohto algoritmu?

Dva parametre: X a k
Tabuľka veľkosti O(X), každé políčko čas O(k). Celkovo O(Xk)

CI03

Dynamické programovanie

Pozri Cvičenia pre biológov

Uvod do proteomiky

Viac informacii: [9], Bafna, Reinert 2004
Pozri tiež prezentáciu k cvičeniu

Gélová elektroforéza (gel electrophoresis) - uvedene pre zaujimavost, nerobili sme

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
- v niektorých prípadoch vieme sekvenciu proteínu určiť priamo z MS/MS, bez databázy proteínov

Sekvenovanie proteinov pomocou MS/MS

Vsetky hmotnosti budeme povazovat za cele cisla

Vstup:

celková hmotnosť peptidu M,
hmotnosti aminokyselín a[1],...,a[20],
spektrum ako tabuľka f[0],...,f[M], ktorá hmotnosti m určí skóre f[m] podľa signálu v okolí príslušného bodu grafu

Označenie:

Uvažujme postupnosť aminokyselín $x=x_{1}\dots x_{k}$
Nech $m(x)=\sum _{{j=1}}^{k}a[x_{j}]$ je hmotnosť x
Nech $M_{P}(x)=\{m(x_{1}\dots x_{j})\mid j=1,\dots ,k\}$ sú hmotnosti prefixov x
Nech $M_{S}(x)=\{m(x_{j}\dots x_{k})\mid j=1,\dots ,k\}$ sú hmotnosti sufixov x

Problém 1

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

Výstup:

postupnosť aminokyselín x taká, že $m(x)=M$ a $\sum _{{m\in M_{P}(x)}}f[m]$ je maximálna možná
Chceme teda najst peptid, ktory maximalizuje sucet skore svojich prefixov

Riešenie

Dynamicke programovanie s podproblemom S[m] je skore najlepsieho prefixu s hmotnostou m
Rekurencia? Zlozitost? Je to polynomialny algoritmus? (Aky velky je vlastne vstup?)

Problém 2

Berme do uvahy aj y-iony, ktore meraju hmotnost sufixu, scitame skore prefixov a sufixov

Výstup:

postupnosť aminokyselín x taká, že $m(x)=M$ a $\sum _{{m\in M_{P}(x)}}f[m]+\sum _{{m\in M_{S}(x)}}f[m]$ je maximálna možná

Riešenie

pouzijeme upravenu skorovaciu tabulku g[m]=f[m]+f[M-m] a algoritmus pre problem 1

Problem tejto formulacie:

jeden signal sa moze ratat dvakrat, raz ako b-ion, raz ako y-ion, algoritmus ma tendenciu pridavat taketo artefakty

Problém 3

Ak hmotnost nejakeho prefixu a nejakeho sufixu su rovnake, zarataj ich skore iba raz (skore peptidu je skore mnoziny hmotnosti jeho prefixov a sufixov)

Výstup:

postupnosť aminokyselín x taká, že $m(x)=M$ a $\sum _{{m\in M_{P}(x)\cup M_{S}(x)}}f[m]$ je maximálna možná

Riesenie:

Ina formulacia: maximalizujeme $\sum _{{m\in M_{p}(x)\cup M_{S}(x),m\leq M/2}}h[m]$
$h[m]=\left\{{\begin{array}{ll}f[m]+f[M-m]&{\mbox{ak }}m<M/2\\f[m]&{\mbox{ak }}m=M/2\end{array}}\right.$
Definuj novy podproblem: S[p,s] je najlepsie skore, ktore moze dosiahnut prefix s hmotnostou p a sufix s hmotnostou s, kde 0<=p,s<=M/2,
Rekurencia

$S[p,s]=\left\{{\begin{array}{ll}\max _{{i=1\dots 20}}S[p,s-a[i]]+h[s]&{\mbox{ak }}p<s\\\max _{{i=1\dots 20}}S[p-a[i],s]+h[p]&{\mbox{ak }}p>s\\\max _{{i=1\dots 20}}S[p-a[i],s]&{\mbox{ak }}p=s\\\end{array}}\right.$

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 (nerobili sme)

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 so 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?
- Zlozitejsi alg. s predspracovanim pomocou FFT (Fast Fourier Transform) Bansal, Cieliebak, Liptak 2004

CB03

Zarovnávanie sekvencií, opakovanie

Uvažujme skórovanie zhoda +3, nezhoda -1, medzera -2
Reťazce TAACGG a CACACT

Globálne zarovnanie

Rekurencia: A[i,j] = max {A[i-1,j]-2, A[i,j-1]-2, A[i-1,j-1]+s(x_i, y_j) }, pričom A[0,i]=-2i, A[i,0]=-2i

        C   A   C   A   C   T
    0  -2  -4  -6  -8  -10  -12
T  -2  
A  -4  
A  -6  
C  -8  
G  -10  
G  -12

Lokálne zarovnanie

Rekurencia: A[i,j] = max {0, A[i-1,j]-2, A[i,j-1]-2, A[i-1,j-1]+s(x_i, y_j) }, pričom A[0,i]=0, A[i,0]=0

        C   A   C   A   C   T
    0   0   0   0   0   0   0
T   0 
A   0 
A   0 
C   0 
G   0 
G   0

Dotploty

Dotplot je graf, ktory ma na kazdej osi jednu sekvenciu a ciarky zobrazuju lokalne zarovnania (cesty v matici)
Niekoľko príkladov dotplotov: pdf
Prvé príklady dotplotov porovnávajú rôzne mitochondriálne 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 prekroci urcenu hodnotu
http://emboss.bioinformatics.nl/cgi-bin/emboss/dotmatcher
Vyskusame protein escargot voci sebe s hodnotami http://pfam.xfam.org/protein/ESCA_DROME window 8 threshold 24
Pomocou YASSu vyskusame kluster zhlukov PRAME z ludskeho genomu

Dynamické programovanie v Exceli

Práca so vzorcami v tabuľkovom procesore (Excel, LibreOffice, ...)

Okrem konkrétnych hodnôt, napr. 0.3, môžu byť aj vzorce, ktoré začínajú =, napr =0.3*0.3 dá do políčka 0.09 (* znamená násobenie)
Vo vzorcoch môžeme používať aj hodnoty z iných políčok, napr. =A2+B2 dáme do políčka C2, zobrazí sa tam súčet
Ak políčko so vzorcom skopírujeme do iného políčka, Excel sa snaží uhádnuť, ako zmeniť vzorec
- Ak sme v C2 mali =A2+B2 a skopírovali sme to do C3, vzorec sa zmení na =A3+B3
Ak niektoré adresy políčok majú zostávať rovnaké aj pri kopírovaní, dáme pred písmeno aj číslo $,
- Ak v C2 máme =A2+$B$2 a skopírujeme to do C3, dostaneme =A3+B2
Dolár môžeme dať aj pred iba jednu súradnicu (stĺpec alebo riadok), tá sa potom nebude pri kopírovaní meniť

Späť k rozmieňaniu mincí

Vráťme sa k príkladu s rozmieňaním mincí a skúsme si ho "naprogramovať" v Exceli, resp. spreadsheet aplikácii v OpenOffice
Vseobecna formulacia:
- Vstup: hodnoty k minci m_1,m_2,...,m_k a cielova suma X (vsetko kladne cele cisla)
- Vystup: najmensi pocet minci, ktore potrebujeme na zaplatenie X
My pouzijeme mince hodnot 1,3,4
Spravime si tabulku, kde si pre kazdu sumu i=0,1,2,...X pamatame A[i]=najmensi pocet minci, ktore treba na vyplatenie sumy i

i      0    1    2    3    4    5    6    7    8    9  
A[i]   0    1    2    1    1    2    2    2    2    3

vzorec A[i] = min { A[i-1]+1, A[i-3]+1, A[i-4]+1 }
aby sme nemuseli zvlast uvazovat hodnoty mensie ako 4, (kde sa neda A[i-4]), urcime si A[-1], A[-2] atd ako nejake velke cislo (napr 100), takze vzorec plati pre vsetky i>0

i      -4  -3  -2  -1  0    1    2    3    4    5    6    7    8    9  
A[i]  100 100 100 100  0    1    2    1    1    2    2    2    2    3

v exceli si najskor spravime horny riadok tabulky
- do nejakeho policka (napr, B4) zapiseme prvu hodnotu (-4)
- do susedneho C4 zapiseme vzorec =B4+1, dostaneme hodnotu -3
  - vzorce zacinaju znamienkom =
  - B4 je suradnica policka o jedno vlavo, k nej pripocitame 1
- policko C4 nakopirujeme do riadku kolkokrat chceme, dostaneme hodnoty -2, -1, 0, 1,...
  - kopirovat sa da tahanim laveho dolneho rohu okienka
  - vzorec sa automaticky posuva na =C4+1, =D4+1, atd
- o riadok nizsie do B5..E5 napiseme hodnotu 100 (okienka A[-4]..A[-1])
- do F5 dame 0 (okienko A[0] nasej tabulky)
- do G5 napiseme vzorec =MIN(F5+1,D5+1,C5+1), t.j. A[1] = min(A[1-1]+1,A[1-3]+1,A[1-4]+1)
- tento vzorec potom nakopirujeme do riadku tabulky
- F5 sa bude posuvat na G5, H5,... a podobne ostatne dva cleny

Cvičenie:

Ako by sme zmenili na inu mincovu sustavu, napr. 1,2,5?
Stiahnite si subor zo stranky predmetu a skuste si tuto zmenu urobit [10]

Zarovnávanie sekvencií v Exceli

skusme si dynamicke programovanie pre globalne zarovnanie naprogramovat v Exceli
budeme postupovat podobne ako pri minciach, ale potrebujeme dve specialne funkcie: MID(text,od,dlzka) z textu vyberie urcitu cast. Pomocou toho si vstupny text rozdelime na jednotlive pismena, ktore si napiseme do zahlavia tabulky
vsimnite si pouzivanie dolarov v nazvoch policok: ak je pred menom stlpca alebo riadku $, tento sa neposuva ked vzorec kopirujem do inych policok
IF(podmienka,hodnota1,hodnota2) vyberie bud hodnotu 1 ak je podmienka splnena alebo hodnotu2 ak nie je. Napr IF(F$8=$B12 ,1,-1) zvoli skore +1 ak sa hodnota v F8 rovna hodnote v B12 a skore -1 ak sa nerovnaju.

Cvičenie:

Zmente tabulku tak, aby skore pre zhody, nezhody a medzery bolo dane bunkami B1, B2 a B3 tabulky. Staci zmenit vzorce a policka D9, C10 a D10 a nakopirovat do zvysku tabulky. Ake bude skore najlepsieho zarovnania sekvencii AACGTA a ACACCTA ak skore nezhody je -2 a medzery -3?
Ako treba zmenit vzorce, aby sme pocitali lokalne zarovnanie?
Subor najdete tu

CI04

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]-1\right\}$

Ako presne by sme implementovali?
Ako spočítame maticu spätných šípok B?
Aká je časová a pamäťová zložitosť?

Reprezentácia pomocou grafu

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

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

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

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

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

Lokálne zarovnanie

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

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

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

Afínne skóre medzier

Napr. otvorenie medzery o=-3, pokracovanie medzery e=-1

A  -  -  -  T  C  G
A  C  G  C  T  C  C
1 -3 -1 -1  1  1  -1

Nesprávne riešenie pomocou dynamického programovania

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

$A[i,j]=\max \left\{A[i-1,j-1]+s(x_{i},y_{j}),A[i-1,j]+c(i-1,j,hore),A[i,j-1]+c(i,j-1,vlavo)\right\}$

c(i,j,s) = o, ak v policku A[i,j] mame sipku s
c(i,j,s) = e, ak v policku A[i,j] mame inu sipku

Preco toto riesenie nefunguje?

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

Toto je obvykla chyba pri dynamickom programovani:

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

Správne riešenie pomocou dynamického programovanania

Riesenie 1:

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

Riesenie 2:

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

Lineárna pamäť: Hirshbergov algoritmus 1975

Klasicke dynamicke programovanie potrebuje cas O(nm)
Trivialna implementacia tiez pouzije pamat O(mn) - uklada si celu maticu A, pripadne maticu B so sipkami naspat
Na vypocet matice A nam z stacia dva riadky tejto matice: riadok i pocitam len pomocou riadku i-1, starsie viem zahodit
Ale ak chcem aj vy[isat zarovnanie, stale potrebujem pamat O(mn) na maticu sipok B
Hirschbergov algoritmus znizi pamat na O(m+n), zhruba zdvojnasobi cas (stale O(mn))

Prejdeme celú maticu a spočítame maticu A. Zapamätáme si, kde moja cesta prejde cez stredný riadok matice
- Nech B_k[i,j] je najväčší index v riadku k, cez ktorý prechádza najkratšia cesta z (0,0) do (i,j)

Ako vieme B_k[i,j] spočítať?
- ak A[i,j] = A[i-1,j-1]+w(S[i],T[j])$, potom B_k[i,j]=B_k[i-1,j-1].
- ak A[i,j]=A[i-1,j]+1, potom B_k[i,j]=B_k[i-1,j].
- ak A[i,j]=A[i,j-1]+1, potom B_k[i,j]=B_k[i,j-1]
- Toto platí, ak i > k. Pre i=k nastavíme B_k[i,j]=j

Ak už poznáme A[i-1,*] a B_k[i-1,*], vieme spočítať A[i,*] a B_k[i,*].
- Stacia nam teda iba dva riadky matice A a B_k
Nech k'=B_k[m,n]. Potom v optimálnom zarovnaní sa S[1..k] zarovná s T[1..k'] a S[k+1..m] s T[k'+1..n].
- Toto použijeme na rekurzívny algoritmus na výpočet zarovnania:

optA(l1, r1, l2, r2) { // align S[l1..r1] and T[l2..r2]
    if(r1-l1 <= 1 ||  r2-l2 <=1) 
        solve using dynamic programming
    else {
        k=(r-l+1)/2;
        for (i=0; i<=k; i++) 
           compute A[i,*] from A[i-1,*]
        for (i=k+1; i<=r-l+1; i++) 
           compute A[i,*], B_k[i,*] from A[i-1,*], B_k[i-1,*]
        k2=B_k[r1-l1-1,r2-l2-1];
        optA(l1, l1+k-1, l2, l2+k2-1); 
        optA(l1+k, r2, l2+k2, r2); 
    }
}

Casova zlozitost:

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

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

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

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

CB04

Skórovacie matice

Chceme určiť skórovaciu schému pre zarovnávanie dvoch DNA sekvencií (bez medzier). Máme dva modely, každý z nich vie vygenerovať 2 zarovnané sekvencie dĺžky n.

Model R (random) reprezentuje nezávislé náhodne sekvencie

Použijeme naše vrece s guličkami označenými A,C,G,T, pričom guličiek označených A je 30%, C 20%, G 20% a T 30%.
Vytiahneme guličku, zapíšeme si písmeno, hodíme ju naspäť, zamiešame a opakujeme s ďalším písmenom atď až kým nevygenerujeme n písmen pre jednu sekvenciu a n písmen pre druhú
Máme jednu sekvenciu ACT a druhú ACC. Aká je pravdepodobnosť, že práve tieto sekvencie vygenerujeme v našom modeli R?
Nezávislé udalosti pre jednotlivé písmená, t.j. Pr(X1=A)*Pr(X2=C)*Pr(X3=T)*Pr(Y1=A)*Pr(Y2=C)*Pr(Y3=C) = 0.3*0.2*0.3*0.3*0.2*0.2 = 0.000216
Spolu máme v modeli $4^{6}=4096$ možností ako vygenerovať 2 DNA sekvencie dĺžky 3

Model H (homolog) reprezentuje zarovnanie vzájomne súvisiacich sekvencií

máme vrece, v ktorom je napr.
- po 21% guličiek označených AA, TT
- po 14% označených CC, GG
- po 2.4% označených AC, AG, CA, CT, GA, GT, TC, TG
- po 3.6% označených AT, TA
- po 1.6% označených CG, GC.
Spolu máme 70% guličiek označených rovnakými písmenami, 30% rôznymi

n krát z vreca vytiahneme guličku a písmená píšeme ako stĺpce zarovnania A1, A2,.., An.
aká je pravdepodobnosť, ze dostaneme ACT zarovnané s ACC?
Pr(A1=AA)*Pr(A2=CC)*Pr(A3=TC) = 0.21*0.14*0.024 = 0.0007056

Skóre zarovnania je log Pr(zarovnania v H)/Pr(zarovnania v R), t.j. log (0.0007056 / 0.000216) = 0.514105 (pre desiatkový logaritmus)

kladné skóre znamená, že model H lepšie zodpovedá dátam (zarovnaniu) ako model R
záporné skóre znamená, že model R lepšie zodpovedá dátam

Cvičenie pri počítači

Stiahnite si súbor nižšie, uložte si ho a otvorte v Exceli / OpenOffice / LibreOffice
V záložke Matica vyplňte do žltej oblasti vzorce na výpočet pravdepodobnosti krátkeho zarovnania, logaritmus pomeru pravdepodobnosti a súčet skóre, pričom vo vzorcoch použijete odkazy na políčka v riadkoch 9-13, stĺpcoch B a E
Súčet skóre by mal zhruba rovný desaťnásobku logaritmu pomeru - prečo vidíme rozdiely?

Potom skúšajte meniť %GC a %identity v horných riadkoch tabuľky a pozrite sa, ako to ovplyvní skórovanie. Výsledné skóre zo stĺpca E ručne prepíšte (bez formúl) do tabuľky v záložke Výsledky. Prečo nastávajú také zmeny ako vidíte?

Praktické cvičenie pri počítači: dotploty

Yass a dotploty

Program Yass hlada lokalne zarovnania v DNA sekvenciach, zobrazuje vo forme dot plotov
V novom okne/tabe si otvorte YASS server na adrese http://bioinfo.lifl.fr/yass/yass.php
V dalsom okne si na stranke UCSC genome browseru si zobrazte oblast chr21:9,180,027-9,180,345 vo verzii hg38 ludskeho genomu [11]
- tento región obsahuje Alu repeat. Tieto opakovania tvoria cca 10% ľudského genómu, viac ako milión kópií
- zobrazte si DNA sekvenciu tohto useku takto: na hornej modrej liste zvolte View, potom v podmenu 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 sa presunte na poziciu chr21:8,552,000-8,562,000 (10kb sekvencie na chromozome 21, s niekolkymi vyskytmi Alu) [12]
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?
For whole genomes, you can try D-GENIES based on minimap2 program
- For example, two yeast genomes:
  - ftp://ftp.ebi.ac.uk/pub/databases/ena/wgs/public/uid/UIDE01.fasta.gz [13]
  - ftp://ftp.ebi.ac.uk/pub/databases/ena/wgs/public/cab/CABVLU01.fasta.gz [14]

Príklady praktických programov

Pozrime sa na niekolko nastrojov, vsimnime si, ake poskytuju nastavenia a co vypisuju na vystupe, dajme to do suvisu s prednaskami

viacnasobne zarovnania neskor

Plné dynamické programovanie

Balicek emboss, obsahuje programy na klasicke dynamicke programovanie (needle - globalne, water - lokalne), najdu sa na stranke EBI http://www.ebi.ac.uk/Tools/psa/
porovnanie lokalneho a globalneho zarovnania
- Dva proteiny z rôznych kvasiniek zarovnáme lokálne, globálne a globálne s tým, že neplatíme za medzery na koncoch
sekvencie a vysledne zarovnania: #CB-aln-dp
vo vysledku si vsimnime, kolko ma kazde z nich %identity, %gaps, a kam sa zarovna sekvencia IRESPLGG ktora je na pozicii 29 v prvom a 30 v druhom proteine

Lokalne zarovnanie
Length: 588
Identity:     170/588 (28.9%)
Similarity:   270/588 (45.9%)
Gaps:         116/588 (19.7%)
Score:  611.0

Globalne zarovnanie
Length: 650
Identity:     178/650 (27.4%)
Similarity:   282/650 (43.4%)
Gaps:         153/650 (23.5%)
Score: 588.5

Globalne zarovnanie s nulovou penaltou za medzeru na koncoch
Length: 651
Identity:     177/651 (27.2%)
Similarity:   282/651 (43.3%)
Gaps:         155/651 (23.8%)
Score: 608.0

NCBI Blast

NCBI BLAST http://blast.ncbi.nlm.nih.gov/ vela roznych nastrojov (porovnavanie DNA vs proteiny, pripadne translacia DNA na protein v 6 ramcoch)
- Heuristicky algoritmus, moze niektore zarovnania vynechat
- rozne nastavenia, vystup E-value

Low complexity masking: nepouzivat pri hladani jadier zarovnania regiony v ktorych sa velakrat opakuje ta ista aminokyselina

Priklad (z ucebnice Zvelebil and Baum):

>sp|P04156|PRIO_HUMAN Major prion protein OS=Homo sapiens GN=PRNP PE=1 SV=1
MANLGCWMLVLFVATWSDLGLCKKRPKPGGWNTGGSRYPGQGSPGGNRYPPQGGGGWGQP
HGGGWGQPHGGGWGQPHGGGWGQPHGGGWGQGGGTHSQWNKPSKPKTNMKHMAGAAAAGA
VVGGLGGYMLGSAMSRPIIHFGSDYEDRYYRENMHRYPNQVYYRPMDEYSNQNNFVHDCV
NITIKQHTVTTTTKGENFTETDVKMMERVVEQMCITQYERESQAYYQRGSSMVLFSSPPV
ILLISFLIFLIVG

Hladajme v databaze Reference sequence (Refseq), organizmus human

Bez maskovania vypise napr aj toto zarovnanie:

>ref|NP_065842.1| serine/threonine-protein kinase TAO1 isoform 1 [Homo sapiens]
Length=1001

 Score = 45.1 bits (105),  Expect = 1e-06, Method: Composition-based stats.
 Identities = 26/61 (43%), Positives = 27/61 (44%), Gaps = 11/61 (18%)

Query  38   YPGQGSPGGNRYPPQGGGG--WGQPHGG---GWGQPHGGG---WGQPHGGGWGQPHGGGWG  90
            YPG     G  + P GG G  WG P GG    WG P  GG   WG P G   G P G   G
Sbjct  904  YPGAS---GWSHNPTGGPGPHWGHPMGGPPQAWGHPMQGGPQPWGHPSGPMQGVPRGSSMG  961

 Score = 40.0 bits (92),  Expect = 4e-05, Method: Composition-based stats.
 Identities = 25/62 (40%), Positives = 25/62 (40%), Gaps = 10/62 (16%)

Query  26   PKPGGW--NTGGSRYPGQGSPGGNRYPPQGGGGWGQPHGGG---WGQPHGGGWGQPHGGGWG  82
            P   GW  N  G   P  G P G   PPQ    WG P  GG   WG P G   G P G  
Sbjct  905  PGASGWSHNPTGGPGPHWGHPMGG--PPQA---WGHPMQGGPQPWGHPSGPMQGVPRGSSMG  961

Ak zapneme maskovanie, toto zarovnanie uz nenajde, v zarovnani sameho so sebou sa objavia male pismena alebo Xka:

>ref|NP_000302.1|major prion protein preproprotein [Homo sapiens]
Length=253

 Score =   520 bits (1340),  Expect = 0.0, Method: Compositional matrix adjust.
 Identities = 253/253 (100%), Positives = 253/253 (100%), Gaps = 0/253 (0%)

Query  1    MANLGCWMLVLFVATWSDLGLCKKRPKPGGWNTGGSRYPGQGSPGGNRYppqggggwgqp  60
            MANLGCWMLVLFVATWSDLGLCKKRPKPGGWNTGGSRYPGQGSPGGNRYPPQGGGGWGQP
Sbjct  1    MANLGCWMLVLFVATWSDLGLCKKRPKPGGWNTGGSRYPGQGSPGGNRYPPQGGGGWGQP  60

Query  61   hgggwgqphgggwgqphgggwgqphgggwgqgggTHSQWNKPSKPKTNMKHMagaaaaga  120
            HGGGWGQPHGGGWGQPHGGGWGQPHGGGWGQGGGTHSQWNKPSKPKTNMKHMAGAAAAGA
Sbjct  61   HGGGWGQPHGGGWGQPHGGGWGQPHGGGWGQGGGTHSQWNKPSKPKTNMKHMAGAAAAGA  120

Query  121  vvgglggymlgsamsRPIIHFGSDYEDRYYRENMHRYPNQVYYRPMDEYSNQNNFVHDCV  180
            VVGGLGGYMLGSAMSRPIIHFGSDYEDRYYRENMHRYPNQVYYRPMDEYSNQNNFVHDCV
Sbjct  121  VVGGLGGYMLGSAMSRPIIHFGSDYEDRYYRENMHRYPNQVYYRPMDEYSNQNNFVHDCV  180

Query  181  NITIKQHtvttttkgenftetDVKMMERVVEQMCITQYERESQAYYQRGSSMVLFSsppv  240
            NITIKQHTVTTTTKGENFTETDVKMMERVVEQMCITQYERESQAYYQRGSSMVLFSSPPV
Sbjct  181  NITIKQHTVTTTTKGENFTETDVKMMERVVEQMCITQYERESQAYYQRGSSMVLFSSPPV  240

Query  241  illisfliflivG  253
            ILLISFLIFLIVG
Sbjct  241  ILLISFLIFLIVG  253

BLAT, chains, nets v UCSC browseri

Program BLAT v UCSC browseri http://genome-euro.ucsc.edu/ rychlo vyhladava sekvencie v genome, ale nevie najst slabsie podobnosti
- Vhodne pouzitie: zarovnanie mRNA ku genomu, presne urcenie suradnic nejakej sekvencie, a pod.
Net tracky v UCSC genome browseri nam umoznuju prechadzat medzi homologickymi oblastami roznych genomov

Praktické cvičenie pri počítači: BLAT vs BLAST

BLAT/BLAST

Sekvencia uvedena nizsie vznikla pomocou RT-PCR na ľudských cDNA knižniciach
Choďte na UCSC genome browser http://genome-euro.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 sliepky sa vyskytuje homologický úsek.
Skusme tu istu sekvenciu zarovnat ku genomu sliepky programom Blat: 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
Aká je dĺžka, identity a E-value najlepšieho zarovnania? Na ktorom je chromozóme?

RT PCR sekvencia z cvičenia vyššie

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

Uniprot pre projekty

pozrime si ľudský proteín PTPRZ1 v databáze Uniprot http://www.uniprot.org/

CI05

Vzorec na výpočet senzitivity jadra

Uvažujme jadro dĺžky k (ako v programe BLAST pre nukleotidy, na prednáške sa dĺžka jadra označovala w, teraz k)
Uvažujme pravdepodobnostný model zarovnania, v ktorom má každá pozícia pravdepodobnosť p, že bude zhoda a (1-p), ze bude nezhoda alebo medzera, zarovnanie ma dlzku L
Nahodna premenna X_i = 1 ak na pozicii i je zhoda, 0 inak
Nahodna premenna Y_i = 1 ak na pozicii i zacina jadro, t.j. ak $X_{i}=1,X_{{i+1}}=1,\dots ,X_{{i+k-1}}=1$
$P(Y_{i}=1)=p^{k}$ , nakolko X_i su navzajom nezavisle
Nech $Y=\sum _{{i=1}}^{{L-k-1}}Y_{i}$
Z linearity strednej hodnoty vieme lahko odhadnut $E(Y)=(L-k+1)p^{k}$
Nas ale zaujima P(Y>0) = 1-P(Y=0)
$P(Y=0)=P(Y_{1}=0\wedge \dots \wedge Y_{{L-k+1}}=0)$
Preco neplati, ?
- Jednotlive Y_i nie su nezavisle, napr. $P(Y_{{i+1}}=1|Y_{i}=1)=p$
- V postupnosti Y_i sa jednotky maju tendenciu zhlukovat spolu
Pravdepodobnost nepritomnosti jadra P(Y=0) ale vieme spocitat dynamickym programovanim
Nech A[n] je pravdepodobnost nepritomnosti jadra v prvých n stlcoch zarovnania (0<=n<=L)
Budeme rozlisovat pripady podla toho, kolko je na konci X_1..X_n jednotiek
- Tento pocet moze byt 0..k-1 (ak by bol >=k, mali by sme vyskyt jadra)
$A[n]=\left\{{\begin{array}{ll}1&{\mbox{ak }}n<k\\\sum _{{i=0}}^{{k-1}}p^{i}(1-p)A[n-i-1]&{\mbox{ak }}n\geq k\\\end{array}}\right.$

V druhom riadku $p^{i}(1-p)$ zodpoveda $P(X_{1}...X_{n}{\mbox{ konci presne }}i{\mbox{ jednotkami}})$ a A[n-i-1] je $P(X_{1}...X_{{n-i-1}}{\mbox{ neobsahuje jadro}})$ , ale to je to iste ako $P(X_{1}...X_{n}{\mbox{ neobsahuje jadro }}|X_{1}...X_{n}{\mbox{ konci presne }}i{\mbox{ jednotkami}})$

Minimizery: ako ušetriť pamäť a čas

k-merom nazveme k za sebou iducich pismen v nejakom retazci (nukleotidov DNA)
Zakladne pouzite jadier: pri porovnavani dvoch sekvencii (alebo mnozin sekvencii) uloz vsetky k-mery jednej sekvencie do slovnika (napr. hash tabulky), potom prechadzaj vsetky k-mery druhej sekvencie a hladaj ich v slovniku

Priklad k=5,
- DB vľavo, k-mery a ich pozície uložené do slovníka,
- query vpravo, k-mery hľadané v slovníku, šípkou vyznačené nájdené k-mery (jadrá)

AGTGGCTGCCAGGCTGG    cGaGGCTGCCtGGtTGG  
AGTGG,0              CGAGG              
 GTGGC,1              GAGGC             
  TGGCT,2              AGGCT <-          
   GGCTG,3              GGCTG <-
    GCTGC,4              GCTGC <-          
     CTGCC,5              CTGCC <-       
      TGCCA,6              TGCCT        
       GCCAG,7              GCCTG        
        CCAGG,8              CCTGG
         CAGGC,9              CTGGT     
          AGGCT,10             TGGTT    
           GGCTG,11             GGTTG   
            GCTGG,12             GTTGG

Trik na znizenie potrebnej pamate (napr. program BLAT): ukladaj iba kazdy s-ty k-mer z prvej sekvencie, potom hladaj vsetky k-mery z druhej
Trochu znizi aj senzitivitu, ale kedze jadra sa zhlukuju, mame sancu aspon jedno jadro zo zhluku najst
Zarucene najdeme jadro, ak mame aspon k+s-1 zhod za sebou
Mohli by sme v query tiež hľadať iba každý s-ty k-mer?
- Čo ak by db a query boli to isté, iba v query ba chýbalo prvé písmeno?

Priklad k=5, s=3, k-mery nalavo sa ulozia, k-mery napravo sa hladaju, najde sa jedno jadro

AGTGGCTGCCAGGCTGG    cGaGGCTGCCtGGtTGG  
AGTGG,0              CGAGG              
   GGCTG,3            GAGGC             
      TGCCA,6          AGGCT
         CAGGC,9        GGCTG <-        
            GCTGG,12     GCTGC          
                          CTGCC         
                           TGCCT        
                            GCCTG       
                             CCTGG
                              CTGGT
                               TGGTT    
                                GGTTG   
                                 GTTGG

Prefikanejsia idea je minimizer: uvazuj vsetky skupiny s po sebe iducich k-merov (sliding window), v kazdej skupine najdi abecedne najmensi k-mer (minimizer) a uloz do slovnika
Pri posune okna o jedno doprava casto najmensi k-mer zostava ten isty a netreba ho znovu ukladat, cim sa usetri pamat (a čas)
Rozdiel je pri hladani: v slovniku nehladame vsetky k-mery druhej sekvencie, ale tiez len minimizery, co moze usetrit cas
Zarucene najdeme jadro, ak mame aspon k+s-1 zhod za sebou

Priklad k=5, s=3, vlavo: do slovnika sa ulozia k-mery s cislom; vpravo: v slovniku sa hladaju k-mery s hviezdickou, najde sa jedna zhoda

AGTGGCTGCCAGGCTGG    cGaGGCTGCCtGGtTGG  
AGTGG,0              CGAGG              
 GTGGC                GAGGC             
  TGGCT                AGGCT*           
   GGCTG,3              GGCTG           
    GCTGC,4              GCTGC          
     CTGCC,5              CTGCC* <--    
      TGCCA                TGCCT        
       GCCAG                GCCTG       
        CCAGG,8              CCTGG*     
         CAGGC,9              CTGGT*    
          AGGCT,10             TGGTT    
           GGCTG                GGTTG*  
            GCTGG                GTTGG

Obzvlast vyhodne ak prva a druha mnozina sekvencii su ta ista, napr. pri hladani prekryvov v citaniach pri skladani genomu. Kazde citanie ma mnozinu minimizerov, ktore sa pouziju ako kluce v slovniku, hodnoty su zoznamy citani. Dvojice citani zdielajuce nejaky minimizer (binning) sa dostanu do jedneho zoznamu a budu uvazovane pri vypocte vzajomneho prekryvu
V praxi sa do slovnika neuklada lexikograficky najmensi k-mer, ale kazdy k-mer sa prehasuje nejakou funkciou f a zoberie sa ten s minimalnou hodnotou
Dovod je, ze sa chceme vyhnut, aby minimizermi boli casto sekvencie typu AAAAA, ktore su v biologickych sekvenciach nadreprezentovane a casto nie su funkcne zaujimave
Priemerna hustota minimizerov v sekvencii pri nahodnom hashovani je cca 2/(s+1)⁠ (v BLATe bola nizsia, 1/s)
Minimizery vyuziva napr. aj minimap2, velmi popularny nastroj na zarovnavanie citani navzajom a ku genomom
- na zarovnanie nanoporovych citani ku genomu pouziva k=15, s=10, prekryvy v nanoporovych citaniach k=15, s=5, porovnanie genomov s identitou nad 80% k=19, s=10

Li, Heng. "Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences." Bioinformatics 32.14 (2016): 2103-2110. [15]

Roberts M, Hayes W, Hunt BR, Mount SM, Yorke JA. Reducing storage requirements for biological sequence comparison. Bioinformatics. 2004 Dec 12;20(18):3363-9. [16]

Marçais, G., Pellow, D., Bork, D., Orenstein, Y., Shamir, R., & Kingsford, C. (2017). Improving the performance of minimizers and winnowing schemes. Bioinformatics, 33(14), i110-i117. [17]

MinHash

Odbočka do analýzy web-stránok: Podobnosť textov

Majme množinu webových stránok (webová stránka je postupnosť slov). Chceme nájsť medzi nimi dvojice podobných. Ako môžeme definovať podobnosť dvoch textov?

Jeden zo spôsobov ako to spraviť je pozrieť sa na množstvo slov, spoločných pre jednotlivé dvojice stránok. Očakávame, že čím viac spoločných slov majú, tým sú podobnejšie. Túto mieru podobnosti formalizuje matematický pojem Jaccardovej miery podobnosti.

Nech $U$ je univerzum slov a nech $A,B\subseteq U$ sú jeho podmnožinami (t.j. množiny slov dvoch textov). Potom Jaccardova miera podobnosti $J(A,B)$ je definovaná nasledovne:

$J(A,B):={\dfrac {|A\cap B|}{|A\cup B|}}$

Táto miera nadobúda hodnotu 0 iba v prípade, ak množiny sú disjunktné, a hodnotu 1 iba v prípade, že množiny sú totožné. Inak sa jej hodnota nachádza na otvorenom intervale $(0,1)$ , a čím viac spoločných slov majú, tým je jej hodnota vyššia.

Potom otázku "Ktoré dvojice textov sú podobné?" môžeme preformulovať napríklad ako "Ktoré dvojice textov majú Jaccardovu mieru podobnosti vyššiu ako $\alpha$ ?", kde $\alpha \in (0,1)$ je nejaká prahová hodnota.

Ako rýchlo by sme vedeli spočítať Jaccardovu mieru pre dve množiny slov, každú s n prvkami?

Exaktný výpočet Jaccardovej miery podobnosti nie je vždy dostatočne rýchly pre účely konkrétnej aplikácie, takže logickým riešením je pokúsiť sa jej hodnotu vypočítať iba približne (t.j. aproximovať).

Prvá idea: aproximácia vzorkovaním

Ak by sme vedeli vzorkovať z $A\cup B$ prvky u1, u2, ..., u_s (rovnomerne, nezavisle), a pre kazdy prvok u_i by sme vedeli rychlo zistit, ci patri do prieniku, mohli by sme odhadnut J(A, B) pomocou nahodnej premennej X

$X={\frac {1}{n}}\sum _{{i=1}}^{s}X_{i}$ kde X_i=1 ak u_i patri do prieniku a X_i=0 inak

Toto sa podoba na zistovanie oblubenosti politika prieskumom verejnej mienky, u1, u2, ..., u_s su "respondenti"

Pre kazde X_i

$E(X_{i})=0\cdot Pr[X_{i}=0]+1\cdot Pr[X_{i}=1]=Pr[X_{i}=1]=J(A,B)$ .

Z linearity strednej hodnoty

$E(X)=E({\frac {1}{n}}\sum _{{i=1}}^{n}X_{i})={\frac {1}{n}}\sum _{{i=1}}^{n}E[X_{i}]=J(A,B)$ .

Z toho vyplýva, že náhodná premenná $X$ je nevychýleným odhadom Jaccardovej miery.

V štatistike základnou mierou kvality nevychýleného odhadu je jeho disperzia, odvodenie pozri nižšie

$Var(X)={\dfrac {J(A,B)-J^{2}(A,B)}{s}}\leq {\dfrac {1}{4s}}$

Pri zvyšujúcej veľkosti vzorky s teda klesá disperzia. Podobná situácia ako pri prieskumoch verejnej mierky, kde pri väčšom súbore respondentov dostaneme dôveryhodnejšie výsledky.

Problémy:

nie je ľahké rovnomerne vzorkovať z $A\cup B$
na zisťovanie, či u_i je v prieniku potrebujeme mat reprezentaciu A a B v pamati, co moze byt problem pri velkej kolekcii dokumentov
idea: chceme dostat nejake ine premenne X_i, ktore budu nezavisle a ich $E(X_{i})=J(A,B)$ a ktore sa budu lahsie pocitat

Aproximácia Jaccardovej miery: MinHash

Budeme mať náhodné hašovacie funkcie $h_{1},h_{2},\dots ,h_{s}$ .

O každej hašovacej funkcii predpokladáme, že ak ju použijeme na nejakú množinu $A=\{a_{1},a_{2},\ldots ,a_{n}\}$ , tak $h(a_{1}),h(a_{2}),\ldots ,h(a_{n})$ bude náhodná permutácia množiny $A$ .

Pre množinu $A=\{a_{1},a_{2},\ldots ,a_{n}\}$ a hašovaciu funkciu h je $minHash_{{h}}(A)$ je definovaný nasledovne:

$minHash_{h}(A):=\min\{h(a_{1}),h(a_{2}),\ldots ,h(a_{n})\}$

Keďže $h$ je náhodná hašovacia funkcia, tak sa na hodnotu $minHash(A)$ môžeme pozerať ako na náhodnú premennú, ktorá reprezentuje rovnomerne náhodný výber prvku z množiny $A$ .

Nech $X_{i}$ je náhodná premenná, ktorá nadobúda hodnotu 1, ak $minHash_{{h_{i}}}(A)=minHash_{{h_{i}}}(B)$ , a inak hodnotu 0. Potom $E[X_{i}]=J(A,B)$ , lebo celkovo máme $|A\cup B|$ prvkov a $J(A,B)=|A\cap B|/|A\cup B|$ značí, aké percento z nich je v prieniku.

Takýmito hodnotami $X_{i}$ teda nahradíme náhodné vzorky diskutované vyššie.

Algoritmus:

Pre kazdy dokument hasuj kazde slovo s funkciami, najdi minHash pre kazdu funkciu a uloz vektor tychto hodnot ako "sketch" dokumentu. Cas vypoctu pre dokument s n slovami je O(ns)
Pre dva dokumenty porovname vektor po zlozkach a ak najdeme x zhod, J(A,B) odhadneme ako x/s. Cas vypoctu O(s)
Vieme si tiez pre kazdu hasovaciu funkciu spravit slovnik, ktory mapuje minHash do zoznamu dokumentov a budeme porovnavat iba dvojice dokumentov, ktore sa niekde dostali do toho isteho zoznamu (t.j ich odhad J(A,B) bude nenulovy)

Alternativa: namiesto s roznych funkcii pouzijeme iba jednu a vezmeme nielen minimum, ale s najmensich prvkov. Potom J(A,B) odhadneme pomocou $J(S_{A},S_{B})$ kde $s_{A}$ je mnozina hodnot v sketchi mnoziny $A$ . To usetri cas pri vypocte sketchu.

Broder AZ. On the resemblance and containment of documents. InProceedings. Compression and Complexity of SEQUENCES 1997 (Cat. No. 97TB100171) 1997 Jun 13 (pp. 21-29). IEEE. [18]

Hľadanie podobných sekvencií

Ako "slová" použijeme všetky k-mery danej sekvencie. Potom na hľadanie dvoch podobných sekvencií z množiny sekvencií môžeme použiť minhash.

Napríklad Mash používa k=21, s=1000 (s najmenších v jednej funkcii) na porovnávanie genómov, sketch má asi 8kb na genóm (genóm má milióny alebo miliardy nukleotidov)

Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, Phillippy AM. Mash: fast genome and metagenome distance estimation using MinHash. Genome biology. 2016 Dec;17(1):1-4. [19]

Vypocet disperzie

Binarna premenna X_i, kde P(X_i=1)=p. Disperzia $Var(X_{i})=E(X_{i}^{2})-(E(X_{i}))^{2}$ . Spočítajme si postupne obe hodnoty.

$E(X_{i}^{2})=0^{2}\cdot Pr[X_{i}=0]+1^{2}\cdot Pr[X_{i}=1]=Pr[X_{i}=1]=p$

$(E(X_{i}))^{2}=p^{2}$

Čiže $Var(X_{i})=p-p^{2}$ . Aká je maximálna možná hodnota disperzie?

Táto otázka je ekvivalentná otázke "Aké je maximum funkcie $f(x)=x-x^{2}$ na intervale $[0,1]$ ?". Pre určenie extrémov hladkých funkcií treba nájsť korene jej prvej derivácie. Derivácia tejto funkcie je $f'(x)=1-2x$ , jej koreň je hodnota $0.5$ . Hodnota funkcie v tomto bode je rovná $0.25$ . Čiže $Var(X_{i})\leq 0.25$ .

$X={\frac {1}{s}}\sum _{{i=1}}^{s}X_{i}$ kde X_i su nezavisle premenne a každá z nich má strednú hodnotu $E(X_{i})=E(X)=p$ a rovnakú disperziu $Var(X_{i})=Var(X)=p-p^{2}$ . T.j. X je ich priemer.

Pozrieme sa na jej varianciu: $Var(X)=Var\left({\frac {X_{1}+X_{2}+\ldots +X_{s}}{s}}\right)={\dfrac {1}{k^{2}}}Var(X_{1}+X_{2}+\ldots +X_{s}){\overset {*}{=}}{\dfrac {1}{s^{2}}}[Var(X_{1})+\ldots Var(X_{s})]={\dfrac {1}{k^{2}}}s\cdot Var(X_{i})={\dfrac {Var(X_{i})}{s}}\leq {\dfrac {1}{4s}}$

(*) tento prechod je možný len kvôli tomu, že premenné $X_{i}$ sú nezávislé.

Vidíme teda, že varianciu (resp. kvalitu) môžeme potlačiť na ľubovoľne malú postupným zvýšením počtu hashov.

Všimnite si, že premenná $sX$ (t.j. nie priemer, ale súčet jednotlivých $X_{i}$ ) je súčtom nezávislých indikátorov s rovnakou distribúciou, a teda má binomické rozdelenie s parametrami $s$ a $p$ .

CB05

Príklady stavových automatov pre HMM

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

Nakreslite HMM (stavovy diagram), 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
toto sa da rozsirit na HMM, ktory reprezentuje ORF, teda nieco, co zacina start kodonom, potom niekolko beznych kodonov, ktore nie su stop kodonom a na koniec stop kodon

Dalsi biologicky priklad HMM: topologia transmembranovych proteinov.

E-hodnota (E-value) zarovnania

Priklady k tejto casti v prezentacii pdf
Mame dotaz dlzky m, databazu dlzky n, skore najlepsieho lokálneho zarovnania S
E-value je ocakavany pocet zarovnani so skore aspon S ak dotaz aj databaza su nahodne
Hrackarsky priklad: dotaz dlzky m=10, databaza dlzky n=300, S=6
Zoberme nas nahodny model s obsahom GC 50%
Mame vrece s gulockami oznacenymi A,C,G,T, z kazdej 25%
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
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

Vypocet strednej hodnoty vzorcom namiesto simulacie (rychlejsie)

zlozita matematicka teoria [20]
E-value sa priblizne da odhadnut vzorcom:

$E=Kmne^{{-\lambda S}}$

n a m su dlzky porovnavanych sekvencii, S je skore, K a lambda su parametre, ktore zavisia od skorovacej schemy a od frekvencii vyskytu jednotlivych baz v nasom modeli nahodnej sekvencie.
Napr blastn pre skorovaci system zhoda 1, nezhoda -1, medzera -2 používa lambda=0.800, K=0.0640
$E=0.0640nm0.45^{S}$ .
Zdvojnásobenie dĺžky databázy alebo dĺžky dotazu zdvojnásobí E-value
Zníženie skóre o 1 tiež zhruba zdvojnásobí E-value (delenie 0.45, t.j. nasobenie 2.2)
- Cislo, ktorym nasobime, zalezi od lambda a teda od skorovacej schemy a frekvencii vyskytu baz

Praktické príklady z predchádzajúcich cvičení

#CB04#Praktick.C3.A9_cvi.C4.8Denie_pri_po.C4.8D.C3.ADta.C4.8Di:_BLAT_vs_BLAST BLAT vs BLAST

CI06

HMM 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,S_{i}}}$ (**)

Algoritmus:

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

for u=1 to m=počet stavov

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

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 zlogaritmované 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}=u}}\Pr(A_{1},A_{2},...,A_{i}\wedge S_{1},S_{2},...,S_{i})$

$F[1,u]=\pi _{u}e_{{u,S_{1}}}$

$F[i,u]=\sum _{v}F[i-1,v]a_{{v,u}}e_{{u,S_{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 zodpovedá 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 prokaryotických génov
Repeaty v intrónoch: multiple path problem
Intrón má dĺžku aspoň 10

Zovšeobecnené HMM

Predstavme si HMM s dvoma stavmi, napr. gén / negén, pričom každý stav má prechod do seba aj do druhého stavu
Úloha: Nech p je pravdepodobnosť, že zostaneme v tom istom stave, (1-p), že prejdeme do druhého stavu. Aká je pravdepodobnosť, že v stave zostaneme presne k krokov (k>=1)?
- Riešenie: $p^{k}(1-p)$
- Toto rozdelenie sa nazýva geometrické a pravdepodbnosť exponenciálne rýchlo klesá s rastúcim k
Keď sa pozrieme na histogram reálny dĺžkov génov / exónov a iných oblastí, väčšinou sa enpodobá na geometrické rozdelnie, môže priponínať napr. normálne rozdelenie s určitou priemenrou dĺžkou a rozptylom okolo
- Jednoduché HMM teda dobre nemodeluje tento fenomén
Zovšeobecnené HMM (semi-Markov) pracuje tak, že v stave má ľubovoľné rozdelenie pravdepodobnosti dĺžok. Model vôjde do stavu, vygeneruje dĺžku k z tohto rozdelenia, potom vygeneruje k znakov z príslušnej emisnej tabuľky a na záver sa rozhodne, ktorým prechodom opustí stav
Úloha: ako spočítame pravdepodobnosť konkrétnej sekvencie a konkrétnej postupnosti stavov aj s dĺžkami? (zaveďme si aj nejaké vhodné označenie)
Úloha: ako treba upraviť Viterbiho algoritmus pre tento model? Aká bude jeho zložitosť?
- Zložitosť bude kvadraticky rásť od dĺžky sekvencie, predtým rástla lineárne
Predstavme si teraz, že rozdelenie dĺžok má hornú hranicu D takú, že všetky dĺžky väčšie ako D majú nulovú pravdepodobnosť.
- Úloha: ako sa toto obmedzenie prejaví v zložitosti Viterbiho algoritmu?
- Uloha: navrhnite, ako modelovať zovšeobecnený HMM s rozdelením dĺžok ohraničeným D pomocou normálneho stavu, kde sa jedne zovšeobecnený stav nahradí vhodnou postupnosťou D obyčajných stavov.

Párové HMM (pair HMM)

Nebrali sme, uvedene pre zaujimavost

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

CB06

Fylogenetické stromy

Terminológia:

zakorenený strom, rooted tree
nezakorenený strom, unrooted tree
hrana, vetva, edge, branch
vrchol, uzol, vertex, node
list, leaf, leaf node, tip, terminal node
vnútorný vrchol, internal node
koreň, root
podstrom, subtree, clade

Zopár faktov o stromoch

Majme zakorenený strom s n listami, v ktorom má každý vnútorný vrchol 2 deti. Takýto strom vždy má n-1 vnútorných vrcholov a 2n-2 vetiev (prečo?)
Majme nezakorenený strom s n listami, v ktorom má každý vnútorný vrchol 3 susedov. Takýto strom vždy má n-2 vnútorných vrcholov a 2n-3 vetiev
Koľkými spôsobmi môžeme zakoreniť nezakorenený strom s n listami?
- koreň môže byť na hociktorej vetve stromu, teda je 2n-3 možností zakorenenia
Ak nakreslíme zakorenený strom obvyklým spôsobom, listy sú usporiadané zhora nadol (alebo zľava doprava). Koľko rôznych poradí listov vieme dostať rôznym zakresľovaním toho istého stromu s n listami?
- máme n-1 vnutornych vrcholov, v kazdom mozeme vymenit lave a prave dieta. Pre kazdu konfiguraciu takychto vymen dostavame ine poradie, celkovy pocet poradi je 2^n-1
Čo vieme zistiť o príbuznosti organizmov z nezakoreneného stromu (napr. kvartet 4 organizmov)?
- skúsime zakoreniť všetkými spôsobmi a vidíme, že o žiadnych dvoch listoch nevieme povedať, že by boli sesterské (evolučne bližšie než ostatné), lebo koreň stromu môže byť zrovna na niektorej hrane vedúcej ku nim
- vieme však zistiť, že niektoré dvojice sesterské nebudú

Iné použitie stromov v informatike:

uvidíme hierarchické zhlukovanie, bayesovské siete, ale tiež efektívne dátové štruktúry

Bootstrap

Náhodne vyberieme niektoré stĺpce zarovnania, zostrojíme strom vybranou metódou
Celé to opakujeme veľa krát
Značíme si, koľkokrát sa ktorá vetva opakuje v stromoch, ktoré dostávame
- Pri nezakorenených stromoch je vetva rozdelenie listov na dve skupiny (bipartícia)
Nakoniec zostavíme strom z celých dát a pozrieme sa ako často sa ktorá jeho vetva vyskytovala
- Môžeme zostaviť aj strom z často sa vyskytujúcich hrán (napr. tých, ktoré sú vo viac ako 50% stromov)
Bootstrap hodnoty nám dajú určitý odhad spoľahlivosti, hlavne ak máme celkovo málo dát (krátke zarovnanie)
Ak však dáta nezodpovedajú vybranej metóde/modelu, tak aj pre zlý strom môžeme dostať vysoký bootstrap

Detail: pri náhodnom výbere stĺpcov pre bootstrap väčšinou používame rovnaký počet stĺpcov ako pôvodné zarovnanie, pričom vzorkujeme z pôvodného zarovnania s povolením opakovania stĺpca.

Príklad

Robili sme 100x bootstrap, 40x sme dostali strom (i) na obrázku, 40x sme dostali strom (ii) a 20x sme dostali strom (iii)
Strom (iii) sme dostali aj spustením metody na celých dátach
Zistite úroveň bootstrap podpory pre jednotlivé vetvy stromu (iii)
Ktoré ďalšie vetvy majú podporu aspoň 20%?
Aký strom by sme dostali, ak by sme chceli nechať iba vetvy s podporou aspoň 80%?

Opakovanie pravdepodobnostných modelov

Keď počítame pravdepodobnosť, rozmýšľame o myšlienkovom experimente, v ktorom hádžeme kockou, ťaháme gulôčky z vreca a pod.

Dôležité je vždy si poriadne uvedomiť, ako tento experiment prebieha
Tieto myšlienkové experimenty však nastavujeme tak, aby odzrkadľovali nejaké aspekty reality, napr. skutočných DNA sekvencií, ich evolúcie a pod.
Takže pravdepodobnosti, ktoré spočítame v idealizovanom svete nám možno niečo povedia o reálnom svete
Slávny citát štatistika Georga Boxa "All models are wrong, but some are useful."

Aké sme doteraz videli modely

Skórovacie matice: porovnavame model nahodnych sekvencii a model nahodnych zarovnani
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 anotacia 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?
- 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)

Zložitejšie evolučné modely

Jukes-Cantorov model uvažuje len dĺžku hrany udanú ako priemerný počet substitúcií (vrátane tých, ktore nevidíme, kvôli tomu, že boli dve na tom istom mieste)
Nie všetky substitúcie sa dejú rovnako často: tranzície (v rámci pyrimidínov T<->C, v rámci purínov A<->G) sú pravdepodobnejšie ako tranzverzie (A,G)<->(C,T)
Nie všetky nukleotidy sa v danom genóme vyskytujú rovnako často (napr. mitochondriálne genómy majú nízky obsah GC)
Tieto javy zachytáva HKY model (Hasegawa, Kishino & Yano)
Matica rýchlostí zmeny (substitution rate matrix)

$\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 riadku matice ma byt 0, t.j. $\mu _{A}=\beta \pi _{C}+\alpha \pi _{G}+\beta \pi _{T}$
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 rychlosti spocitaju pravdepodobnost, ze baza X zmutuje na bazu Y Pr(Y|X,t)
Napr. pre 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

Fitchov algoritmus

Nerobili sme, uvedene pre zaujimavost

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 rovnaku 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 mutacia na e3 ale urcite ubudne mutacia 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

CI07

Substitučné modely - odvodenie

Nech $P(b|a,t)$ je pravdepodobnosť, že ak začneme s bázou a, tak po čase t budeme mať bázu b.
Pre dané t môžeme také pravdepodobnosti usporiadať do matice 4x4 (ak študujeme DNA):

$S(t)=\left({\begin{array}{cccc}P(A|A,t)&P(C|A,t)&P(G|A,t)&P(T|A,t)\\P(A|C,t)&P(C|C,t)&P(G|C,t)&P(T|C,t)\\P(A|G,t)&P(C|G,t)&P(G|G,t)&P(T|G,t)\\P(A|T,t)&P(C|T,t)&P(G|T,t)&P(T|T,t)\\\end{array}}\right)$

Riadky zodpovedaju povodnej baze a, stlpce novej baze b
Sucet kazdeho riadku je 1

Požiadavky na S(t)

Intuitivne cim vacsie t, tym vacsia pravdepodobnost zmeny, pre nulovy cas este ziadna zmena nemala kedy nastat, mame teda $S(0)=I$ (jednotkova matica)
Naopak ked t ide do nekonecna, kazda baza velakrat zmutovala a teda uz prilis nezalezi, co to bolo na zaciatku. S(t) ma teda v limite pre velke t vsetky riadky rovnake.
$\lim _{{t\rightarrow \infty }}S(t)=\left({\begin{array}{cccc}\pi _{A}&\pi _{C}&\pi _{G}&\pi _{T}\\\pi _{A}&\pi _{C}&\pi _{G}&\pi _{T}\\\pi _{A}&\pi _{C}&\pi _{G}&\pi _{T}\\\pi _{A}&\pi _{C}&\pi _{G}&\pi _{T}\\\end{array}}\right)$
Rozdelenie pravdepodobnosti $\pi$ nazyvame limitne (equilibrium)
Predpokladame tiez, ze pravdepodobnost mutacie zavisi len od aktualnej bazy, nie od minulych stavov a ze charakter procesu mutacii sa v case nemeni. Teda 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})$ . Preto takyto model nazyvame multiplikativny.
Ak by sme uvazovali iba diskretne (celociselne) 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.
Tento typ modelu sa nazyva Markovov retazec so spojitym casom (continuous-time Markov chain)

Jukes-Cantorov substitučný model

Tento model predpoklada, ze vsetky substitucie su rovnako pravdepodobne, matica teda musí vyzerať nejako takto:

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

Matica rýchlostí pre J-C model

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

Pre velmi maly cas t je s(t) velmi male cislo (blizke 0) a pre velmi male s(t) su kvadraticke cleny $s(t)^{2}$ ovela mensie ako linearne cleny s(t) a teda

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

Aj pre ine rozumne male nasobky

$S(c\epsilon )\approx \left({\begin{array}{cccc}1-3cs(\epsilon )&cs(\epsilon )&cs(\epsilon )&cs(\epsilon )\\\dots \end{array}}\right)$

Tento vztah dava zmysel: pre velmi male t mozeme zanedbat moznost, ze baza zmutovala viackrat a pravdepodobnost jednej mutacie linearne rastie s dlzkou casu.

Ak oznacime $t=c\epsilon$ a $\alpha =s(\epsilon )/\epsilon$ dostaneme pre velmi male t

$S(t)\approx \left({\begin{array}{cccc}1-3\alpha t&t\alpha &t\alpha &t\alpha \\\dots \end{array}}\right)$

Vytvorme si teraz maticu rychlosti (intenzit) (transition rate matrix, substitution 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(t)\approx I+Rt$
$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$
V limite dostaneme $S(t)R=\lim _{{\epsilon \rightarrow 0}}(S(t+\epsilon )-S(t))/\epsilon =S'(t)$
Dostali sme diferencialnu rovnicu S(t)R = S'(t), pociatocny stav $S(0)=I$ ).
Nasobenim matic S(t) a R dostavame, ze diagonalny prvok $S(t)R$ je $-3\alpha +12\alpha s(t)$ a nediagonalny $\alpha -4\alpha s(t)$ . Takze dostavame diferencialnu rovnicu $s'(t)=\alpha -4\alpha s(t)$ z rovnosti nediagonalnych prvkov (z rovnosti diagonalnych prvkov dostavame tu istu rovnicu len prenasobenu konstantou -3).

Poriadnejšie odvodenie diferenciálnej rovnice

$s'(t)=\lim _{{\epsilon \to 0}}{\frac {s(t+\epsilon )-s(t)}{\epsilon }}$ z definicie limity
$s(t+\epsilon )=(1-3s(t))s(\epsilon )+s(t)(1-3s(\epsilon ))+s(t)s(\epsilon )+s(t)s(\epsilon )$ podla multiplikativnosti S(t)
po úprave $s(t+\epsilon )=s(\epsilon )+s(t)-4s(t)s(\epsilon )=s(t)+s(\epsilon )(1-4s(t))$
$s'(t)=\lim _{{\epsilon \to 0}}{\frac {s(\epsilon )(1-4s(t))}{\epsilon }}=(1-4s(t))\lim _{{\epsilon \to 0}}{\frac {s(\epsilon )}{\epsilon }}=(1-4s(t))s'(0)$
oznacme $\alpha =s'(0)$ (alfa je konstanta, nezavisi od epsilon ani t)
$s'(t)=\alpha (1-4s(t))$ (finalna diferencialna rovnica, rovnaka ako predtym pre matice)
Riesenie diferencialnej rovnice $s(t)=1/4+ce^{{-4\alpha t}}$ pre kazdu konstantu c
Mozeme overit dosadenim do rovnice, pricom $s'(t)=-4c\alpha e^{{-4\alpha t}}$
c=-1/4 dopocitame z pociatocnej podmienky s(0)=0
overime tiez, ze $s'(0)=\alpha$

Vlastnosti riešenia

Takže máme 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)$

Ked $t\rightarrow \infty$ , dostávame všetky prvky matice rovné 1/4, t.j. $\lim _{{t\to \infty }}s(t)=\lim _{{t\to \infty }}1-3s(t)=1/4$ .
$\alpha$ je teda pravdepodobnosť konkrétnej zmeny za jednotku času, ak uvažujeme veľmi krátke časy alebo presnejšie derivácia prvku s(t) vzhľadom na t v bode 0
Aby sme nemali naraz aj $\alpha$ aj $t$ , zvykneme maticu R normalizovat tak, aby priemerný počet substitúcii za jednotku času bol 1. V prípade Jukes-Cantorovho modelu je to keď $\alpha =1/3$ .

Substitučné matice, zhrnutie

S(t): matica 4x4, kde políčko $S(t)_{{a,b}}=P(b|a,t)$ je pravdepodobnosť, že ak začneme s bázou a, tak po čase t budeme mať bázu b.
Jukes-Cantorov model predpokladá, že táto pravdepodobnosť je rovnaká pre každé dve bázy $a\neq b$
Pre daný čas t máme teda všade mimo diagonály s(t) a na diagonále 1-3s(t)
Matica rýchlostí R: pre Jukes-Cantorov model všade mimo diagonály $\alpha$ , na diagonále $-3\alpha$
Pre veľmi malý čas t je S(t) zhruba I-Rt
Rýchlost alpha je teda pravdepodobnosť zmeny za jednotku casu, ak uvažujeme veľmi krátke časy, resp. derivácia s(t) vzhľadom na t v bode 0
Riešením diferenciálnych rovníc pre Jukes-Cantorov model dostávame $s(t)=(1-e^{{-4\alpha t}})/4$
Matica rýchlostí sa zvykne normalizovať tak, aby na jednotku času pripadla v priemere jedna substitúcia, čo dosiahneme ak $\alpha =1/3$

Použitie na odhad evolučnej vzdialenosti

V case $t$ je pravdepodobnost, ze uvidime zmenenu bazu $D(t)={\frac {3}{4}}(1-e^{{-4\alpha t}})$
V realnom pouziti (vypocet matice vzdialenosti pre metodu spajania susedov) mame dve zarovnane sekvencie, medzi ktorymi vidime zmenenych baz, chceme odhadnut t
- Spatne teda zratame t, ktore by hodnote $D(t)=d$ prinalezalo.
Dostavame teda vzorec pre vzdialenost, ktory sme videli na prednaske $t=-{\frac {3}{4}}\log \left(1-{\frac {4}{3}}d\right)$
Ak $d\rightarrow 0.75$ , dostavame $t\rightarrow \infty$
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.

Zložitejšie modely

V praxi sa používajú komplikovanejsie substitučné modely, ktoré majú všeobecnejšiu maticu rýchlostí R

$R=\left({\begin{array}{cccc}.&\mu _{{AC}}&\mu _{{AG}}&\mu _{{AT}}\\\mu _{{CA}}&.&\mu _{{CG}}&\mu _{{CT}}\\\mu _{{GA}}&\mu _{{GC}}&.&\mu _{{GT}}\\\mu _{{TA}}&\mu _{{TC}}&\mu _{{TG}}&.\\\end{array}}\right)$

Hodnoty na diagonále matice sa dopočítavajú aby súčet každého riadku bol 0.

Hodnota $\mu _{{xy}}$ v tejto matici vyjadruje rýchlosť, akou sa určitá báza x mení na inú bázu y.
Presnejšie $\mu _{{xy}}=\lim _{{t\rightarrow 0}}{\frac {\Pr(y\,|\,x,t)}{t}}$ .

Kimurov model napr. zachytáva, ze puríny sa častejšie menia na iné puríny (A a G) a pyrimidíny na ine pyrimidíny (C a T).

Má dva parametre: rýchlosť tranzícií alfa, transverzií beta

$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) tiež umožnuje rôzne pravdepodobnosti A, C, G a T v ekvilibriu.
Ak nastavíme čas v evolučnom modeli na nekonečno, nezáleží na tom, z ktorej bázy sme začali, frekvencia výskytu jednotlivých báz sa ustáli v tzv. ekvilibriu.
V Jukes-Cantorovom modeli je pravdepodobnosť ľubovoľnej bázy v ekvilibriu 1/4.
V HKY si zvolime aj frekvencie jednotlivých nukleotidov v ekvilibriu $\pi _{A},\pi _{C},\pi _{G},\pi _{T}$ so súčtom 1
Parameter kapa: pomer tranzícií a transverzií (alfa/beta)
Matica rýchlostí:
- $\mu _{{x,y}}=\kappa \pi _{y}$ ak mutácia x->y je tranzícia,
- $\pi _{y}$ ak mutácia x->y je transverzia

Pre zložité modely nevieme odvodiť explicitný vzorec na výpočet S(t), ako sme mali pri Jukes-Cantorovom modeli
Ale vo všeobecnosti pre maticu rýchlostí dostávame .
- Exponenciálna funkcia matice A sa definuje ako $e^{A}=\sum _{{k=0}}^{\infty }{1 \over k!}A^{k}.$
- Ak maticu rychlosti R diagonalizujeme (určite sa dá pre symetrické R) $R=UDU^{{-1}}$ , kde D je diagonálna matica (na jej diagonále budu vlastné hodnoty R), tak $e^{{Rt}}=Ue^{{Dt}}U^{{-1}}$ , t.j. exponenciálnu funkciu uplatníme iba na prvky na uhlopriečke matice D.

CB07

Praktická ukážka 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-data/cb06/cb06-aln.fa
- UCSC browseri si pozrieme usek ludskeho genomu (verzia hg38) chr6:135,851,998-136,191,840 s genom PDE7B (phosphodiesterase 7B) [21]
- Na modrej liste zvolime Tools, Table browser. V nastaveniach tabuliek Group: Genes and Gene Predictions, Track: GENCODE v 32., zaklikneme Region: position, a Output fomat: CDS FASTA alignment a stlacime Get output
- Na dalsej obrazovke zaklikneme show nucleotides, zvolime MAF table multiz100way a vyberieme si, ktore organizmy chceme. V nasom pripade z primatov zvolime chimp, rhesus, bushbaby, z inych cicavcov mouse, rat, rabbit, pig, cow, dog, elephant a z dalsich organizmov opposum, platypus, chicken, stlacime Get output.
- Vystup ulozime do suboru, nechame si iba prvu formu genu (ENST00000308191.11_hg38), z mien sekvencii zmazeme spolocny zaciatok (ENST00000308191.11_hg38), pripadne celkovo prepiseme mena na anglicke nazvy

Skusme zostavit strom na stranke http://www.ebi.ac.uk/Tools/phylogeny/clustalw2_phylogeny/
- Distance correction: ako na prednáške, z počtu pozorovaných mutácií na evolučný čas
- Exclude gaps: vynechať všetky stĺpce s pomlčkami
- Clustering method: UPGMA predpokladá molekulárne hodiny, spájanie susedov nie
- P.I.M. vypíš aj maticu vzdialeností (% identity, pred korekciou)
- Vo vyslednom strome by sme mali zmenit zakorenenie, aby sme mali sliepku (chicken) ako outgroup

Výsledky z programu http://www.phylogeny.fr/alacarte.cgi , ktorý podporuje aj bootstrap:
- Vysledok s povodnym zakorenenim
- Vysledok so spravnym zakorenenim (chicken = outgroup)

"Spravny strom" [22] v nastaveniach Conservation track-u v UCSC browseri (podla clanku Murphy WJ et al Resolution of the early placental mammal radiation using Bayesian phylogenetics. Science 2001 [23])
Nas strom ma dost zlych hran: zle postavenie hlodavcov, ale aj slona a psa. Zle postavenie hlodavcov môže byť spôsobené long branch attraction.
Ak chcete skusit zostavit aj zarovnania, treba zacat z nezarovnanych sekvencii: [24]

Ďalšia zaujímavá webstránka s veľa nástrojmi https://usegalaxy.eu/

na tvorbu stromov sa dá použiť IQ-TREE
- modely vid tu: https://github.com/Cibiv/IQ-TREE/wiki/Substitution-Models
- vysledok #CB:phylo
viewer napr. http://phylotree.hyphy.org/

Gény, evolúcia a komparatívna genomika v UCSC genome browseri (cvičenie pri počítači)

Zobrazme si gén CLCA4 [25]
Zapnite si štandardnú sadu track-ov
Po kliknutí na gén si môžete prečítať o jeho funkcii, po kliknutí na ľavú lištu alebo na názov tracku v zozname na spodku stránky si môžete prečítať viac o tracku a meniť nastavenia
V tracku RefSeq genes si všimnite, že v tejto databáze má tento gén dve formy zostrihu, jedna z nich sa považuje za nekódujúcu, pretína sa aj s necharakterizovanou nekódujúcou RNA na opačnom vlákne
- Track RefSeq a jeho subtrack RefSeq Curated treba zapnut na pack
Nižšie vidíte track H3K27Ac Mark (Often Found Near Regulatory Elements) on 7 cell lines from ENCODE, kde bola táto histónová modifikácia v okolí génu detegovaná?
Všimnite si aj track DNase I Hypersensitivity, ktorý zobrazuje otvorený chromatin, prístupný pre viazanie transkripčných faktorov. Všimnite si jeho súvis s H3K27Ac trackom
Obidva tracky sú súčasťou tracku ENCODE regulation, v ktorom si môžete zapnúť aj ďalšie pod-tracky

Vsimnime si track Vertebrate Multiz Alignment & Conservation (100 Species)
- v spodnej casti tracku vidime zarovnania s roznymi inymi genomami
- v nastaveniach tracku zapnite Element Conservation (phastCons) na full a Conserved Elements na dense
- v tomto tracku vidíme PhyloP, co zobrazuje uroven konzerovanosti danej bazy len na zaklade jedneho stlpca zarovnania a dva vysledky z phyloHMM phastCons, ktory berie do uvahy aj okolite stlpce
Konkretne cast Conserved elements zobrazuje konkretne useky, ktore su najvac konzervovane
- Ak chceme zistit, kolko percent genomu tieto useky pokryvaju, ideme na modrej liste do casti Tools->Table browser, zvolime group Comparative genomics, track Conservation, table 100 Vert. El, region zvolime genome (v celom genome) a stlacime tlacidlo Summary/statistics, dostaneme nieco taketo:

item count	10,350,729
item bases	162,179,256 (5.32%)
item total	162,179,256 (5.32%)
smallest item	1
average item	16
biggest item	3,732
smallest score	186
average score	333
biggest score	1,000

- Ak by nas zaujimali iba velmi dlhe "conserved elements", v Table browser stlacime tlacidlo Filter a na dalsej obrazovke do policka Free-form query dame chromEnd-chromStart>=1500
- Potom mozeme skusit Summary/Statistics alebo vystup typu Hyperlinks to genome browser a Get output - dostaneme zoznam tychto elementov a kazdy si mozeme jednym klikom pozriet v browseri, napr. taketo
  - lod=24051 at chr1:50201403-50203312
  - lod=1899 at chr1:55663689-55667047 atd

Pozrime si teraz ten isty gen CLCA4 v starsej verzii genomu hg18 [26]
- 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 detegovany pozitivny vyber
- Po priblizeni do jedneho z exonov [27] vidite dosledky nesynonymnych mutacii

Poznamka: Existuju aj webservery na predikciu pozitivneho vyberu, napriklad tieto dva:

Selecton, clanok
Data monkey clanok
Skusili sme na Selecton poslat CLCA4 zo 7 cicavcov, subor tu: [28]
- vysledky [29] a [30] (metoda ale odporuca aspon 10 homologov)

Objavenie génu HAR1 pomocou komparatívnej 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 všetky regióny dĺžky aspoň 100bp s > 96% podobnosťou medzi šimpanzom a myšou/potkanom (35,000)
Porovnali s ostatnými cicavcami, zistili, ktoré majú veľa mutáci v človeku, ale málo inde (pravdepodobnostný model)
49 štatisticky významných regiónov, 96% nekódujúcich oblastiach
Najvýznamnejší HAR1: 118nt, 18 substitúcii u človeka, očakávali by sme 0.27. Iba 2 zmeny medzi šimpanzom a sliepkou (310 miliónov rokov), ale nebol nájdený v rybách a žabe.
Nezdá sa byť polymorfný u človeka
Prekrývajúce sa RNA gény HAR1A a HAR1B
HAR1A je exprimovaný v neokortexe u 7 a 9 týždenných embrií, neskôr aj v iných častiach mozgu (u človeka aj iných primátov)
Všetky substitúcie v človeku A/T->C/G, stabilnejšia RNA štruktúra (ale tiež sú blízko k telomére, kde je viacej takýchto mutácii kvôli rekombinácii a biased gene conversion)

Cvičenie pri počítači

Môžete si pozrieť tento region v browseri: chr20:63102114-63102274 (hg38), pricom ak sa este priblizite, uvidite zarovnanie aj s bazami a mozete vidiet, ze vela zmien je specifickych pre cloveka

Praktické príklady z predchádzajúcich cvičení

#CB01#QUAST:_program_na_.C5.A1tatistiky_o_kvalite_poskladania_gen.C3.B3mu Quast + Galaxy
#CB04#Pr.C3.ADklady_praktick.C3.BDch_programov lokálne vs globálne zarovnanie

Hľadanie génov

K hladaniu genov pozri aj prezentacie pdf

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 [31]
Glimmer at NCBI [32]
GeneMark at NCBI [33]
tRNAscan-SE [34]

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
Landscape of transcription in human cells. Djebali et al (ENCODE), Nature 2012

CI08

Felsensteinov algoritmus 1981

Mame dany strom T s dlzkami hran a bazy v listoch (jeden stlpec zarovnania) a model substitucii (zadany napr. maticou rychlosti R). Spocitajme pravdepodobnost, ze z modelu dostaneme prave tuto kombinaciu baz v listoch.
Oznacenie:
- 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, pricom koren je 2n-1.
- Nech p_v je rodic vrchola v a 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, vid minule cvicenia).
  - 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$
- Nech q_a je pravdepodobnost bazy a v koreni (ekvilibrium matice R)
  - Napr. v Jukes-Cantorovom modeli q_a = 1/4

Ak by sme poznali bazy vo vsetkych vrcholoch, mame $P(X_{1}=x_{1}\dots X_{{2n-1}}=x_{{2n-1}}|T,R)=q_{{x_{{2n-1}}}}\prod _{{v=1}}^{{2n-2}}P(x_{v}|x_{{p_{v}}},t_{v})$

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}}\in \{A,C,G,T\}^{{n-1}}}}P(X_{1}=x_{1}\dots X_{{2n-1}}=x_{{2n-1}}|T,R)$

Pocitat sucet cez exponencialne vela dosadeni hodnot za vnutorne vrcholy je neefektivne, spocitame rychlejsie dynamickym programovanim.
Nech A[v,a] je pravdepodobnost 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 y a z mame $A[v,a]=\sum _{{b,c}}A[y,b]A[z,c]P(b|a,t_{y})P(c|a,t_{z})$
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[y,b]P(b|a,t_{y}))(\sum _{c}A[z,c](c|a,t_{z}))$
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 (nerobili sme)

Co ak chceme spocitat pravdepodobnost P(X_v=a|X_1=x_1, X_2=x_2,\dots X_n=x_n,T,R)? Zaujimaju nas teda sekvencie genomov predkov.
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)$

CB08

Na týchto cvičeniach sa budeme venovať trom štatistickým témam súvisiacim s komparatívnou genomikou a s analýzou expresie génov. Tieto techniky sa však využívajú aj v iných oblastiach a môžete sa s nimi často stretnúť v genomických článkoch.

Zhlukovanie

Máme vstupné dáta, väčšinou ako vektory dĺžky n
Snažíme sa ich rozdeliť do skupín tak, aby dáta v rámci skupiny boli podobné a medzi skupinami rôzne

Využitie:

hľadanie génov s podobným profilom expresie
hľadanie skupín pacientov s podobným profilom expresie génov (objavovanie podtypov nejakej choroby)
hľadanie rodín podobných proteínov
automatická segmentácia obrázkov (napríklad rozlíšiť jednotlivé políčka microarray alebo gelu od pozadia)

Na prednáške sme videli hierarchické zhlukovanie, ktoré z dát vytvorilo strom. Teraz si ukážeme zhlukovanie, ktoré sa snaží dáta rozdeliť na k skupín, kde k je vopred daný parameter.

K-Means

Príklad vstupu pre zhlukovanie

Príklad 3 zhlukov a ich centier nájdených k-means algoritmom

pozri tiež prezentáciu pdf

Vstup: n-rozmerné vektory $x_{1},x_{2},...,x_{t}$ a počet zhlukov k

Výstup: Rozdelenie vektorov do k zhlukov takéto:

$c_{1},c_{2},...,c_{t}\;(1\leq c_{i}\leq k)$ - priradenie vektoru k zhluku

n-rozmerné vektory $\mu _{1},\mu _{2},...,\mu _{k}$ - centrá každého zhluku

Úloha: minimalizovať súčet štvorcov vzdialeností od každého vektoru k centru jeho zhluku:

$J(\mu ,c)=\sum _{{i=1}}^{t}{\big \|}x_{i}-\mu _{{c_{i}}}{\big \|}_{2}^{2}$

${\big \|}x_{i}-\mu _{{c_{i}}}{\big \|}_{2}^{2}$ je druhá mocnina vzdialenosti vektora xi od centra jeho zhluku

Algoritmus

Heuristika, ktorá nenájde vždy najlepšie zhlukovanie. Začne z nejakého zhlukovania a postupne ho zlepšuje. Pozri aj clanok na Wikipedii

inicializácia: náhodne vyber k centier $\ \mu _{1},\mu _{2},...,\mu _{k}$
opakuj kým sa niečo mení:

priraď každý bod najbližšiemu centru: $c_{i}=\arg \min _{j}{\big \|}x_{i}-\mu _{j}{\big \|}_{2}$

vypočítaj nové centroidy: $\mu _{j}=\operatorname {avg_{{i:c_{i}=j}}}x_{i}$ (spriemerujeme všetky body v jednom zhluku)

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.
Kontingencna tabulka

	Pozitivny vyber	Bez poz. vyberu	Sucet
Imunitne	8 (nip)	62	70 (ni)
Ostatne	392	16067	16459
Sucet	400 (np)	16129	16529 (n)

Nulova hypoteza: geny v nasom zozname boli nahodne vybrane z celeho genomu, t.j. ak v celom genome je frekvencia imunitnych genov ni/n (cca 0.4%), vo vzorke velkosti np (geny s pozitivnym vyberom) ocakavame cca np * (ni / n) imunitnych genov.
- aj v nulovej hypoteze vsak vzorka velkosti ni cisto nahodou moze obsahovat viac alebo menej takych genov.
- presnejsie mame urnu so ni (70) bielymi a n-ni (16459) ciernymi gulickami, vytiahneme nahodne np (400) guliciek, kolko bude medzi nimi bielych, nazvime tuto nahodnu premennu Xip
- v nasom priklade by sme ocakavali 1.7 genu s innate immune response, ale mame 8 (4.7xviac)
Rozdelenie pravdepodobnosti Xip je hypergeometricke, t.j. $\Pr(X_{{ip}}=n_{{ip}})={n_{i} \choose n_{{ip}}}{n-n_{i} \choose n_{p}-n_{{ip}}}/{n \choose n_{p}}$
Aka je pravdepodobnost, ze v nulovej hypoteze bude Xip 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

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.
- Reimand, Jüri, et al. "Pathway enrichment analysis and visualization of omics data using g: Profiler, GSEA, Cytoscape and EnrichmentMap." Nature protocols 14.2 (2019): 482. [35]

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);
# pr pocet bielych>7 (t.j. >=8) ak taham 400 z vreca so 70 bielymi a 16529-70 ciernymi
# sucet cisiel z tabulky od 2.332580e-04 az po koniec

d = dhyper(0:15, 70, 16529-70, 400);
plot(0:15,d)

# test pre danu tabulku
a=matrix(c(8,62,392,16067),nrow=2, ncol=2)
fisher.test(a,alternative = "greater")

Multiple testing correction

V mnohych situaciach robime vela testov toho isteho typu, kazdy ma urcitu p-value
Napr. testujeme 1000 genov v genome na pozitivny vyber, zvolime tie, kde p-value <= 0.05
Alebo testujeme obohatenie 1000 funkcnych kategorii v nejakej vzorke genov, zvolime tie, kde p-value <= 0.05
Problem: ak kazda z 1000 kategorii ma 5% sancu tam byt len nahodou, ocakavali by sme 50 cisto nahodnych pozitivnych vysledkov. Ak sme napr. nasli 100 pozitivnych vysledkov (obohatenych kategorii), cca polovica z nich je zle
Preto potrebujeme pri velkom mnozstve testov umelo znizit prah na p-value tak, aby nahodny sum netvoril velke percento nasich vysledkov
Toto sa vola multiple testing correction, je viac technik, napr. FDR (false discovery rate)

Nadreprezentácia, Uniprot (cvičenie pri počítači)

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

Chodte na genome browser http://genome-euro.ucsc.edu/
Zvolte Tools->Gene Sorter, sort by nechajme Expression (GTEx), a do okienka search zadajme identifikator genu PTPRZ1
- Dostane tabulku genov s podobnym profilom expresie ako PTPRZ1 (červená je vysoká expresia, zelená nízka)
- Zoznam tychto genov v textovom formate najdete tu
http://biit.cs.ut.ee/gprofiler/ mena genov skopirujme do policka Query, stlacte g:Profile!
- Ak by výpočet dlho trval, nájdete ho aj tu
- Vo výslednej tabuľke je každý riadok jedna funkcna kategoria, v ktorej su geny s tymto profilom expresie nadreprezentovane, kazdy stlpec jeden gen.
- V spodnej casti tabuly su aj asociacie k chorobam a k transkripcnym faktorom, ktore by mohli prislusne geny regulovat
Co by sme na zaklade nadreprezentovanych kategorii usudzovali o gene PTPRZ1?

Najdite tento gen v Uniprote (http://www.uniprot.org/), potvrdzuje nase domnienky?
- O mnohých údajoch na stránke sme sa rozprávali na prednáške (GO kategórie, domény, sekundárna a 3D štruktúra)
- na veľa miestach na stránke je uvedené aj odkiaľ jednotlivé údaje pochádzajú
- Všimnime si Pfam domény a pozrime si ich stránku

Vratme sa do genome browsera, najdime si PTPRZ1 gen v genome [36]
V browseri su rozne tracky tykajuce sa expresie, napr. GTEx. 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 uvidite zase prehlad expresie v roznych tkanivach (podla GTEx)

CI09

Hľadanie motívov zadefinovaných pravdepodobnostnou maticou

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

EM algoritmus

Iterativne zlepsuje W, pricom berie vsetky O vahovane podla ich pravdepodobnosti vzhladom na W z minuleho kola
Videli sme na prednaske, tu je trochu prepisany:

Inicializácia: priraď každej pozícii j v sekvencii $S_{i}$ nejaké skóre $p_{{i,j}}$
Iteruj:
- Spočítaj W zo všetkých možných výskytov v $S_{1},\dots ,S_{k}$ váhovaných podľa $p_{{i,j}}$
- Prepočítaj všetky skóre $p_{{i,j}}$ tak, aby zodpovedali pomerom pravdepodobností výskytu W na pozícii j v $S_{i}$ , t.j. $p_{{i,j}}$ je umerne $Pr(S_{i}|W,o_{i}=j)$ , pricom hodnoty normalizujeme tak, aby sucet v riadku bol 1

Gibbsovo vzorkovanie (Gibbs sampling)

Inicializácia: Vezmi náhodné pozície výskytov O
Iteruj:
- Spočítaj W z výskytov O
- Vyber náhodne jednu sekvenciu $S_{i}$
- Pre každú možnú pozíciu j v $S_{i}$ spočítaj skóre $p_{{i,j}}$ (ako v EM) výskytu W na tejto pozícii
- Zvoľ $o_{i}$ náhodne s váhovaním podľa $s_{{i,j}}$

Takto dostavame postupnost vzoriek $O^{{(0)}},O^{{(1)}},...$ .
Za sebou iduce vzorky sa podobaju (lisia sa len v jednej zlozke $o_{i}$ ) nie su teda nezavisle
Pre kazdu vzorku $O^{{(t)}}$ najdeme najlepsie $W^{{(t)}}$ a spocitame vierohodnost $\Pr(S|W^{{(t)}},O^{{(t)}})$ . Nakoniec vyberieme O a W, kde bola vierohodnost najvyssia.
Tento algoritmus (s malymi obmenami) bol pouzity v clanku Lawrence, Charles E., et al. (1993) "Detecting subtle sequence signals: a Gibbs sampling strategy for multiple alignment." Science.
- V clanku v kazdej iteracii maticu W rataju zo vsetkych sekvencii okrem $S_{i}$
- Obcas robia krok, kde nahodne skusaju posunut vsetky vyskyty o jedna dolava alebo doprava
- Tento algoritmus nie je uplne matematicky korektne Gibbsovo vzorkovanie (nema ani poradne zadefinovane rozdelenie, z ktoreho vzorkuje). Na spodku stranky pre informaciu uvadzame algoritmus Gibbsovho vzorkovanie pre hladanie motivov z ineho clanku.

Vzorkovanie z pravdepodobnostného modelu vo všeobecnosti

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

Pouzitie vzoriek

spomedzi ziskanych vzoriek zvolime tu, pre ktoru je vierohodnost Pr(D|X) najvacsia (iny pristup k maximalizovaniu vierohodnosti)
ale vzorky nam daju aj informaciu o tom, aka je velka neurcitost v odhade X
- mozeme odhadovat stredne hodnoty a odchylky roznych velicin
- napr. pri hladani motivov mozeme sledovat ako casto je ktora pozicia vyskytom motivu

generovat nezavisle vzorky z Pr(X|D) moze byt tazke
metoda Markov chain Monte Carlo (MCMC) generuje postupnost zavislych vzoriek $X^{{(0)}},X^{{(1)}},\dots$ , konverguje v limite k cielovej distribucii Pr(X|D)
Gibbsovo vzorkovanie je specialnym pripadom MCMC

Markovove reťazce

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

Matica

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

Stacionarne rozdelenie

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

Priklady Markovovskych retazcov v bioinformatike

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

Ergodické Markovove reťazce

Vravime ze matica je ergodicka, ak $P^{t}$ pre nejake t>0 ma vsetky polozky nenulove
Priklady neergodickych matic

1 0          0.5 0.5          0 1             0.5 0.5
0 1          0   1            1 0             1   0
nesuvisla    slabo suvisla    periodicka      ergodicka

V HMM stavy tvoria Markovov retazec; hladanie genov ergodicky stavovy priestor, profilove HMM nie

Markov chain Monte Carlo MCMC

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

Gibbsovo vzorkovanie

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

Dôkaz správnosti Gibbsovho vzorkovania

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

Poriadnejšie Gibbsovo vzorkovanie pre motívy

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

Pravdepodobnostny model

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

Gibbsovo vzorkovanie

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

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

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

Dalsie detaily vypoctu $\Pr(S|O)$ :

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

CB09

Bezkontextové gramatiky

Na modelovanie štruktúry RNA sa používajú stochastické bezkontextové gramatiky (bude na ďalšej prednáške)
My si teraz ukážeme bezkontextové gramatiky, ktoré nemajú pravdepodobnosti
Zaviedol Noam Chomsky v lingvistike 50-te roky 20. storočia, tiež dôležité v informatike

Gramatika

Príklad: S->aSb, S->epsilon (píšeme aj skrátene S->aSb|epsilon)
Dva typy symbolov: terminály (malé písmená), neterminály (veľké písmená)
Pravidlá prepisujúce neterminál na reťazec terminálov a neterminálov (môže byť aj prázdny reťazec, ktorý označujeme epsilon)
Neterminál S je "štartovací"

Použitie gramatiky na generovanie reťazcov

Začneme so štartovacím neterminálom S
V každom kroku prepíšeme najľavejší neterminál podľa niektorého pravidla
Skončíme, keď nezostanú žiadne neterminály
Príklad: S->aSb->aaSbb->aaaSbbb->epsilon
Aké všetky slová vie táto gramatika generovať?
- V tvare aa...abb...b s rovnakým počtom á-čok a b-čiek (informatici píšu $a^{k}b^{k}$ )

Cvičenia

Zostavte gramatiku na slová typu aa..abb..b kde á-čok je rovnako alebo viac ako b-čok, pre
- S->aSb|aS|epsilon
Zostavte gramatiku pre slová toho istého typu, kde á-čok je viac ako b-čok, t.j. i>j
- S->aSb|aT T->aT|epsilon (alebo S->aSb|aS|a)
Zostavte gramatiku pre dobre uzátvorkované výrazy zo zátvoriek (,),[,]. Napr. [()()([])] je dobre uzátvorkovaný, ale [(]) nie je.
- S->SS|(S)|[S]|epsilon
- príklad odvodenia v tejto gramatike: S->[S]->[SS]->[SSS]->[(S)SS]->[()SS]->[()(S)S]->[()()S]->[()()(S)]->[()()([S])]->[()()([])]

Parsovanie reťazca pomocou gramatiky: určiť, ako mohol byt reťazec vygenerovaný pomocou pravidiel

Gramatika pre dobre uzátvorkované výrazy nám pomôže určiť, ktorá zátvorka patrí ku ktorej: tie, ktoré boli vygenerované v jednom kroku

Ďalšie cvičenia

Zostavte gramatiku na DNA palindromy, t.j. sekvencie, ktore zozadu po skomplementovani baz daju to iste, ako napr. GATC
- S->gSc|cSg|aSt|tSa|epsilon
Vlasenky RNA s lubovolne dlhou sparovanou castou a 3 nesparovanymi nukleotidmi v strede
- S->gSc|cSg|aSu|uSa|aaa|aac|aag|aau|...|uuu

Tazsi priklad: Zostavte gramatiku na slova s rovnakym poctom acok a bcok v lubovolnom poradi
- S->epsilon|aSbS|bSaS
- ako bude generovat aababbba?
- preco vie vygenerovat vsetky take retazce?

CI10

Ú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 roznych databazach
- Pubmed: databaza clankov, napr. najdime phastcons
- Gene: najdime DNA polymerase
- 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, % zhodnych 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 známy koronavírusový proteín Spike
- Nájdime ho na stránke http://www.uniprot.org/ pod názvom SPIKE_SARS2
- Pozrime si podrobne jeho stránku, ktoré časti boli predpovedané bioinformatickými metódami z prednášky?
- Všimnime si niektorú Pfam doménu a pozrime si jej stránku

UCSC genome browser

http://genome-euro.ucsc.edu/
On-line grafický nástroj na prezeranie genómov
Konfigurovateľný, veľa možností, prijemne pouzivatelske rozhranie
Moznost stiahnut data vhodne na dalsie spracovanie alebo zobrazit vlastne data
Pomerne málo organizmov
- doraz hlavne na ludsky genom

Základy

Adresa http://genome-euro.ucsc.edu/
Hore v modrom menu zvoľte Genomes, potom zvoľte ľudský genóm. Do okienka search term zadajte HOXA2. Vo výsledkoch hľadania (UCSC genes) zvoľte gén homeobox A2 na chromozóme 7.
- Pozrime si spolu túto stránku
- V hornej časti sú ovládacie prvky na pohyb vľavo, vpravo, približovanie, vzďaľovanie
- Pod tým schéma chromozómu, červeným vyznačená zobrazená oblasť
- Pod tým obrázok vybranej oblasti, rôzne tracky
- Pod tým zoznam všetkých trackov, dajú sa zapínať, vypínať a konfigurovať
- Po kliknutí na obrázok sa často zobrazí ďalšia informácia o danom géne alebo inom zdroji dát
- V génoch exony hrubé, UTR tenšie, intróny vodorovné čiary
Po kliknutí na gén alebo inú časť nejakého tracku väčšinou o ňom dostaneme viac informácií. Kliknutim na listu ku tracku (lavy okraj obazku) sa dozviete viac o tracku a mozete nastavovat parametre zobrazenia

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 (hg19 a hg38)
Na tej istej stránke dole nájdete stručný popis zvolenej verzie genómu. Pre ktoré oblasti genómu máme v hg38 najviac alternatívnych verzií?
Zadajte región chr21:31,250,000-31,300,000 v hg19 [37]
Zapnite si tracky Mapability a RepeatMasker 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
Ako a prečo sa pri rôznych dĺžkach readov líšia? (Keď kliknete na linku "Mapability", môžete si prečítať bližšie detaily.)
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 v hg19. [38] 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ĺžka aj tak nebola známa, nakoľko nebola osekvenovaná).
Cez menu položku View, In other genomes si pozrite, ako zobrazený úsek vyzerá vo verzii hg38. Ako sa zmenila dĺžka z pôvodných 300kb?
Prejdite na genóm Rhesus, verzia rheMac2, región chr7:59,022,000-59,024,000 [39], zapnite si tracky Contigs, Gaps, Quality scores
Aké typy problémov v kvalite sekvencie v tomto regióne vidíte?

Komparativna genomika

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 na uroven "base", t.j. zobrazenych cca 100bp, v obdlzniku multiz alignment uvidite zarovnanie s homologickym usekom v inych genomoch.
V casti conservation by PhyloP vidime graf toho, ako silne su zachovane jednotlive stlpce zarovnania
Da sa zapnut track Placental Chain/Net a pozriet sa na ktorych chromozomoch je ortologicky usek v inych genomoch

Blat

Choďte na UCSC genome browser (http//genome.ucsc.edu/), na modrej lište zvoľte BLAT, zadajte DNA sekvenciu vyssie 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*

Úvod do bezkontextových gramatík

Na modelovanie štruktúry RNA sa používajú stochastické bezkontextové gramatiky (bude na ďalšej prednáške)
Tie sú založené na bezkontextových gramatikách, ktoré mnohí poznáte z bakalárskeho štúdia

Gramatika

Príklad: S->aSb, S->epsilon (píšeme aj skrátene S->aSb|epsilon)
Dva typy symbolov: terminály (malé písmená), neterminály (veľké písmená)
Pravidlá prepisujúce neterminál na reťazec terminálov a neterminálov (môže byť aj prázdny reťazec, ktorý označujeme epsilon)
Neterminál S je "štartovací"

Použitie gramatiky na generovanie reťazcov

Začneme so štartovacím neterminálom S
V každom kroku prepíšeme najľavejší neterminál podľa niektorého pravidla
Skončíme, keď nezostanú žiadne neterminály
Príklad: S->aSb->aaSbb->aaaSbbb->epsilon
Aké všetky slová vie táto gramatika generovať?
- V tvare aa...abb...b s rovnakým počtom á-čok a b-čiek (informatici píšu $a^{k}b^{k}$ )

Cvičenia

Zostavte gramatiku na slova typu aa..abb..b kde acok je rovnako alebo viac ako bcok, pre
- S->aSb|aS|epsilon
Zostavte gramatiku pre slova toho isteho typu, kde acok je viac ako bcok, t.j. i>j
- S->aSb|aT T->aT|epsilon (alebo S->aSb|aS|a)
Zostavte gramatiku pre dobre uzatvorkovane vyrazy zo zatvoriek (,),[,]. Napr. [()()([])] je dobre uzatvorkovany, ale [(]) nie je.
- S->SS|(S)|[S]|epsilon
- priklad odvedenia v tejto gramatike: S->[S]->[SS]->[SSS]->[(S)SS]->[()SS]->[()(S)S]->[()()S]->[()()(S)]->[()()([S])]->[()()([])]

Parsovanie retazca pomocou gramatiky: urcit, ako mohol byt retazec vygenerovany pomocou pravidiel

Gramatika pre dobre uzatvrokovane vyrazy nam pomoze urcit, ktora zatvorka patri ku ktorej: tie, ktore boli vygenerovane v jednom kroku

Dalsie cvicenia

Zostavte gramatiku na DNA palindromy, t.j. sekvencie, ktore zozadu po skomplementovani baz daju to iste, ako napr. GATC
- S->gSc|cSg|aSt|tSa|epsilon
Vlasenky RNA s lubovolne dlhou sparovanou castou a 3 nesparovanymi nukleotidmi v strede
- S->gSc|cSg|aSu|uSa|aaa|aac|aag|aau|...|uuu

Tazsi priklad: Zostavte gramatiku na slova s rovnakym poctom acok a bcok v lubovolnom poradi
- S->epsilon|aSbS|bSaS
- ako bude generovat aababbba?
- preco vie vygenerovat vsetky take retazce?

CB10

RNA štruktúra

Znama databaza rodin RNA genov je Rfam: http://rfam.xfam.org/
Najdite si v nej rodinu RF00015 (U4 spliceosomal RNA)
V casti Secondary structure si mozete pozriet obrazky farebne kodovane podla roznych kriterii
- Skuste pochopit, co jednotlive obrazky a ich farby znamenaju
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 GENCODE genes, conservation, RepeatMasker v jej okoli
Vo verzii hg19 (kam sa viete z inej verzii dostat cez horne menu View->In Other Genomes) je track "CSHL Sm RNA-seq" ktory obsahuje RNASeq kratkych RNA z roznych casti buniek, zapnite si v jeho nastaveniach aj zobrazenie RNA z jadra (nucleus)
Zadajte sekvenciu na RNAfold serveri [40]
Ak vypocet dlho trva, pozrite si vysledok tu
Podoba sa na strukturu zobrazenu v Rfame? v com sa lisi?

RNA dizajn: mozete sa skusit zahrat na stranke http://www.eternagame.org/web/

PSI BLAST

Toto cvičenie je z časti inšpirované stránkou [41]
Budeme uvažovať vzdialene podobné enzýmy
- Bis(5'-adenosyl)-triphosphatase (Uniprot)
- Galactose-1-phosphate uridylyltransferase (GALT/GAL7) (Uniprot)
- Ich domény patria v databáze Pfam do toho istého klanu
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, spustime program PSI-BLAST
V prvom kole PSI-BLAST spúšťa bežný BLASTP
GAL gén (konkrétne GAL7_HAEIN, accession P31764) sa nachádza medzi výsledkami, ale má príliš vysokú E-value
Spustíme teraz druhú iteráciu PSI-BLAST, ktorá zostaví profil z proteínov s nízkou E-value v prvej iterácii
Aká je E-value nájdeného zarovnania?
Ak by výpočet dlho trval, výsledky sú tu: 1. kolo, 2. kolo

Nadreprezentácia, Uniprot (cvičenie pri počítači)

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

Chodte na genome browser http://genome-euro.ucsc.edu/
Zvolte Tools->Gene Sorter, sort by nechajme Expression (GTEx), a do okienka search zadajme identifikator genu PTPRZ1
- Dostane tabulku genov s podobnym profilom expresie ako PTPRZ1 (červená je vysoká expresia, zelená nízka)
- Zoznam tychto genov v textovom formate najdete tu
http://biit.cs.ut.ee/gprofiler/ mena genov skopirujme do policka Query, stlacte g:Profile!
- Ak by výpočet dlho trval, nájdete ho aj tu
- Vo výslednej tabuľke je každý riadok jedna funkcna kategoria, v ktorej su geny s tymto profilom expresie nadreprezentovane, kazdy stlpec jeden gen.
- V spodnej casti tabuly su aj asociacie k chorobam a k transkripcnym faktorom, ktore by mohli prislusne geny regulovat
Co by sme na zaklade nadreprezentovanych kategorii usudzovali o gene PTPRZ1?

Najdite tento gen v Uniprote (http://www.uniprot.org/), potvrdzuje nase domnienky?
- O mnohých údajoch na stránke sme sa rozprávali na prednáške (GO kategórie, domény, sekundárna a 3D štruktúra)
- na veľa miestach na stránke je uvedené aj odkiaľ jednotlivé údaje pochádzajú
- Všimnime si Pfam domény a pozrime si ich stránku

Vratme sa do genome browsera, najdime si PTPRZ1 gen v genome [42]
V browseri su rozne tracky tykajuce sa expresie, napr. GTEx. 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 uvidite zase prehlad expresie v roznych tkanivach (podla GTEx)

Objavenie génu HAR1 pomocou komparatívnej 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 všetky regióny dĺžky aspoň 100bp s > 96% podobnosťou medzi šimpanzom a myšou/potkanom (35,000)
Porovnali s ostatnými cicavcami, zistili, ktoré majú veľa mutáci v človeku, ale málo inde (pravdepodobnostný model)
49 štatisticky významných regiónov, 96% nekódujúcich oblastiach
Najvýznamnejší HAR1: 118nt, 18 substitúcii u človeka, očakávali by sme 0.27. Iba 2 zmeny medzi šimpanzom a sliepkou (310 miliónov rokov), ale nebol nájdený v rybách a žabe.
Nezdá sa byť polymorfný u človeka
Prekrývajúce sa RNA gény HAR1A a HAR1B
HAR1A je exprimovaný v neokortexe u 7 a 9 týždenných embrií, neskôr aj v iných častiach mozgu (u človeka aj iných primátov)
Všetky substitúcie v človeku A/T->C/G, stabilnejšia RNA štruktúra (ale tiež sú blízko k telomére, kde je viacej takýchto mutácii kvôli rekombinácii a biased gene conversion)

Cvičenie pri počítači

Môžete si pozrieť tento region v browseri: chr20:63102114-63102274 (hg38), pricom ak sa este priblizite, uvidite zarovnanie aj s bazami a mozete vidiet, ze vela zmien je specifickych pre cloveka

Expresia génov

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

Databaza gene expression dat na NCBI
Do Search okienka zadajme GDS2925
Mali by sme dostat dataset 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
- napr. K=4 a K=5 pre Pearsonovu korelaciu
- mozeme is pozriet aj hierarchicke zhlukovanie

Sekvenčné motívy, program MEME

Vazobne miesta transkripcnych faktorov sa casto reprezentuju ako sekvencne motivy
Ak mame skupinu sekvencii, mozeme hladat motiv, ktory maju spolocny
Znamy program na tento problem je MEME
Chodte na stranku http://meme-suite.org/
Zvolte nastroj MEME a v casti Input the primary sequences zvolte Type in sequences a zadajte tieto sekvencie
Pozrite si ostatne nastavenia. Co asi robia?
Ak server pocita dlho, mozete si pozriet vysledky tu

Kvasinkové transkripčné faktory v SGD

Yeast genome database SGD obsahuje pomerne podrobne stranky pre jednotlive transkripcne faktory
Pozrime si stranku pre transkripcny faktor GAL4 [43]

Uniprot

Prehladnejsi pohlad na proteiny, vela linkov na ine databazy, cast vytvarana rucne
- Pozrieme sa na známy koronavírusový proteín Spike
- Nájdime ho na stránke http://www.uniprot.org/ pod názvom SPIKE_SARS2
- Pozrime si podrobne jeho stránku, ktoré časti boli predpovedané bioinformatickými metódami z prednášky?
- Všimnime si niektorú Pfam doménu a pozrime si jej stránku

Nussinovovej algoritmus (nerobili sme)

Z cvičných príkladov na skúšku

Vyplňte maticu dynamického programovania (Nussinovovej algoritmus) pre nájdenie najväčšieho počtu dobre uzátvorkovaných spárovaných báz v RNA sekvencii GAACUUCACUGA (dovoľujeme len komplementárne páry A-U, C-G) a nakreslite sekundárnu štruktúru, ktorú algoritmus našiel.

 G A A C U U C A C U G A
 0 0 0 1 1 2 3 3 3 4 4 4  G
   0 0 0 1 2 2 2 2 3 4 4  A
     0 0 1 1 1 2 2 2 3 4  A
       0 0 0 0 1 1 1 2 3  C
         0 0 0 1 1 1 2 3  U
           0 0 1 1 1 2 3  U
             0 0 0 1 2 2  C
               0 0 1 1 1  A
                 0 0 1 1  C
                   0 0 1  U
                     0 0  G
                       0  A

Ako by sme algoritmus upravili, aby dlzka slucky na konci helixu bola vzdy aspon 3?

CI11

RNA struktura

Opakovanie Nussinovovej algoritmu

Z cvičných príkladov na skúšku

Vyplnte maticu dynamického programovania (Nussinovovej algoritmus) pre nájdenie najväčšieho počtu dobre uzátvorkovaných spárovaných báz v RNA sekvencii GAACUUCACUGA (dovoľujeme len komplementárne páry A-U, C-G) a nakreslite sekundárnu štruktúru, ktorú algoritmus našiel.

Rozsirenia Nussinovovej algoritmu

lahke: kazdy par i,j musi mat vzdialenost |i-j|>=3 (RNA sa na kratsom useku nevie ohnut o 180 stupnov)
tazsie (bolo s hintom na skuske): chceme davat skore iba "stackovanym parom", t.j. ak i a j aj i+1 a j-1 su sparovane, dostaneme +1, osamotene pary nedostavaju ziadne skore. Úlohou je opäť pre danú sekvenciu nájsť dobre uzátvorkovanú štruktúru s maximálnym skóre.
- pomocka: pouzijeme dve tabulky A a B, pričom A[i,j] obsahuje maximálne skóre pre podreťazec X[i...j] a B[i...j] obsahuje maximálne skóre pre podreťazec X[i...j], za predpokladu, že X[i] a X[j] sú spárované v štruktúre (táto hodnota je definovaná iba pre i a j, kde sú X[i] a X[j] komplementárne).

Stochasticke bezkontextove gramatiky

Ako asi funguje algoritmus, ktory hlada najpravdepodobnejsie odvodenie?
- rozsirme Nussinovovej algoritmus o dalsi rozmer - neterminal, z ktoreho je podretazec X[i...j] vygenerovany
Je najpravdepodobnejsie odvodenie to iste ako najpravdepodobnejsia sekundarna struktura pri gramatike z prednasky?
- S->aSu|uSa|cSg|gSc|aS|cS|gS|uS|Sa|Sc|Sg|Su|SS|epsilon
- jednu strukturu vieme vyjadrit pomocou viacerych odvodeni, napr. v jednoduchej strukture nizsie vieme slucku ccu generovat zlava aj sprava (cS vs Su), tiez hocikde vieme spravit S->SS a potom jedno S znicit

acgccucgu
(((...)))

Viete zmenit gramatiku tak, aby najlavejsie odvodenia zodpovedali 1 k 1 sekundarnym strukturam?
- napr. S->aS|cS|gS|tS|aSuS|uSaS|cSgS|gScS|epsilon
- vid clanok Dowell RD, Eddy SR. Evaluation of several lightweight stochastic context-free grammars for RNA secondary structure prediction. BMC bioinformatics. 2004 Jun 4;5(1):71. [44]

CB11

Ukážka práce v Linuxe

Prvá časť - príprava

Prihláste sa na server podľa pokynov.
Potom spúšťajte jednotlivé príkazy podľa pokynov nižšie.
Odporúčame príkazy kopírovať myšou (v internetovom prehliadači vysvietiť, stlačiť Ctrl-C, v konzole Ctrl-Shift-V)

# riadky začínajúce mrežou # sú komentáre, netreba ich spúšťať

# Dôležité: v príkazoch nižšie xx nahraďte vašim číslom skupiny, napr. 01
mkdir xx
cd xx
# príkaz mkdir (make directory) vytvoril priečinok
# príkaz cd (change directory) zmenil váš aktuálny priečinok na tento nový

# v konzole by ste mali mať user@server:~/xx$
# kde xx je číslo vašej skupiny, napr. 01

# stiahneme si súbor s dátami zo stránky
wget http://compbio.fmph.uniba.sk/vyuka/mbi-data/cb12.zip
# rozzipujeme ho
unzip cb12.zip

Druhá časť - skladanie genómov, mapovanie čítaní, zarovnanie

# prejdeme na priečinok s prvou časťou ohľadom sekvenovania
cd 1-seq

# ls vypíše zoznam súborov v priečinku
ls
# ls -l vypíše dlhšiu informáciu (long)
ls -l
# ls -lSh usporiada súbory podľa veľkosti (Size) a veľkosti vypíše priateľskejšie pre ľudí (human)
ls -lSh

# mali by sme vidieť kúsok sekvencie z E.coli (prípona .fasta)
# a 2 súbory zo sekvenovania prístrojom Illumina Miseq  (prípona .fastaq.gz)
# tieto súbory obsahujú čítania z vyššie uvedeného kúsku genómu


# ideme skladať genóm, bude to trvať dlho, preto to chceme spustiť na pozadí
# aby sme mohli medzitým robiť niečo iné
screen # stlačte Enter
# spustite skladanie programom spades
spades.py -t 1 -m 1 --pe1-1 miseq_R1.fastq.gz --pe1-2 miseq_R2.fastq.gz -o spades > spades.log
# stlačte naraz Ctrl-a potom d
# spades teraz beží na pozadí

# príkaz top zobrazí bežiace procesy
# ukončíte ho stlačením q (quit)
top

# príkaz less umožňuje prezerať si obsah textového súboru
# aj príkaz less ukončíte stlačením q, šípkami sa pohybujete po súbore
less ref.fasta
# čítania sú komprimované, preto namiesto less použijeme zless
zless miseq_R1.fastq.gz
# tieto príkazy spočítajú počet riadkov - ako z toho zistíme počet čítaní?
zcat miseq_R1.fastq.gz | wc -l 
zcat miseq_R2.fastq.gz | wc -l 

# keď spades skončí, vrátime sa do screen a ukončíme ho
screen -r
# exit ukončí screen
exit

# spades dal výstup do podpriečinku spades, pozrime si ho
ls spades
# skopírujeme si hlavný výsledok do nášho priečinka (cp = copy)
cp -ip spades/contigs.fasta spades.fasta
less spades.fasta
# pozrime si hlavičky jednotlivých sekvencií vo fasta súbore
grep '>' spades.fasta

# programom last si spravíme dotplot referencia vs. naše skladanie
# 1) vytvorenie indexu pre referenciu
lastdb ref.fasta ref.fasta 
# 2) samotné zarovnanie
lastal -f TAB ref.fasta spades.fasta > aln.tab
# 3) vytvorenie obrázku s dotplotom
last-dotplot aln.tab aln.png

# a ešte dotplot referencia vs. referencia
# 2) samotné zarovnanie (index už máme)
lastal -f TAB ref.fasta ref.fasta > aln2.tab
# 3) vytvorenie obrázku s dotplotom
last-dotplot aln2.tab aln2.png

# pozrieme si dotploty programom eog
eog aln.png &
eog aln2.png &


# zarovnajme čítania k referenčnému genómu v 4 krokoch
# 1) indexovanie fasta súboru
bwa index ref.fasta
# 2) samotné zarovnávanie čítaní programom bwa
bwa mem ref.fasta miseq_R1.fastq.gz miseq_R2.fastq.gz > ref-miseq.sam
# 3) zmeníme textový sam formát na binárny bam formát
samtools view -S -b ref-miseq.sam | samtools sort - -o ref-miseq.bam
# 4) vytvoríme index bam súboru
samtools index ref-miseq.bam

# pozrime sa na zoznam súborov od najnovšieho po najstarší
ls -lth
# sam súbor so zarovnaniami sa dá pozrieť, ale nie je veľmi prehľadný
less ref-miseq.sam


# vytvoríme aj zarovnanie nášho poskladaného genómu k referencii vo formáte bam
samtools faidx ref.fasta
lastal ref.fasta spades.fasta -E1e-20 | maf-convert sam > ref-spades.sam
samtools view -S -b -t ref.fasta.fai ref-spades.sam | samtools sort - -o ref-spades.bam
samtools index ref-spades.bam

# výsledky si zobrazíme v grafickom prehliadači igv 
# obdoba genome browsera, ktorú si môžete nainštalovať na vašom počítači
# POZOR: POTREBUJE VEĽA PAMÄTE, SPUSTÍME IBA JEDEN NARAZ
igv -g ref.fasta
# pomocou Menu->File->Load from File otvorte ref-spades.bam a ref-miseq.bam
# pozrime si región ecoli-frag:224,000-244,000
#   Vidíte jednotlivé kontigy? Sedí tento pohľad s dotplotom? 
# a potom bližšie ecoli-frag:227,300-227,600
#   Všimnite si sekvenačné chyby rozdiely medzi referenciou a kontigmi

Tretia časť - hľadanie génov, RNA-seq

# v druhom cvičení si vyskúšame hľadanie génov
# najskôr sa presuňme do druhého priečinku
cd ../2-genes

# pozrime si, aké máme súbory
ls -lSh
# mali by sme mať kúsok referenčného genómu huby Aspergillus nidulans 
# fastq súbor s čítaniami z RNA-seq pre tento kúsok referencie
# gff súbor s anotáciou génov z databázy

# spustíme hľadač génov Augustus 2x:
# raz s parametrami priamo pre A.nidulans a raz s parametrami pre ľudský genóm
augustus --species=anidulans ref2.fasta > augustus-anidulans.gtf
augustus --species=human ref2.fasta > augustus-human.gtf

# RNA-seq zarovnáme k sekvencii nástrojom tophat2 (podporuje intróny)
bowtie2-build ref2.fasta ref2.fasta
tophat2 -i 10 -I 10000 --max-multihits 1 --output-dir rnaseq ref2.fasta rnaseq.fastq
samtools sort rnaseq/accepted_hits.bam rnaseq
samtools index rnaseq.bam

# predikcie génov a RNA-seq si pozrieme v igv
igv -g ref2.fasta
# v igv si otvorte annot.gff, augustus-anidulans.gtf, augustus-human.gtf, rnaseq.bam
# - ktoré parametre Augustusu dali presnejšie predpovede (za predpokladu, že anotácia je správna)
# - pozrite si zblízka niektorý gén s vysokou expresiou (napr. druhy gen sprava), 
#   mali by ste vidieť čítania podporujúce intróny

CI12

Protein threading

Praktické programy na NP ťažké problémy

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)

najdu optimalne riesenie, mnohe instancie zrataju v rozumnom case, ale mozu bezat aj velmi dlho
CPLEX [45] a Gurobi [46] komercne baliky na ILP, akademicka licencia zadarmo
SCIP [47] nekomercny program pre ILP
SYMPHONY v projekte COIN-OR [48]
Minisat [49] open source SAT solver, tiež Lingeling, glucose, CryptoMiniSat, painless
Concorde TSP solver [50] - riesi problem obchodneho cestujuceho so symetrickymi vzdialenostami, zadarmo na akademicke ucely
- Pre zaujimavost: TSP art [51]

ILP

Lineárny program:

Mame reálne premenné 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
Da sa riesit v polynomialnom case

Integer linear program

Program, v ktorom vsetky/vybrane premenne musia mat celociselne hodnoty, alebo dokonca povolime iba hodnoty 0 a 1.
NP uplny problem

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 najmensi pocet zo vstupnych mnozin tak, aby ich zjednotenie bola cela mnozina {1..m}
Binarne premenne x_i=1 ak vyberieme i-tu mnozinu
Chceme minimalizovat $\sum _{{i=1}}^{n}x_{i}\,$
za podmienky, ze pre kazde j z {1..m} plati $\sum _{{i: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 interagovali, 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 zvolit jadra a skorovacie tabulky, co je modelovaci, nie algoritmicky problem (mozeme skusit napr. nejake pravdepodobnostne 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}}-1$ 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?

Zdroj:

Jinbo Xu, Ming Li, Dongsup Kim, and Ying Xu. "RAPTOR: optimal protein threading by linear programming." Journal of bioinformatics and computational biology 1, no. 01 (2003): 95-117. [52]

CB12

Zhrnutie semestra a úvod do teórie grafov

vid prezentacie k cviceniu

Populacna genomika v UCSC genome browseri

Zopar zaujimavych polymorfizmov v ludskom genome

SNP rs1815739 CC: SNPedia, genome browser
SNP rs12255372 GT: SNPedia, genome browser
SNP rs2472297 TT: SNPedia, genome browser
Ďalšie zaujímavé SNPy: rs10427255 CC, rs671 GG, rs713598 GG, rs17822931 CT, rs4988235 CC, rs1042725 CC, rs7495174 AA, rs1426654 AA, rs4481887 AG
V browseri si vsimnite tracky (specificke pre verziu genomu hg19):
- HGDP Allele Freq s mapou sveta s distribuciou alel
- Genome Variants obsahuje genomy niekolkych ludi, napr Jima Watsona
- Takisto sa da pozriet genom ludi z jaskyne Denisova a Neandertalcov

UCSC genome browser ma aj dalsie tracky tykajuce sa populacnej genomiky a polymorfizmov

Pozrime si napriklad region chr2:46,570,000-46,630,000 v hg38
V casti Phenotype and Disease Associations
- napr. ClinVar Variants obsahuju asociacie variantov k chorobam
- GWAS Catalog sú výsledky GWAS štúdií

V starsej verzii ludskeho genomu hg18 je aj trojuholnikovy graf vazbovej nerovnovahy

region chr2:164,862-426,468 v hg18
zapnite "HapMap LD Phased" na Full (cast Variation and Repeats)
vsimnite si, ze miery LD sa medzi ludskymi podpopulaciami lisia (YRI: Nigeria; CEU: Europa; JPT+CHB: Japonsko, Cina)

Sekvenčné motívy, program MEME

Vazobne miesta transkripcnych faktorov sa casto reprezentuju ako sekvencne motivy
Ak mame skupinu sekvencii, mozeme hladat motiv, ktory maju spolocny
Znamy program na tento problem je MEME
Chodte na stranku http://meme-suite.org/
Zvolte nastroj MEME a v casti Input the primary sequences zvolte Type in sequences a zadajte tieto sekvencie
Pozrite si ostatne nastavenia. Co asi robia?
Ak server pocita dlho, mozete si pozriet vysledky tu

Gény, evolúcia a komparatívna genomika v UCSC genome browseri (cvičenie pri počítači)

Zobrazme si gén CLCA4 [53]
Zapnite si štandardnú sadu track-ov
Po kliknutí na gén si môžete prečítať o jeho funkcii, po kliknutí na ľavú lištu alebo na názov tracku v zozname na spodku stránky si môžete prečítať viac o tracku a meniť nastavenia
V tracku RefSeq genes si všimnite, že v tejto databáze má tento gén dve formy zostrihu, jedna z nich sa považuje za nekódujúcu, pretína sa aj s necharakterizovanou nekódujúcou RNA na opačnom vlákne
- Track RefSeq a jeho subtrack RefSeq Curated treba zapnut na pack
Nižšie vidíte track H3K27Ac Mark (Often Found Near Regulatory Elements) on 7 cell lines from ENCODE, kde bola táto histónová modifikácia v okolí génu detegovaná?
Všimnite si aj track DNase I Hypersensitivity, ktorý zobrazuje otvorený chromatin, prístupný pre viazanie transkripčných faktorov. Všimnite si jeho súvis s H3K27Ac trackom
Obidva tracky sú súčasťou tracku ENCODE regulation, v ktorom si môžete zapnúť aj ďalšie pod-tracky

Vsimnime si track Vertebrate Multiz Alignment & Conservation (100 Species)
- v spodnej casti tracku vidime zarovnania s roznymi inymi genomami
- v nastaveniach tracku zapnite Element Conservation (phastCons) na full a Conserved Elements na dense
- v tomto tracku vidíme PhyloP, co zobrazuje uroven konzerovanosti danej bazy len na zaklade jedneho stlpca zarovnania a dva vysledky z phyloHMM phastCons, ktory berie do uvahy aj okolite stlpce
Konkretne cast Conserved elements zobrazuje konkretne useky, ktore su najvac konzervovane
- Ak chceme zistit, kolko percent genomu tieto useky pokryvaju, ideme na modrej liste do casti Tools->Table browser, zvolime group Comparative genomics, track Conservation, table 100 Vert. El, region zvolime genome (v celom genome) a stlacime tlacidlo Summary/statistics, dostaneme nieco taketo:

item count	10,350,729
item bases	162,179,256 (5.32%)
item total	162,179,256 (5.32%)
smallest item	1
average item	16
biggest item	3,732
smallest score	186
average score	333
biggest score	1,000

- Ak by nas zaujimali iba velmi dlhe "conserved elements", v Table browser stlacime tlacidlo Filter a na dalsej obrazovke do policka Free-form query dame chromEnd-chromStart>=1500
- Potom mozeme skusit Summary/Statistics alebo vystup typu Hyperlinks to genome browser a Get output - dostaneme zoznam tychto elementov a kazdy si mozeme jednym klikom pozriet v browseri, napr. taketo
  - lod=24051 at chr1:50201403-50203312
  - lod=1899 at chr1:55663689-55667047 atd

Pozrime si teraz ten isty gen CLCA4 v starsej verzii genomu hg18 [54]
- 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 detegovany pozitivny vyber
- Po priblizeni do jedneho z exonov [55] vidite dosledky nesynonymnych mutacii

Poznamka: Existuju aj webservery na predikciu pozitivneho vyberu, napriklad tieto dva:

Selecton, clanok
Data monkey clanok
Skusili sme na Selecton poslat CLCA4 zo 7 cicavcov, subor tu: [56]
- vysledky [57] a [58] (metoda ale odporuca aspon 10 homologov)

@@ Riadok 10: / Riadok 10: @@
 * [https://www.youtube.com/playlist?list=PLU2XVjShDFwVeTDmo9Uv2NHWz3ijUjy5p playlist on youtube]
+'''Presentations'''
+* [[Media:prednasky-2022.pdf|All lectures in one pdf]]
+* [[Media:cb-2022.pdf|All tutorials for biology in one pdf]]
+* [[Media:ci-2022.pdf|All tutorials for computer science in one pdf]]
-'''Notes and presentations'''
+'''Notes'''
 L: lecture (everybody), TI: tutorial for computer science/informatics students,
@@ Riadok 28: / Riadok 33: @@
 | || colspan="3" | '' BV chapter 1'', [https://www.youtube.com/watch?v=3HI-wKmTuSs video 1] [https://www.youtube.com/watch?v=iHt5WPVFd7E video 2]
 |-
-| ''' TI:''' ||  Introduction to biology   ||    ||  [[CI01|notes]]
+| ''' TI:''' ||  Introduction to biology   ||    ||  [[#CI01|notes]]
 |-
 <!--za sablonu sa da pridat linka na video, ktory nie je parameter a tym padom nema problem s = -->
@@ Riadok 34: / Riadok 39: @@
 | || colspan="3" | '' ZB chapter 1'', [https://www.youtube.com/watch?v=EtPqUOHTJD4 video]
 |-
-| ''' TB:''' ||  Introduction to computer science, UCSC genome browser  ||    ||   [[CB01|notes]]
+| ''' TB:''' ||  Introduction to computer science, UCSC genome browser  ||    ||   [[#CB01|notes]]
 |-
 |-
@@ Riadok 47: / Riadok 52: @@
 | || colspan="3" | '' '' [https://www.youtube.com/watch?v=YvZlTL1qyUg video]
 |-
-| ''' TI:''' ||  Introduction to probability, genome coverage by sequencing reads  ||    ||  [[CI02|notes]]
+| ''' TI:''' ||  Introduction to probability, genome coverage by sequencing reads  ||    ||  [[#CI02|notes]]
 |-
 |-
-| ''' TB:''' ||  Introduction to dynamic programming, introduction to probability  ||    ||   [[CB02|notes]]
+| ''' TB:''' ||  Introduction to dynamic programming, introduction to probability  ||    ||   [[#CB02|notes]]
 |-
 |-
@@ Riadok 63: / Riadok 68: @@
 | || colspan="3" | '' BV chapter 2, DEKM chapter 2.1-2.4, 2.8, ZB chapter 4.1-4.4, 5.1-5.2'', [https://www.youtube.com/watch?v=0GkhkRiqbl4 video]
 |-
-| ''' TI:''' || Introduction to dynamic programming, proteomics  ||    ||  [[CI03|notes]]
+| ''' TI:''' || Introduction to dynamic programming, proteomics  ||    ||  [[#CI03|notes]]
 |-
 |-
-| ''' TB:''' || Dynamic programming for sequence alignment, dotplots  ||    ||  [[CB03|notes]]
+| ''' TB:''' || Dynamic programming for sequence alignment, dotplots  ||    ||  [[#CB03|notes]]
 |-
 |-
@@ Riadok 79: / Riadok 84: @@
 | || colspan="3" | '' BV chapter 2, DEKM chapter 2.5, 2.7, 6.1-6.3; ZB chapter 4.5-4.7, 5.3-5.5'', [https://www.youtube.com/watch?v=jkQGXiqKbxM video]
 |-
-| ''' TI:''' || Advanced algorithms for sequence alignment  ||   ||   [[CI04|notes]]
+| ''' TI:''' || Advanced algorithms for sequence alignment  ||   ||   [[#CI04|notes]]
 |-
 |-
-| ''' TB:''' || Programs for sequence alignment, scoring schemes ||   ||  [[CB04|notes]]
+| ''' TB:''' || Programs for sequence alignment, scoring schemes ||   ||  [[#CB04|notes]]
 |-
 |-
@@ Riadok 95: / Riadok 100: @@
 | || colspan="3" | '' BV chapter 4, DEKM chapter 3; ZB chapter 9.3, 10.4-10.7'', [https://www.youtube.com/watch?v=jkQGXiqKbxM video]
 |-
-| ''' TI:''' || Fast similarity search, BLAST, MinHash  ||    ||  [[CI05|notes]]
+| ''' TI:''' || Fast similarity search, BLAST, MinHash  ||    ||  [[#CI05|notes]]
 |-
 |-
-| ''' TB:''' || Hidden Markov models, E-value  ||      ||  [[CB05|notes]]
+| ''' TB:''' || Hidden Markov models, E-value  ||      ||  [[#CB05|notes]]
 |-
 |-
@@ Riadok 111: / Riadok 116: @@
 | || colspan="3" | '' BV chapter 3, DEKM chapter 7,8; ZB chapter 7, 8.1-8.2'', [https://youtu.be/RzfNHvZH4l8 video]
 |-
-| ''' TI:''' || Algorithms for HMM  ||     ||  [[CI06|notes]]
+| ''' TI:''' || Algorithms for HMM  ||     ||  [[#CI06|notes]]
 |-
 |-
-| ''' TB:''' || Substitution models, bootstrap, tree rooting ||    ||  [[CB06|notes]]
+| ''' TB:''' || Substitution models, bootstrap, tree rooting ||    ||  [[#CB06|notes]]
 |-
 |-
@@ Riadok 127: / Riadok 132: @@
 | || colspan="3" | '' BV chapter 5, ZB chapter 9.8, 10.8'', [https://youtu.be/1WM4QI2qx8A video]
 |-
-| ''' TI:''' || Substitution models  ||     ||  [[CI07|notes]]
+| ''' TI:''' || Substitution models  ||     ||  [[#CI07|notes]]
 |-
 |-
-| ''' TB:''' || Phylogenetic trees  ||    ||  [[CB07|notes]]
+| ''' TB:''' || Phylogenetic trees  ||    ||  [[#CB07|notes]]
 |-
 |-
@@ Riadok 143: / Riadok 148: @@
 | || colspan="3" | '' DEKM chapter 5.1, 11.5, ZB chapter 6.6,15.1,16.1-16.5,17.1'', [https://youtu.be/GFJ_oDV1KGU SK video]
 |-
-| ''' TI:''' || Felsenstein algorithm, algorithms for HMM and phyloHMM  ||    ||  [[CI08|notes]]
+| ''' TI:''' || Felsenstein algorithm, algorithms for HMM and phyloHMM  ||    ||  [[#CI08|notes]]
 |-
 |-
-| ''' TB:''' || K-means clustering, enrichment, multiple testing correction ||     ||  [[CB08|notes]]
+| ''' TB:''' || K-means clustering, enrichment, multiple testing correction ||     ||  [[#CB08|notes]]
 |-
 |-
@@ Riadok 163: / Riadok 168: @@
 | || colspan="3" | '' DEKM chapter 5; ZB chapter 4.8-4.10, 6.1-6.2, 13.1-13.2'', [https://youtu.be/ugMM81jZRpc SK video]
 |-
-| ''' TI:''' || Motif finding by EM and Gibbs sampling  ||    ||   [[CI09|notes]]
+| ''' TI:''' || Motif finding by EM and Gibbs sampling  ||    ||   [[#CI09|notes]]
 |-
 |-
-| ''' TB:''' || Introduction to context-free grammars  ||   ||  [[CB09|notes]]
+| ''' TB:''' || Introduction to context-free grammars  ||   ||  [[#CB09|notes]]
 |-
 |-
@@ Riadok 179: / Riadok 184: @@
 | || colspan="3" | '' DEKM chapter 10, ZB chapter 11.9'', [https://youtu.be/_Hh03Khsr9k SK video]
 |-
-| ''' TI:''' ||  Examples of biological databases, HW1 model solutions, introduction to context-free grammars  ||    ||  [[CI10|notes]]
+| ''' TI:''' ||  Examples of biological databases, HW1 model solutions, introduction to context-free grammars  ||    ||  [[#CI10|notes]]
 |-
 |-
-| ''' TB:''' ||  Genes, comparative genomics, Pfam   ||   ||  [[CB10|notes]]
+| ''' TB:''' ||  Genes, comparative genomics, Pfam   ||   ||  [[#CB10|notes]]
 |-
 |-
@@ Riadok 195: / Riadok 200: @@
 | || colspan="3" | '' '' [https://youtu.be/7vPt1vQX21M SK video]
 |-
-| ''' TI:''' ||  RNA structure   ||    ||  [[CI11|notes]]
+| ''' TI:''' ||  RNA structure   ||    ||  [[#CI11|notes]]
 |-
 |-
-| ''' TB:''' ||  Example of command-line tools   ||    ||  [[CB11|notes]]
+| ''' TB:''' ||  Example of command-line tools   ||    ||  [[#CB11|notes]]
 |-
 |-
@@ Riadok 208: / Riadok 213: @@
 |-
 |-
-| ''' TI:''' ||  Course summary, protein threading via integer linear programming  ||    ||   [[CI12|notes]]
+| ''' TI:''' ||  Course summary, protein threading via integer linear programming  ||    ||   [[#CI12|notes]]
 |-
 |-
-| ''' TB:''' ||  PSI-BLAST, microarray data, RNA structure, MEME, transcription factors in SGD, population genetics, course summary, graphs  ||     ||  [[CB12|notes]]
+| ''' TB:''' ||  PSI-BLAST, microarray data, RNA structure, MEME, transcription factors in SGD, population genetics, course summary, graphs  ||     ||  [[#CB12|notes]]
 |-
 |-
@@ Riadok 517: / Riadok 522: @@
 ==Použitie pravdepodobnosti na analýzu potrebného pokrytia pri sekvenovaní==
 <!-- Nerobili sme, uvedené pre zaujímavosť,-->
-Pozri [[CI02|cvičenia pre informatikov]]
+Pozri [[#CI02|cvičenia pre informatikov]]
 ==Dynamické programovanie==
@@ Riadok 1 043: / Riadok 1 048: @@
 * porovnanie lokalneho a globalneho zarovnania
 ** Dva proteiny z rôznych kvasiniek zarovnáme lokálne, globálne a globálne s tým, že neplatíme za medzery na koncoch
-* sekvencie a vysledne zarovnania: [[CB-aln-dp]]
+* sekvencie a vysledne zarovnania: [[#CB-aln-dp]]
 * vo vysledku si vsimnime, kolko ma kazde z nich %identity, %gaps, a kam sa zarovna sekvencia IRESPLGG ktora je na pozicii 29 v prvom a 30 v druhom proteine
@@ Riadok 1 421: / Riadok 1 426: @@
 ==Praktické príklady z predchádzajúcich cvičení==
-* [[CB04#Praktick.C3.A9_cvi.C4.8Denie_pri_po.C4.8D.C3.ADta.C4.8Di:_BLAT_vs_BLAST]] BLAT vs BLAST
+* [[#CB04#Praktick.C3.A9_cvi.C4.8Denie_pri_po.C4.8D.C3.ADta.C4.8Di:_BLAT_vs_BLAST]] BLAT vs BLAST
 =CI06=
 ==HMM opakovanie==
@@ Riadok 1 861: / Riadok 1 866: @@
 * na tvorbu stromov sa dá použiť IQ-TREE
 ** modely vid tu: https://github.com/Cibiv/IQ-TREE/wiki/Substitution-Models
-** vysledok [[CB:phylo]]
+** vysledok [[#CB:phylo]]
 * viewer napr. http://phylotree.hyphy.org/
@@ Riadok 1 923: / Riadok 1 928: @@
 ==Praktické príklady z predchádzajúcich cvičení==
-* [[CB01#QUAST:_program_na_.C5.A1tatistiky_o_kvalite_poskladania_gen.C3.B3mu]] Quast +  Galaxy
+* [[#CB01#QUAST:_program_na_.C5.A1tatistiky_o_kvalite_poskladania_gen.C3.B3mu]] Quast +  Galaxy
-* [[CB04#Pr.C3.ADklady_praktick.C3.BDch_programov]] lokálne vs globálne zarovnanie
+* [[#CB04#Pr.C3.ADklady_praktick.C3.BDch_programov]] lokálne vs globálne zarovnanie
 ==Hľadanie génov==

MBI 2022/2023: Rozdiel medzi revíziami

Aktuálna revízia z 20:38, 18. september 2023

Obsah

CI01

Úvod do biológie pre informatikov

CB01

Používanie počítačov v M 217

UCSC genome browser

Sekvenovanie v UCSC genome browseri

QUAST: program na štatistiky o kvalite poskladania genómu

Prehľad systému Galaxy

CI02

Úvod do pravdepodobnosti

Počítanie pokrytia genómov

Zhrnutie

CB02

Úvod do pravdepodobnosti

Pravdepodobnostný model náhodnej sekvencie

Použitie pravdepodobnosti na analýzu potrebného pokrytia pri sekvenovaní

Dynamické programovanie

CI03

Dynamické programovanie

Uvod do proteomiky

Sekvenovanie proteinov pomocou MS/MS

Problém 1

Problém 2

Problém 3

Detekcia znamych proteinov pomocou MS (nerobili sme)

CB03

Zarovnávanie sekvencií, opakovanie

Dotploty

Dynamické programovanie v Exceli

Práca so vzorcami v tabuľkovom procesore (Excel, LibreOffice, ...)

Späť k rozmieňaniu mincí

Zarovnávanie sekvencií v Exceli

CI04

Opakovanie dynamického programovania pre globálne zarovnanie

Reprezentácia pomocou grafu

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

Lokálne zarovnanie

Afínne skóre medzier

Nesprávne riešenie pomocou dynamického programovania

Správne riešenie pomocou dynamického programovanania

Lineárna pamäť: Hirshbergov algoritmus 1975

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

CB04

Skórovacie matice

Cvičenie pri počítači

Praktické cvičenie pri počítači: dotploty

Yass a dotploty

Príklady praktických programov

Plné dynamické programovanie

NCBI Blast

BLAT, chains, nets v UCSC browseri

Praktické cvičenie pri počítači: BLAT vs BLAST

BLAT/BLAST

RT PCR sekvencia z cvičenia vyššie

Uniprot pre projekty

CI05

Vzorec na výpočet senzitivity jadra

Minimizery: ako ušetriť pamäť a čas

MinHash

Odbočka do analýzy web-stránok: Podobnosť textov

Prvá idea: aproximácia vzorkovaním

Aproximácia Jaccardovej miery: MinHash

Hľadanie podobných sekvencií

Vypocet disperzie

CB05

Príklady stavových automatov pre HMM

E-hodnota (E-value) zarovnania

Praktické príklady z predchádzajúcich cvičení

CI06

HMM opakovanie

Inferencia pomocou najpravdepodobnejšej cesty, Viterbiho algoritmus

Inferencia - dopredný algoritmus

Inferencia - posterior decoding

Trénovanie HMM

Tvorba stavového priestoru modelu

Zovšeobecnené HMM

Párové HMM (pair HMM)