Bioinformatikai Algoritmusok 9 GY Overlap Layout Consensus Assembly

Bioinformatikai Algoritmusok 9. GY Overlap Layout Consensus Assembly Nyíri Tamás nytuaai@gmail. com http: //people. inf. elte. hu/nytuaai (Ben Langmead diasora alapján)

Assembly módszerek a gyakorlatban OLC: Overlap-Layout-Consensus assembly DBG: De Bruijn gráf alapú assembly Mindkettő kihagyja a kezelhetetlen ismétléseket. A kezelhetetlen ismétlések darabokra törik az assembly-t. Ezeket a darabokat contig-oknak hívjuk (a contiguous szó rövidítése). a_long_long_time Assembly mohó SCS használatával a_long_time Assembly OLC vagy DBG használatával a_long_time

Assembly alternatívák 1. alternatíva: Overlap-Layout-Consensus (OLC) assembly 2. alternatíva: de Bruijn gráf alapú (DBG) assembly Overlap Error correction Layout de Bruijn graph Consensus Refine Scaﬀolding

Overlap Layout Consensus Overlap Átfedési gráf felépítése Layout Az átfedési gráf összefüggő részeinek contig-okba gyűjtése Consensus Mindegyik contig-ra a legvalószínűbb nukleotid kiválasztása

Átfedések keresése Lehetnénk ennél naivabbak? Mondjuk l = 3 Terjesszük ki balra; itt láthatjuk, hogy Y 6 hosszú prefixe egyezik X szuffixével Keressük ezt Y-ban, jobbról balra haladva X: CTCTAGGCC Y: TAGGCCCTC Megtaláltuk Ezt végigcsináljuk minden input string párra

Átfedések keresése Tudunk szuffix fákat használni az átfedésekre? Feladat: Adott string-ek egy S halmaza, mindegyik x eleme S string-re keressük meg az összes olyan átfedést amelyben szerepel x egy prefix-e és egy másik y string egy szuffix-e Segítség: Építsünk egy általános szuffix fát S string-jeire

Átfedések keresése szuffix fával Általános szuffix fa { “GACATA”, “ATAGAC” }-ra A $0 C 6 $0 C TA 5 ATA$ 0 $ 1 1 11 $0 3 GAC$ 1 7 GAC TA 13 GAC$ 1 9 $1 GACATA$0 ATAGAC$1 ATA$ 0 $ 1 2 12 ATA$ 0 $ 1 0 10 $0 4 GAC$ 1 8 Legyen a lekérés = GACATA. A gyökértől indulva kövessük a lekéréssel címkézett utat. A zöld él azt mutatja, hogy a 2. string 3 hosszú szuffix-e egyenlő a lekérés 3 hosszú prefixével. ATAGAC ||| GACATA

Átfedések keresése szuffix fával Általános szuffix fa { “GACATA”, “ATAGAC”}-ra A $0 C 6 $0 C TA 5 ATA$ 0 $ 1 1 11 $0 3 GAC$ 1 7 GAC TA 13 GAC$ 1 9 $1 GACATA$0 ATAGAC$1 ATA$ 0 $ 1 2 12 ATA$ 0 $ 1 0 10 $0 4 GAC$ 1 8 Stratégia: (1) Építsük fel a fát (2) Mindegyik string-re: Sétáljunk le a gyökértől és jegyezzünk fel minden szeparátorral elválasztott kimenő élt. Minden ilyen él megfelel egy olyan prefix/szuffix egyezésnek amelyben a lekért string prefixe egyezik egy szeparátorral végződő string szuffixével.

Átfedések keresése szuffix fával Általános szuffix fa { “GACATA”, “ATAGAC” }-re GACATA | ATAGAC $0 A 6 C TA 5 1 GAC$ 1 9 ATA$ 0 $ 1 11 $0 C $0 3 GAC$ 1 7 GACATA ||| ATAGAC GACATA$0 ATAGAC$1 $1 GAC TA 13 ATA$ 0 $ 1 2 12 ATA$ 0 $ 1 0 10 $0 4 ATAGAC ||| GACATA Most legyen a lekérés: ATAGAC GAC$ 1 8

Átfedések keresése szuffix fával A $0 C 6 $0 C TA 5 GAC$ 1 9 ATA$ 0 $ 1 1 11 $0 3 $1 GAC TA 13 ATA$ 0 $ 1 2 12 ATA$ 0 $ 1 0 10 $0 4 GAC$ 1 8 GAC$ 1 7 Legyen d darab n hosszú read, teljes hossz: N = dn, és a = egymást fedő read párok száma Adott string párra csak a leghosszabb szuffix/prefix találatot jegyezzük fel Általános szuffix fa építése: O(N) . . . A piros utak végigjárása: . . . Átfedések megtalálása (zöld): O(N) Teljes: O(N + a) O(a) d 2 nem jelenik meg explicit, de a = O(d 2) legrosszabb esetben

Átfedések keresése Mi van ha meg akarjuk engedni az eltéréseket és a réseket az átfedésekben? Így hogy tudjuk megtalálni X egy szuffixének és Y egy prefixének legjobb átfedését? X: CTCGGCCCTAGG ||||| Y: GGCTCTAGGCCC Dinamikus programozás De úgy kell megfogalmaznunk a feladatot, hogy csak az X szuffixét és Y prefixét tartalmazó backtrace-ek engedélyezettek

Átfedések keresése dinamikus programozással Keressük meg X egy szuffixének Y egy prefixére való legjobb illeszkedését X: CTCGGCCCTAGG ||||| Y: GGCTCTAGGCCC global alignment recurrence-t és score függvényt használunk s(a, b) D[i, j] = min D[i — 1, j] + s(x[i — 1], —) D[i, j — 1] + s(—, y[j — 1]) D[i — 1, j — 1] + s(x[i — 1], y[j — 1]) De hogyan tudjuk rávenni, hogy megtalálja a prefix / szuﬃx illeszkedéseket? A A 0 C 4 G 2 T 4 ‐‐‐ 8 C 4 0 4 2 8 G 2 4 0 4 8 T 4 2 4 0 8 ‐‐‐ 8 8

Átfedések keresése dinamikus programozással Keressük meg X egy szuffixének Y egy prefixére való legjobb illeszkedését D[i, j] = min D[i — 1, j] + s(x[i — 1], —) D[i, j — 1] + s(—, y[j — 1]) D[i — 1, j — 1] + s(x[i — 1], y[j — 1]) s(a, b) A A C G T ‐‐‐ 0 4 2 4 8 C 4 0 4 2 8 G 2 4 0 4 8 T 4 2 4 0 8 ‐‐‐ 8 8 Y ‐‐‐ G G C T A G G C C C Hogyan inicializáljuk az első sort és oszlopot, hogy X szuﬃxe illeszkedjen Y prefixére? Az első oszlop 0 -kat kap (X bármilyen szuffixe lehetséges) Az első sor ∞-ket kap (prefixe kell hogy legyen Y-nak) Haladjunk visszafelé az utolsó sorból ‐‐‐ 0 ∞ ∞ ∞ X C T C G G C C C T A G G 0 0 0 4 12 20 2 8 36 44 52 60 68 76 84 9 2 CTCGGCCCTAGG 0 28 36 44 52 60 68 76 8 4 4 8 14 2 X: 4 8 8 6 20 28||| ||||| 36 44 52 60 68 7 6 2 20 24 GGCTCTAGGCCC 30 36 44 52 60 68 0 4 12 1 Y: 0 0 8 16 16 24 26 30 36 44 52 60 4 4 0 8 16 18 26 30 34 36 44 52 4 8 4 2 8 16 22 30 34 34 36 44 4 8 8 6 2 10 18 26 34 34 34 36 4 8 10 8 8 2 10 18 26 34 36 36 2 6 12 14 12 10 18 26 34 40 0 2 10 16 18 16 10 0 10 18 26 34 0 0 6 14 20 22 18 10 2 10 18 26

Átfedések keresése dinamikus programozással Keressük meg X egy szuffixének Y egy prefixére való legjobb illeszkedését D[i, j] = min D[i — 1, j] + s(x[i — 1], —) D[i, j — 1] + s(—, y[j — 1]) D[i — 1, j — 1] + s(x[i — 1], y[j — 1]) s(a, b) A A C G T ‐‐‐ 0 4 2 4 8 C 4 0 4 2 8 G 2 4 0 4 8 T 4 2 4 0 8 ‐‐‐ 8 8 Y ‐‐‐ Probléma: a nagyon rövid találatok C véletlenül nagyon magas T pontszámokat kaphatnak. . . C G. . . ami elfedheti a relevánsabb G találatokat X C Mondjuk, hogy a minimum átfedési hossz l = 5 C C T A G G ‐‐‐ G G C T A G G C C C 0 0 0 0 ∞ ∞∞∞∞∞ 4 12 20 28 36 44 52 60 68 76 84 92 4 8 14 20 28 36 44 52 60 68 76 84 4 8 8 16 20 28 36 44 52 60 68 76 0 4 12 12 20 24 30 36 44 52 60 68 0 0 8 16 16 24 26 30 36 44 52 60 4 4 0 8 16 18 26 30 34 36 44 52 4 8 4 2 8 16 22 30 34 34 36 44 4 8 8 6 2 10 18 26 34 34 34 36 4 8 10 8 8 2 10 18 26 34 36 36 2 6 12 14 12 10 18 26 34 40 0 2 10 16 18 16 10 0 10 18 26 34 0 0 6 14 20 22 18 10 2 10 18 26

Átfedések keresése dinamikus programozással Keressük meg X egy szuffixének Y egy prefixére való legjobb illeszkedését D[i, j] = min s(a, b) A D[i — 1, j] + s(x[i — 1], —) D[i, j — 1] + s(—, y[j — 1]) D[i — 1, j — 1] + s(x[i — 1], y[j — 1]) A C G T ‐‐‐ C 4 0 4 2 8 0 4 2 4 8 G 2 4 0 4 8 T 4 2 4 0 8 ‐‐‐ 8 8 Y Oldjuk meg úgy, hogy néhány más cellát is ∞-re inicializálunk Azokat a cellákat amiknek így megváltozott az értéke, pirossal jelöljük Így a releváns találatok a legjobbak ‐‐‐ G G C T A G G C C C 0 0 0 0 0 ∞∞∞∞∞∞ 4 4 4 0 0 4 4 12 20 12 8 8 4 0 4 8 8 8 6 12 20 14 8 12 8 0 4 8 10 12 10 28 20 16 12 16 8 2 6 8 14 16 ∞ ∞∞∞∞∞ 36 28 20 20 16 16 8 2 8 12 18 20 44 36 28 24 24 18 16 10 2 10 16 22 52 44 36 30 26 26 22 18 10 2 10 18 60 52 44 36 30 30 30 26 18 10 0 10 68 60 52 44 36 34 34 34 26 18 10 2 76 68 60 52 44 36 34 34 34 26 18 10 84 76 68 60 52 44 36 34 26 18 92 84 76 68 60 52 44 36 36 40 34 26

Átfedések keresése dinamikus programozással Legyen d n hosszú read, a teljes hossz N = dn, és a az átfedéssel rendelkező párok száma Kipróbálandó átfedések száma: O(d 2) Egy dinamikus programozási mátrix mérete: O(n 2) Teljes költség: O(d 2 n 2) = O(N 2) Hasonlítsuk össze O(N 2)-t a szuffix fával: O(N + a), de ahol a a legrosszabb esetben O(d 2) De a dinamikus programozás flexibilisebb, megenged hibákat és réseket A valódi átfedő programok keverik a kettőt, és indexeket használnak, hogy kiszűrjék a nagyrészét a nem-átfedő pároknak, aztán dinamikus programozást használnak a maradék párokra.

Átfedések keresése Az átfedések keresése általában az assembly leglassabb része Képzeljünk el egy második generációs szekvenciálásból származó adathalmazt 10^8 -10^9 read-el! Az illesztési technikák is hasznosak lehetnek az átfedések keresésére Láttunk példákat a naive exact matching, szuﬃx-fa-segített exact matching, és a dinamikus programozás módosításaira Felhasználhattuk volna az eﬃcient exact matching-et, az approximate string matching-et, a co-traversal-t, stb…

Átfedések keresése A Celera Assembler átfedő program talán a legjobbak dokumentált: Invertált substring indexek read batch-ekre építve Csak olyan read-ek közti átfedéseket keressünk amelyek egy vagy több tetszőlegesen hosszú substring-ben megegyeznek http: //wgs-assembler. sourceforge. net/wiki/index. php/Run. CA#Overlapper

Overlap Layout Consensus Overlap Átfedési gráf felépítése Layout Az átfedési gráf összefüggő részeinek contig-okba gyűjtése Consensus Mindegyik contig-ra a legvalószínűbb nukleotid kiválasztása

4 6 5 ry_thin 6 y_thing 5 4 6 5 thing_t 4 6 hing_tu 5 6 5 4 ing_tur 4 6 ng_turn 4 5 6 5 g_turn_ 6 _turn_t 4 45 urn_tur 5 6 6 rn_turn 6 54 6 4 4 4 n_turn_ 5 54 turn_tu 4 6 5 urn_the 4 6 rn_ther 5 6 5 4 n_there 4 6 5 _there_ 6 4 there_i 5 6 4 5 here_is 4 6 ere_is_ 5 6 4 re_is_a 4 6 5 4 5 e_is_a_ 6 4 6 5 turn_th 4 5 Az átfedési gráfok nagyok és bonyolultak. Nem látni egyértelműen a contig-okat. 5 Lent: a következő string-hez tartozó overlap gráf egy része: _thing_ to_every_thing_turn_turn_there_is_a_season 4 l = 4, k = 7 6 Layout 4

Layout to_ever 6 Látunk bármi redundánsat az átfedési gráf ezen részében? Néhány élet következtetni lehet (tranzitivitás miatt) más élek alapján 5 o_every 6 4 _every_ 5 6 4 5 Pl. a zöld él kikövetkeztethető a kék-ből every_t 4 6 very_th 5 6 4 ery_thi 4 5 6 5 4 ry_thin 6 y_thing 5 6 _thing_ 4 5

4 6 5 ry_thin 6 y_thing 5 _thing_ 5 4 6 5 thing_t 4 6 hing_tu 5 6 5 4 ing_tur 4 6 ng_turn 4 5 6 5 g_turn_ 6 _turn_t 4 45 urn_tur 5 6 6 rn_turn 6 54 6 4 4 4 n_turn_ 5 54 turn_tu 4 6 5 urn_the 4 6 rn_ther 5 6 5 4 n_there 4 6 5 _there_ 6 4 there_i 5 6 4 5 here_is 4 6 ere_is_ 5 6 4 re_is_a 4 6 5 4 5 e_is_a_ 6 4 6 5 turn_th 4 5 Távolítsuk el a tranzitivitás miatt kikövetkeztethető éleket, kezdve az 1 csúcsot kihagyó élekkel: x 4 Előtte: 6 Layout 4

y_thing 6 thing_t 6 4 hing_tu 6 ing_tur 4 6 ng_turn 6 g_turn_ 4 urn_tur 6 6 rn_turn 4 6 n_turn_ 4 6 _turn_t 4 6 6 turn_th 4 turn_tu 6 urn_the 4 6 rn_ther 6 n_there 6 _there_ 4 6 there_i 6 here_is 6 4 ere_is_ 6 Távolítsuk el a tranzitivitás miatt kikövetkeztethető éleket, kezdve az 1 csúcsot kihagyó élekkel: x Utána: 6 Layout _thing_

6 ry_thin 6 _thing_ 6 thing_t 6 hing_tu 6 ing_tur 6 ng_turn 6 g_turn_ x 4 urn_tur 6 4 rn_turn 6 6 n_turn_ 6 4 _turn_t 6 turn_tu Távolítsuk el a tranzitivitás miatt kikövetkeztethető éleket, kezdve az 1 csúcsot kihagyó élekkel: x Utána: Még egyszerűbb 6 6 turn_th 6 urn_the 6 rn_ther 6 n_there 6 _there_ 6 there_i 6 here_is 6 ere_is_ Layout y_thing

to_every_ 6 every_t 6 Contig 1 6 ery_thi 6 ry_thin 6 y_thing 6 _thing_ 6 thing_t 6 hing_tu 6 ing_tur 6 ng_turn 6 g_turn_ 4 urn_tur 6 4 rn_turn 6 6 n_turn_ 6 4 Contig 2 turn_there_is_a_season Feloldhatatlan ismétlés to_every_thing_turn_ very_th _turn_t 6 turn_tu 6 turn_th 6 urn_the 6 rn_ther 6 n_there 6 _there_ 6 there_i 6 here_is 6 ere_is_ 6 re_is_a 6 e_is_a_ 6 _is_a_s 6 is_a_se 6 s_a_sea 6 _a_seas 6 a_seaso 6 _season Vegyük az egyértelműen meghatározható contig-okat 6 Layout 6 o_every

Layout A gyakorlatban a layout lépésnek a hamis szubgráfokkal is törődnie kell, pl. szekvenciálási hiba miatt Lehetséges ismétlési határ a vágás b Eltérés b a . . . Az eltérések szekvenciálási hiba vagy ismétlés miatt is keletkezhetnek. Mivel a b-n keresztül menő út hirtelen megszakad, ezt lehet hogy hibának vesszük és vágunk.

Overlap Layout Consensus Átfedési gráf felépítése Az átfedési gráf összefüggő részeinek contig-okba gyűjtése Mindegyik contig-ra a legvalószínűbb nukleotid kiválasztása

Consensus TAGATTACACAGATTACTGA TTGATGGCGTAA CTA TAGATTACACAGATTACTGACTTGATGGCGTAAACTA TAG TTACACAGATTATTGACTTCATGGCGTAA CTA TAGATTACACAGATTACTGACTTGATGGCGTAA CTA Soroljuk fel a contigot alkotó read-eket Fogjuk a consensus-t (pl. többségi szavazat) Mindegyik pozíciónál kérdezzük a következőt: milyen nukleotid (és/vagy rés) van itt? Komplikációk: (a) szekvenciálási hiba, (b) ploiditás Mondjuk hogy a valódi genotipus AG, de sok szekvenciálási hiba van csak kb. 6 read fedi a pozíciót.

Overlap Layout Consensus Overlap Átfedési gráf felépítése Layout Az átfedési gráf összefüggő részeinek contig-okba gyűjtése Consensus Mindegyik contig-ra a legvalószínűbb nukleotid kiválasztása OLC hátrányai Az átfedési gráf építése lassú. Láttuk az O(N + a) és O(N 2) –es algoritmusokat. Az átfedési gráf nagy; egy csúcs read-enként és gyakorlatban az élek száma szuperlineárisan nő a read-ek számának függvényében A 2. generációs szekvenciálásból származó adathalmazok nagyságrendileg 10^8 -10^9 read-et és 10^11 nukleotidot tartalmaznak

Assembly alternatívák 1. alternatíva: Overlap-Layout-Consensus (OLC) assembly 2. alternatíva: de Bruijn gráf alapú (DBG) assembly Overlap Error correction Layout de Bruijn graph Consensus Refine Scaﬀolding

Házi Feladat (5 pont + ? bónusz) 1. Generáljunk egy 100 karakter hosszú random (25%-25%-25%) DNS szekvenciát. (1 pont) 2. Mintavételezzünk a darab random b hosszú read-eket ebből a generált stringből. Figyeljünk oda, hogy ha kétszer egymás után ne vételezzünk ugyanolyan string-et. A paramétereket úgy válasszuk meg hogy lehetőleg ne fusson órákig a program (emlékezzünk: SCS NP-teljes), de mégis közel optimális eredményt adjon. (1 pont) 3. Futtassuk le a Shortest Common Superstring és a Greedy Shortest Common Superstring algoritmusokat a kreált read-ekre. Jegyezzük le az eredményeket illetve mérjük (https: //docs. python. org/2/library/time. html) és jegyezzük le a futási időt is a két algoritmusnál. A Greedy SCS paraméterét megint csak próbáljuk megtippelni, kikísérleztezni. (1 pont) 4. Mérjük meg az eltérést a SCS-el generált eredmény és az eredeti DNS string között Edit Distance használatával. Mérjük meg az eltérést a Greedy SCS-el generált eredmény és az eredeti string között a már régi gyakorlatokról ismert Edit Distance használatával. (1 pont) 5. Egy loop-ban ismételjük meg a folyamatot rögzített DNS-re és read-ekre más és más paraméterekkel, illetve különböző DNS-ekre és read-ekre. A kapott futási időket, Edit Distance-okat, stb. a különböző futtatásoknál grafikonokkal szemléltessük (https: //matplotlib. org) és próbáljuk megmagyarázni a kapott eredményeket írásban, hivatkozva az előállított grafikonokra. Törekedjünk a minél tömörebb de mégis minél informatívabb ábrákra, magyarázatokra, paraméter választásokra. (1 pont a minimális erőfeszítésért + bónusz pontok extra erőfeszítésért) • • Az eredményeket Ipython notebook (. ipynb) formátumban várom e-mail címemre (nytuaai@gmail. com) Kérdéseket, megjegyzéseket a házi feladattal vagy gyakorlatokkal kapcsolatban is ide lehet küldeni.

Update: Követelmények A félév végi jegy a házi feladatokra kapott pontszámokból fog összeállni. Ez két részre lesz bontva: 1. minden az Assembly-ig (előző 6 házi – 2 pont / házi) 2. minden az Assembly-től (maradék 2 házi ehetivel együtt – 5 pont / házi) • • Mindkét részből minimum 50%-ot el kell érni a 2 -esért! • • • Egyébként a pontozás a következőképpen alakul: 11 -13 pont: 2 14 -16 pont: 3 17 -19 pont: 4 20 -22 pont: 5 • A házi feladatok határideje a szorgalmi időszak utolsó hete! • Elegendő e-mail-ben elküldeni őket, nem kell személyesen bemutatni.