A genom összeszerelések egyik legfontosabb hozzávalója a minél hosszabb szekvenciák felhasználása. Ebben a cikkben összeszedték a különböző összeszerelési stratégiákat, és az első ábrán is jól látszik, hogy a hosszú readek felhasználásával az N50 nem fog számottevően változni a rövid read-es össszeszereléshez képest, de a hibás összeszerelések száma jóval alacsonyabb lesz.
Hosszú szekvenciákat Nanopore-ből vagy PacBio-ból szerezhetünk. Habár a Nanopore jó megoldás baktériumok, vírusok esetén, nagyobb genomokat csak korlátozottan lehet összerakni vele. Tehát az egyedüli megoldásnak a PacBio tekinthető.
Nézzük, hogyan is néz ki egy genom összeszerelés ezzel a technológiával. Előrebocsátom, hogy az itt bemutatott programok aktív fejlesztés alatt állnak, ezért nem biztos, hogy az itt leírtak egy év múlva is érvényesek lesznek. A dokumentáció pedig szokás szerint előbb avul el, mint a kód.
A PacBio denovo programját Falconnak hívják és a teljes kód fent van a GitHubon, csak ember legyen a talpán, aki kiigazodik a sok félbehagyott repo, átnevezett programok és verziók között. Elég csak megnézni ezt a listát. Az minden esetre látszik, hogy a program három nagyobb részegységből áll: nyers denovo (ez lenne a FALCON, vagy falcon3. A kettes verzió valahol eltűnt egy fekete lyukban.), haplotípus meghatározás (FALCON-Unzip, falcon_unzip3), végső simítások (FALCON-polish).
Az egészet összefogja a pb-assembly, ami egy conda csomag, hogy legalább a telepítéssel ne legyen gond. Igazából máshogy nem is lehet telepíteni. Megpróbálkoztam a fordítással is, csak, hogy lássam, képes vagyok rá, de nem voltam. Nem érdemes időt vesztegetni rá.
A programok vezérlését egyetlen konfigurációs fájllal lehet elvégezni. Felépítésük a régi Windows-os INI állományokra hajaz. Szögletes zárójelbe megnevezünk egy szekciót, majd kulcs-érték párokat definiálunk. A program saját munkafolyamat kezelővel rendelkezik, ami azt jelenti, hogy hiba esetén nem kell előről kezdeni a futtatást.
Képes több feladat ütemezőt is maga alá gyűrni, de ezt nem teszteltem. Ha csak egy gépen futtatjuk, akkor a job_type értéket local-ra kell állítani, és hiába mondja a dokumentáció, hogy ennyi elég, a submit opciót is be kell állítani. Ez valahol érthető is, hiszen nem mindenki használ Bash-t parancsértelmezőnek.
A genom mérete sajnos kötelező paraméter, ezért ha nincs ötletünk, mekkora lehet a végleges genom mérete, akkor valószínűleg több futtatást is kell végeznünk, ami tekintve a futások hosszát, tetemes időt fog igénybe venni.
Az első lépés tehát a nyers összeszerelés:
fc_run config.file
Ez a lépés három fő fázisból áll. Az elsőben a szekvenciákat egymáshoz illeszti, hogy kiküszöbölje az esetleges szekvenálási hibáka. A háttérben egy adatbázist épít a szekvenciákból. A másodikban ismét egymáshoz illeszti a szekvenciákat, de akkor már megpróbál minél hosszabb szekvenciákat képezni. Ez a két lépés nagyon sokáig tart, érdemes minél több párhozamos folyamatot indítani (njobs paraméter a konfigurációs állományban), mert szálakból csak négyet használ, teljesen mindegy, mennyit definiálunk. Az utolsó lépés a tényleges assembly, ez viszont csak egy program futásából áll, nem számít, mennyi párhuzamos folyamatot állítunk be.
Minden fázisnak egy könyvtár felel meg, ezek a 0-rawreads, 1-preads_ovl, 2-asm-falcon. Az ütemező az első két könyvtárban létrehoz egy build könyvtárat, ami az adott munkafázis mindegyik lépésének további könyvtárakat hoz létre. Sőt, ha a munkafázis párhuzamos folyamatokra osztható, akkor további alkönyvtárak készülnek. Ezekben a könyvtárak mindegyikében van egy run.sh szkript, ami beállítja a környezeti változókat, és elindít egy task.sh szkriptet. Az igazi feladat itt van definiálva. Ezért lehet nyomon követni a program futásának állapotát könyvtárak számolásával. Amennyiben a feladat sikeresen lefutott, létrejön egy run.sh.done fájl.
Az egyetlen problémám a második fázisnál volt, mert a fő szkript hibásan adta át az egyik paramétert. Egészen pontosan az történt, hogy a daligner nevű programnak nem lehet 0 hosszúságú szekvencia vágást engedélyezni. Hiába volt 1 beállítva a fő konfigurációs fájlban, ez nem jutott el valami miatt a hívási sor végén álló bináris programnak. Végül ezt úgy oldottam meg, hogy a build könyvtárban kézzel átírtam az állományt (a beszédes nevű length_cutoff fájlt) a megfelelő értékre.
Ennek a szkriptek-futnak-különböző-könyvtárakban megközelítésnek volt egy hátránya is. Nevezetesen, ha szimbolikus linkeket tartalmaz a könyvtár. Ekkor ugyanis az alkönyvtárban lévő szkript, ha el akar érni valamit egy felsőbb könyvtárban, akkor nem fogja megtalálni, hiszen nem látja, hogy egy szimlinkelt könyvtárban van. Viszont mi nem fogjuk érteni a hibaüzenetet, mert mi meg a szimlinken keresztül ellenőrizzük az elérést, és azt látjuk, hogy minden rendben van.
Az első, nyers összeszerelés a 2-asm-falcon könyvtárban lesz p_ctg.fasta néven. Ha szeretnénk némi statisztikát, mennyire lett jó a munkánk, akkor egy másik PacBio repo, a pb-assembly kell nekünk, abból is a get_asm_stats.py szkript.
python ~/pb-assembly/scripts/get_asm_stats.py 2-asm-falcon/p_ctg.fasta
Ez még csak egy nyers szekvencia lesz. Nyers, mert egy diploid genomnál az apai és anyai kromoszóma készlet között könnyen előfordulhat annyi eltérés, ami megzavarjon egy programot. Ezért a következő lépés ezen eltérések azonosítása, amit a FALCON-unzip végez.
Futtatása ugyancsak egy konfigurációs fájl megadásával történik. A lehetséges paraméterek száma viszont kevesebb, nagy részüket az előző programban már megadtuk. Egyetlen új elem a BAM fájlok elérési útjának beállítása. A dokumentáció elég szűk szavú ezeknek a BAM fájloknak az eredetéről, de ez nem más, mint a szekvenálás eredményeként megkapott subreads.bam (itt lehet még olvasgatni a subreads-ről). A FALCON-unzip futtatása tehát csak ennyi lesz:
fc-unzip.py unzip.cfg
Itt is hiába állítjuk be a szálak számát, mert a program csak egyet fog használni, ezért csak a jobok számával tudjuk szabályozni, mennyi CPU-t hasznosítunk. Az eredmények a 4-polish/cns-output/ könyvtárban lesznek. A cns_p_ctg.fasta fájl fogja tartalmazni a javított szekvenciákat, míg a haplotípusok a cns_h_ctg.fasta fájlban lesznek. Ha megnézzük a statisztikát, azt fogjuk látni, hogy azok lényegesen nem lesznek jobbak, ellenben a szekvenciák megbízhatóbbak lesznek.
Amennyiben van Hi-C adatsorunk is, egy következő lépést is lefuttathatunk, amelynek során tovább javíthatjuk eredményeinket. Ez a falcon-phase, és ugyan csak egy konfigurációs fájllal szabályozhatjuk. Itt kell meadni a Hi-C szekvenciák elérési útját. Én személy szerint nem használtam, és a dokumentáció is elég szegényes, ezért arról nem írnék.
