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

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


CI05: Rozdiel medzi revíziami

Z MBI
Prejsť na: navigácia, hľadanie
(Návrat do porovnávania sekvencií)
(31 intermediate revisions by the same user not shown)
Riadok 1: Riadok 1:
==Vzorec na vypocet senzitivity jadra==
+
==Vzorec na výpočet senzitivity jadra==
* Uvazujme jadro dlzky w (ako v programe BLAST pre nukleotidy)
+
* Uvažujme jadro dĺžky k (ako v programe BLAST pre nukleotidy, na prednáške sa dĺžka jadra označovala w, teraz k)
* Uvazujme pravdepodobnostny model zarovnania, v ktorom ma kazda pozicia pravdepodobnost p, ze bude zhoda a (1-p), ze bude nezhoda alebo medzera, zarovnanie ma dlzku L
+
* 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 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 <math>X_i=1, X_{i+1}=1, \dots, x_{i+w-1}=1</math>
+
* Nahodna premenna Y_i = 1 ak na pozicii i zacina jadro, t.j. ak <math>X_i=1, X_{i+1}=1, \dots, X_{i+k-1}=1</math>
* <math>P(Y_i = 1) = p^w</math>, nakolko X_i su navzajom nezavisle
+
* <math>P(Y_i = 1) = p^k</math>, nakolko X_i su navzajom nezavisle
* Nech <math>Y = \sum_{i=0}^{L-w} y_i</math>
+
* Nech <math>Y = \sum_{i=1}^{L-k-1} Y_i</math>
* Z linearity strednej hodnoty vieme lahko odhadnut E(Y) = (L-w+1)p^w
+
* Z linearity strednej hodnoty vieme lahko odhadnut <math>E(Y) = (L-k+1)p^k</math>
 
* Nas ale zaujima P(Y>0) = 1-P(Y=0)
 
* Nas ale zaujima P(Y>0) = 1-P(Y=0)
* <math>P(Y=0) = P(Y_0=0 \wedge \dots \wedge Y_{L-w}=0)</math>
+
* <math>P(Y=0) = P(Y_1=0 \wedge \dots \wedge Y_{L-k+1}=0)</math>
* Preco neplati, <math>P(Y=0) = P(Y_i = 0)^{L-w+1}</math>?
+
* Preco neplati, <math>P(Y=0) = P(Y_i = 0)^{L-k+1}</math>?
** Jendotlive Y_i nie su nezavisle, napr. <math>P(Y_{i+1}=1|Y_i=1)=p</math>
+
** Jednotlive Y_i nie su nezavisle, napr. <math>P(Y_{i+1}=1|Y_i=1)=p</math>
** V postupnosti Y_i sa jendotky maju tendenciu klastrovat spolu
+
** V postupnosti Y_i sa jednotky maju tendenciu zhlukovat spolu
* P(Y>0) ale vieme spoctat dynamickym programovanim
+
* Pravdepodobnost nepritomnosti jadra P(Y=0) ale vieme spocitat dynamickym programovanim
* Nech A[L] je pravdepodobnost vyskytu jadra v zarovnani dlzky L
+
* 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)
 
* <math>
 
* <math>
A[L] = \left\{\begin{array}{ll}
+
A[n] = \left\{\begin{array}{ll}
0 & \mbox{ak } L < w\\  
+
1 & \mbox{ak } n < k\\  
p^w+\sum_{i=0}^{w-1} p^i (1-p)A[L-i-1] & \mbox{ak } L \ge w\\  
+
\sum_{i=0}^{k-1} p^i (1-p)A[n-i-1] & \mbox{ak } n \ge k\\  
 
\end{array}\right.</math>
 
\end{array}\right.</math>
 +
V druhom riadku <math>p^i(1-p)</math> zodpoveda  <math>P(X_1...X_n\mbox{ konci presne }i\mbox{ jednotkami})</math> a A[n-i-1] je <math>P(X_1...X_{n-i-1}\mbox{ neobsahuje jadro})</math>, ale to je to iste ako <math>P(X_1...X_n\mbox{ neobsahuje jadro }| X_1...X_n\mbox{ konci presne }i\mbox{ jednotkami})</math>
  
==Minimizery: ako usetrit pamat==
+
==Minimizery: ako ušetriť pamäť a čas==
  
* k-merom nazveme k za sebou iducich pismen (nukleotidov DNA)
+
* 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
 
* 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á)
 +
<pre>
 +
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 
 +
</pre>
  
 
* 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
 
* 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 klastruju, mame sancu aspon jedno jadro z klastra najst
+
* 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
 +
<pre>
 +
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 
 +
</pre>
 +
 
 +
 
 +
 
 +
* 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
  
* 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
+
* 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
* Pri posune okna o jedno doprava casto najmensi k-mer zostava ten isty a netreba ho znovu ukladat, cim sa usetri pamat
+
* Priklad k=5, s=4
+
 
<pre>
 
<pre>
AGTGGCTGCCAGGCTGG    cGaGGCTGCCaGGtTGG  
+
AGTGGCTGCCAGGCTGG    cGaGGCTGCCtGGtTGG  
AGTGG*              CGAGG               
+
AGTGG,0              CGAGG               
 
  GTGGC                GAGGC             
 
  GTGGC                GAGGC             
 
   TGGCT                AGGCT*           
 
   TGGCT                AGGCT*           
   GGCTG               GGCTG           
+
   GGCTG,3              GGCTG           
     GCTGC*              GCTGC           
+
     GCTGC,4              GCTGC           
     CTGCC*              CTGCC*      
+
     CTGCC,5              CTGCC* <--   
 
       TGCCA                TGCCT         
 
       TGCCA                TGCCT         
 
       GCCAG                GCCTG       
 
       GCCAG                GCCTG       
         CCAGG*              CCTGG*     
+
         CCAGG,8              CCTGG*     
         CAGGC*              CTGGT*     
+
         CAGGC,9              CTGGT*     
           AGGCT*              TGGTT     
+
           AGGCT,10            TGGTT     
           GGCTG                GGTTG  
+
           GGCTG                GGTTG
 
             GCTGG                GTTGG   
 
             GCTGG                GTTGG   
 
</pre>
 
</pre>
  
* Rozdiel je pri hladani: v slovniku nehladame vsetky k-mery druhej sekvencie, ale tiez len minimizery, co moze usetrit cas
+
 
* 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 zdilajuce nejaky minimizer (binning) sa dostanu do jedneho zoznamu a budu uvazovane pri vypocte vzajomneho prekryvu
+
* 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
 
* 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 minimazermi boli casto sekvencie typu AAAAA, ktore su v biologickych sekvenciach nadreprezentovane a casto nie su funkcne zaujimave
+
* 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
 
* 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. [https://academic.oup.com/bioinformatics/article-abstract/32/14/2103/1742895]
 
* Li, Heng. "Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences." Bioinformatics 32.14 (2016): 2103-2110. [https://academic.oup.com/bioinformatics/article-abstract/32/14/2103/1742895]
Riadok 58: Riadok 109:
 
* 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. [https://academic.oup.com/bioinformatics/article/20/18/3363/202143]
 
* 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. [https://academic.oup.com/bioinformatics/article/20/18/3363/202143]
  
 +
* 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. [https://doi.org/10.1093/bioinformatics/btx235]
  
 
==MinHash==
 
==MinHash==
Riadok 74: Riadok 126:
 
Potom otázku "Ktoré dvojice textov sú podobné?" môžeme preformulovať napríklad ako "Ktoré dvojice textov majú Jaccardovu mieru podobnosti vyššiu ako <math>\alpha</math>?", kde <math>\alpha \in (0, 1)</math> je nejaká prahová hodnota.
 
Potom otázku "Ktoré dvojice textov sú podobné?" môžeme preformulovať napríklad ako "Ktoré dvojice textov majú Jaccardovu mieru podobnosti vyššiu ako <math>\alpha</math>?", kde <math>\alpha \in (0, 1)</math> je nejaká prahová hodnota.
  
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ť).
+
Ako rýchlo by sme vedeli spočítať Jaccardovu mieru pre dve množiny slov, každú s ''n'' prvkami?
  
=== Aproximácia Jaccardovej miery: MinHash ===
 
  
Nech je daná množina <math>A = \{a_1, a_2, \ldots, a_n\} \subseteq U</math>. Nech <math>h:U \to \mathbb{R}</math> je injektívna náhodná hashovacia funkcia (t.j. bez kolízií). Potom minimálny hash množiny <math>minHash_h(A)</math> je definovaný nasledovne:
+
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ť).
:::<math>minHash_h(A) := \min\{h(a_1), h(a_2), \ldots, h(a_n)\} =  \min_{1 \leq k \leq n} h(a_k)</math>
+
  
Keďže <math>h</math> je náhodná hashovacia funkcia, tak sa na hodnotu <math>minHash(A)</math> môžeme pozerať ako na náhodnú premennú, ktorá reprezentuje rovnomerne náhodný výber prvku z množiny <math>A</math>.
+
=== Prvá idea: aproximácia vzorkovaním ===
  
Nech <math>X</math> je náhodná premenná, ktorá nadobúda hodnotu 1, ak <math>minHash_h(A) = minHash_h(B)</math>, a inak hodnotu 0. Potom <math>Pr[X = 1] = J(A, B)</math>
+
* Ak by sme vedeli vzorkovať z <math>A \cup B</math> 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
 +
<math>X = \frac{1}{n}\sum_{i=1}^s X_i</math>
 +
kde X_i=1 ak u_i patri do prieniku a X_i=0 inak
  
Potom
+
* Toto sa podoba na zistovanie oblubenosti politika prieskumom verejnej mienky, u1, u2, ..., u_s  su "respondenti"
:::<math>E(X) = 0 \cdot Pr[X = 0] + 1 \cdot Pr[X = 1] = Pr[X = 1] = J(A, B)</math>.
+
  
Z toho vyplýva, že náhodná premenná <math>X</math> je nevychýleným odhadom Jaccardovej miery. Je to však veľmi nepohodlný odhad, lebo namiesto celej škály hodnôt od 0 po 1 vracia len dve možnosti.
+
Pre kazde X_i
 +
:::<math>E(X_i) = 0 \cdot Pr[X_i = 0] + 1 \cdot Pr[X_i = 1] = Pr[X_i = 1] = J(A, B)</math>.
 +
Z linearity strednej hodnoty
 +
:::<math>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)</math>.
 +
Z toho vyplýva, že náhodná premenná <math>X</math> je nevychýleným odhadom Jaccardovej miery.
  
V štatistike základnou mierou kvality nevychýleného odhadu slúži jeho variancia <math>Var(X) = E(X^2) - (E(X))^2</math>. Spočítajme si postupne obe hodnoty.
+
V štatistike základnou mierou kvality nevychýleného odhadu je jeho disperzia, odvodenie pozri nižšie
:::<math>E(X^2) = 0^2 \cdot Pr[X = 0] + 1^2 \cdot Pr[X = 1] = Pr[X = 1] = J(A, B)</math>
+
::::<math>Var(X) = \dfrac{J(A, B) - J^2(A, B)}{s} \le \dfrac{1}{4s}</math>
:::<math>(E(X))^2 = (J(A, B))^2</math>
+
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.
  
Čiže <math>Var(X) = J(A, B) - J^2(A, B)</math>. Aká je maximálna možná hodnota variancie?
+
'''Problémy:'''
 +
* nie je ľahké rovnomerne vzorkovať z <math>A \cup B</math>
 +
* 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 <math>E(X_i)=J(A,B)</math> a ktore sa budu lahsie pocitat
  
Táto otázka je ekvivalentná otázke "Aké je maximum funkcie <math>f(x) = x - x^2</math> na intervale <math>[0, 1]</math>?". Pre určenie extrémov hladkých funkcií treba nájsť korene jej prvej derivácie. Derivácia tejto funkcie je <math>f'(x) = 1 - 2x</math>, jej koreň je hodnota <math>0.5</math>. Hodnota funkcie v tomto bode je rovná <math>0.25</math>. Čiže <math>Var(X) \leq 0.25</math>.
+
=== Aproximácia Jaccardovej miery: MinHash ===
  
Ako môžeme tento odhad zlepšiť?
+
Budeme mať náhodné hašovacie funkcie <math>h_1, h_2, \dots, h_s</math>.
  
Jedna z možností je zobrať viacero nezávislých hashovacích funkcií <math>h_1, h_2, \ldots, h_k</math>, a spočítať <math>minHash_{h_1}, \ldots, minHash_{h_k}</math> pre obidve množiny. Označme si príslušné náhodné premenné ako <math>X_1, X_2, \ldots, X_k</math>. Každá z nich má strednú hodnotu <math>E(X_i) = E(X) = J(A, B)</math> a rovnakú varianciu <math>Var(X_i) = Var(X) = J(A,B) - J^2(A, B)</math>.  
+
O každej hašovacej funkcii predpokladáme, že ak ju použijeme na nejakú množinu <math>A = \{a_1, a_2, \ldots, a_n\}</math>, tak <math>h(a_1), h(a_2), \ldots, h(a_n)</math> bude náhodná permutácia množiny <math>A</math>.
  
Nech náhodná premenná <math>Y_k := \dfrac{X_1+X_2+\ldots+X_k}{k}</math> je ich priemer. Potom jej stredná hodnota je rovná <math>E(Y_k) = E\left(\frac{X_1+X_2+\ldots+X_k}{k}\right) = k^{-1} E(X_1+X_2+\ldots+X_k) = k^{-1} [ E(X_1)+E(X_2)+\ldots+E(X_k)] = k^{-1} k E(X) = E(X) = J(A, B)</math>. Čiže aj <math>Y_k</math> je nevychýleným odhadom Jaccardovej miery.
+
Pre množinu <math>A = \{a_1, a_2, \ldots, a_n\}</math> a hašovaciu funkciu ''h'' je <math>minHash_{h}(A)</math> je definovaný nasledovne:
 +
:::<math>minHash_h(A) := \min\{h(a_1), h(a_2), \ldots, h(a_n)\}</math>
  
Pozrieme sa na jej varianciu: <math>Var(Y_k) = Var\left(\frac{X_1+X_2+\ldots+X_k}{k}\right) = \dfrac{1}{k^2} Var(X_1+X_2+\ldots+X_k) \overset{*}{=} \dfrac{1}{k^2} [Var(X_1) + \ldots Var(X_k)] = \dfrac{1}{k^2} k \cdot Var(X) = \dfrac{Var(X)}{k} \leq \dfrac{1}{4k}</math>  
+
Keďže <math>h</math> je náhodná hašovacia funkcia, tak sa na hodnotu <math>minHash(A)</math> môžeme pozerať ako na náhodnú premennú, ktorá reprezentuje rovnomerne náhodný výber prvku z množiny <math>A</math>.
  
''(*) tento prechod je možný len kvôli tomu, že premenné <math>X_i</math> sú nezávislé.''
+
Nech <math>X_i</math> je náhodná premenná, ktorá nadobúda hodnotu 1, ak <math>minHash_{h_i}(A) = minHash_{h_i}(B)</math>, a inak hodnotu 0. Potom <math>E[X_i] = J(A, B)</math>, lebo celkovo máme <math>|A\cup B|</math> prvkov a <math>J(A,B) = |A\cap B|/|A\cup B|</math> značí, aké percento z nich je v prieniku.
  
Vidíme teda, že varianciu (resp. kvalitu) môžeme potlačiť na ľubovoľne malú postupným zvýšením počtu hashov.
+
Takýmito hodnotami <math>X_i</math> teda nahradíme náhodné vzorky diskutované vyššie.
  
Všimnite si, že premenná <math>k Y_k</math> (t.j. nie priemer, ale súčet jednotlivých <math>X_i</math>) je súčtom nezávislých indikátorov s rovnakou distribúciou, a teda má binomické rozdelenie s parametrami <math>n=k</math> a <math>p=J(A, B)</math>.
+
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 <math>J(S_A, S_B)</math> kde <math>s_A</math> je mnozina hodnot v sketchi mnoziny <math>A</math>. To usetri cas pri vypocte sketchu.
  
Druhá možnosť zlepšenia je nechať jednu hashovaciu funkciu, ale porovnávať nie 1, ale k najmenších hashov dvoch množín. Vedie to ku podobnému asymptotickému správaniu.
 
  
 
* 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. [https://www.cs.princeton.edu/courses/archive/spring13/cos598C/broder97resemblance.pdf]
 
* 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. [https://www.cs.princeton.edu/courses/archive/spring13/cos598C/broder97resemblance.pdf]
  
 +
=== Hľadanie podobných sekvencií ===
  
=== Návrat do porovnávania 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)
  
Ako "slová" použijeme všetky súvislé podreťazce fixnej dĺžky (dĺžka sa tradične označuje ako <math>k</math>) danej sekvencie. Tieto podreťazce sa tradične označujú ako <math>k</math>-mery.
 
Potom na hľadanie dvoch podobných čítaní z množiny čítaní môžeme použiť minhash.
 
  
 
* 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. [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x]
 
* 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. [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x]
 +
 +
===Vypocet disperzie===
 +
 +
Binarna premenna X_i, kde P(X_i=1)=p.
 +
Disperzia <math>Var(X_i) = E(X_i^2) - (E(X_i))^2</math>. Spočítajme si postupne obe hodnoty.
 +
:::<math>E(X_i^2) = 0^2 \cdot Pr[X_i = 0] + 1^2 \cdot Pr[X_i = 1] = Pr[X_i = 1] = p</math>
 +
:::<math>(E(X_i))^2 = p^2</math>
 +
 +
Čiže <math>Var(X_i) = p - p^2</math>. Aká je maximálna možná hodnota disperzie?
 +
 +
Táto otázka je ekvivalentná otázke "Aké je maximum funkcie <math>f(x) = x - x^2</math> na intervale <math>[0, 1]</math>?". Pre určenie extrémov hladkých funkcií treba nájsť korene jej prvej derivácie. Derivácia tejto funkcie je <math>f'(x) = 1 - 2x</math>, jej koreň je hodnota <math>0.5</math>. Hodnota funkcie v tomto bode je rovná <math>0.25</math>. Čiže <math>Var(X_i) \leq 0.25</math>.
 +
 +
<math>X =\frac{1}{s}\sum_{i=1}^s X_i</math> kde X_i su nezavisle premenne a každá z nich má strednú hodnotu <math>E(X_i) = E(X) = p</math> a rovnakú disperziu <math>Var(X_i) = Var(X) = p - p^2</math>.
 +
T.j. X je ich priemer.
 +
 +
Pozrieme sa na jej varianciu: <math>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}</math>
 +
 +
''(*) tento prechod je možný len kvôli tomu, že premenné <math>X_i</math> 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á <math>s X</math> (t.j. nie priemer, ale súčet jednotlivých <math>X_i</math>) je súčtom nezávislých indikátorov s rovnakou distribúciou, a teda má binomické rozdelenie s parametrami <math>s</math> a <math>p</math>.

Verzia zo dňa a času 15:55, 20. október 2022

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, P(Y=0)=P(Y_{i}=0)^{{L-k+1}}?
    • 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. [1]
  • 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. [2]
  • 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. [3]

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. [4]

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. [5]

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.