Kraken: ultragyors metagenomikus szekvencia-osztályozás pontos igazítások alkalmazásával

szekvencia-osztályozási algoritmus

egy DNS-szekvencia s osztályozásához az adott szekvencián belüli összes k-Mer-t egy halmazba gyűjtjük, amelyet K(S) – nek jelölünk. Ezután feltérképezzük az egyes k-mer-eket K(S) – ben, az alábbiakban leírt algoritmus segítségével, az összes Genom LCA taxonjára, amely tartalmazza ezt a k-mer-t. Ezek az LCA taxonok és őseik a taxonómiai fában alkotják azt, amit osztályozási fának nevezünk, egy metszett részfát, amelyet az S osztályozására használnak. Ezután az osztályozási fa minden gyökér-levél (RTL) útvonalát az útvonal mentén lévő összes csomópont súlyának kiszámításával pontozzák. Az osztályozási fában a maximális pontozási RTL útvonal az osztályozási útvonal, az S pedig a levélének megfelelő címkét kapja (ha több maximálisan pontozási út van, akkor az összes út levelének LCA-ja van kiválasztva). Ez az 1. ábrán bemutatott algoritmus lehetővé teszi Kraken számára, hogy egy szekvencián belül minden egyes k-mer-t külön bizonyítéknak tekintsen, majd szükség esetén megpróbálja megoldani az ellentmondó bizonyítékokat. Vegye figyelembe, hogy a megfelelő választáshoz k, a legtöbb k-Mer egyedülállóan egyetlen fajra térképez fel, jelentősen leegyszerűsítve az osztályozási folyamatot. Azok a szekvenciák, amelyeknél a K(S)-ben lévő k-merek egyike sem található meg egyetlen genomban sem, ezt az algoritmust nem osztályozzák.

az RTL útvonal pontozás használata az osztályozási fában szükséges, tekintettel az osztályozandó szekvenciák és a bármely genomkönyvtárban jelen lévő szekvenciák közötti elkerülhetetlen különbségekre. Ilyen különbségek lehetnek, még nagy értékek esetén is k, olyan k-mer-t eredményez, amely jelen van a könyvtárban, de a valódi forrásfajtól távol eső fajhoz kapcsolódik. A különböző RTL útvonalak pontozásával az osztályozási fában kompenzálhatjuk ezeket a különbségeket, és helyesen osztályozhatjuk a szekvenciákat akkor is, ha egy szekvenciában lévő k-Mer kis hányada azt jelzi, hogy a szekvenciához helytelen taxonómiai címkét kell rendelni.

adatbázis létrehozása

a Kraken osztályozási algoritmusának hatékony megvalósítása megköveteli, hogy a K-Mer-ek taxonokhoz történő leképezése egy előre kiszámított adatbázis lekérdezésével történjen. Kraken ezt az adatbázist egy többlépcsős folyamat során hozza létre, kezdve a genomi szekvenciák könyvtárának kiválasztásával. Kraken tartalmaz egy alapértelmezett könyvtárat, amely a Nemzeti Biotechnológiai Információs Központ (NCBI) RefSeq adatbázisában elkészült mikrobiális genomokon alapul, de a könyvtár testreszabható az egyes felhasználók igényei szerint .

a könyvtár kiválasztása után a medúza többszálú k-mer számlálóval létrehozunk egy adatbázist, amely a könyvtár minden különálló 31-mer-jét tartalmazza. Miután az adatbázis elkészült, a 4 bájtos terek medúza tárolására használt k-mer számít az adatbázis fájlban ehelyett Kraken tárolja a taxonómiai azonosító számok a k-Mer’ LCA értékek. Miután a medúza létrehozta az adatbázist, a könyvtárban lévő genomi szekvenciákat egyenként dolgozzák fel. Minden szekvencia esetében a hozzá társított taxon a szekvencia összes k-Mer tárolt LCA-értékének beállítására szolgál. A szekvenciák feldolgozása során, ha egy szekvencia k-mer-jének LCA-értékét korábban beállították, akkor kiszámítják a tárolt érték és az aktuális szekvencia taxonjának LCA-ját, és ezt az LCA-t tárolják a k-mer számára. A taxoninformációkat az NCBI taxonómiai adatbázisából nyerik.

adatbázis-struktúra és keresési algoritmus

mivel a Kraken nagyon gyakran használ k-mer-t adatbázis-lekérdezésként közvetlenül a szomszédos k-mer lekérdezése után, és mivel a szomszédos k-Mer-ek jelentős mennyiségű szekvenciát osztanak meg, a minimalizáló koncepciót használjuk a hasonló k-Mer-ek csoportosítására. Ennek a koncepciónak a magyarázata, itt meghatározzuk az S DNS-szekvencia kanonikus ábrázolását, mint az S lexikográfiailag kisebbét, az S fordított komplementerét. Ahhoz, hogy meghatározzuk az M hosszúság k-mer minimalizálóját, figyelembe vesszük az összes m-Mer kanonikus ábrázolását a k-mer-ben, és a lexikográfiailag legkisebbet választjuk k-mer minimalizálójának. A gyakorlatban a szomszédos k-Mer-ek gyakran ugyanazzal a minimalizálóval rendelkeznek.

Kraken adatbázisában az összes k-Mer azonos minimalizátorral egymás után tárolódik, és kanonikus ábrázolásaik lexikográfiai sorrendjébe vannak rendezve. A K-mer r lekérdezése ezután feldolgozható úgy, hogy egy indexben megkeresi azokat a pozíciókat az adatbázisban, ahol az R minimizerrel rendelkező k-Mer-eket tárolják, majd bináris keresést hajtanak végre az adott régióban (5.ábra). Mivel a szomszédos k-Mer – ek gyakran ugyanazzal a minimalizálóval rendelkeznek, a keresési tartomány gyakran azonos két egymást követő lekérdezés között, és az első lekérdezésben végzett keresés gyakran olyan adatokat hozhat a CPU gyorsítótárába, amelyeket a második lekérdezésben használnak. Azáltal, hogy a későbbi lekérdezések memóriahozzáférései a RAM helyett a CPU gyorsítótárában lévő adatokhoz férhetnek hozzá, ez a stratégia sokkal gyorsabbá teszi a későbbi lekérdezéseket, mint egyébként lennének. Az adatbázisban szereplő k-Mer-csoportok eltolásait tartalmazó indexhez 8 db 4m bájt szükséges. Alapértelmezés szerint a Kraken 15 bp-os minimalizálókat használ, de a felhasználó módosíthatja ezt az értéket; például a minikraken létrehozásakor 13 bp-os minimalizálókat használtunk annak biztosítására, hogy az adatbázis teljes mérete 4 GB alatt maradjon.

Kraken adatbázis szerkezete. Minden lekérdezendő k-mer-nek van egy meghatározott részstringje, amely a minimalizálója. K-mer kereséséhez az adatbázisban megvizsgáljuk az adatbázisban azokat a pozíciókat, amelyek azonos minimalizátorral rendelkező k-Mer-eket tartalmaznak. Ezeket a pozíciókat gyorsan meg lehet találni, ha megvizsgáljuk a minimalizáló eltolási tömböt a rekordok kezdő pozícióihoz a K-mer minimalizálójával (narancssárga) és a következő lehetséges minimalizálóval (kék). Az adott minimalizálóhoz társított rekordok tartományán belül a rekordokat lexikográfiai sorrendjük szerint rendezik k-mers, lehetővé téve a lekérdezés befejezését a bináris keresés ezen a tartományon.

a Kraken implementálása során további optimalizálásokat végeztünk a fent leírt struktúrára és keresési algoritmusra. Első, amint azt Roberts et al. , az M-mers egyszerű lexikográfiai rendezése a minimalizátorok ferde eloszlását eredményezheti, amely túlreprezentálja az alacsony komplexitású M-Mer-eket. Krakenben egy ilyen elfogultság sok nagy keresési tartományt hozna létre, amelyek több időt igényelnének a kereséshez. A minimalizátorok egyenletesebb eloszlásának megteremtéséhez (és ezáltal a keresések felgyorsításához) az exclusive-or (XOR) művelettel váltjuk az egyes m-mer kanonikus reprezentációinak bitjeinek felét, mielőtt az M-Mer-eket lexikográfiai rendezéssel összehasonlítanánk egymással. Ez az XOR művelet hatékonyan megzavarja a standard sorrendet, és megakadályozza a nagy torzítást az alacsony komplexitású minimalizálók felé.

azt is kihasználjuk, hogy a keresési tartomány gyakran megegyezik a lekérdezések között, hogy a Kraken lekérdezései gyorsabbak legyenek. Ahelyett, hogy kiszámítja a minimalizáló minden alkalommal, amikor végre egy lekérdezést, először keresni az előző tartományban. Ha a lekérdezett k-mer ebben a tartományban található, a lekérdezés azonnal visszatérhet. Ha a k-mer nem található, akkor a minimalizáló kiszámításra kerül; ha a k-mer minimizere megegyezik az utoljára lekérdezett k-merével, akkor a lekérdezés sikertelen, mivel a minimalizáló Keresési terében kimutatták, hogy nincs k-mer. Csak akkor, ha a minimalizáló megváltozott, a Krakennek módosítania kell a keresési tartományt, majd újra meg kell keresnie a k-mer-t.

szimulált metagenomok építése

a Hiseq és MiSeq metagenomokat 20 baktérium teljes genomú Vadászpuska leolvasásával építették. Ezeket az olvasásokat vagy a GAGE-B projekt részeként, vagy az NCBI Sequence Read archívumban találták meg. Minden metagenom tíz genom szekvenciáját tartalmazza (további 1.fájl: S1 táblázat). Mind a 10 000, mind a 10 millió leolvasott minta esetében ezeknek a metagenomoknak a szekvenciáinak 10% – át a tíz komponens Genom adatkészletéből választottuk ki (azaz minden genomnak azonos szekvenciabősége volt). Minden szekvenciát levágtunk, hogy eltávolítsuk az alacsony minőségű bázisokat és az adapter szekvenciákat.

e két metagenom összetétele bizonyos kihívásokat jelent osztályozóink számára. Például a Hiseq metagenómunkban található Pelosinus fermentans nem azonosítható helyesen nemzetségi szinten Kraken (vagy bármely más korábban leírt osztályozó), mert a RefSeq teljes genom adatbázisában nincsenek Pelosinus genomok; a Kraken-GB adatbázisában azonban hét ilyen Genom található, köztük a P. fermentans hat törzse. Hasonlóképpen, MiSeq metagenómunkban a Proteus vulgarist gyakran helytelenül osztályozzák a nemzetség szintjén, mert Kraken adatbázisában az egyetlen Proteus Genom egyetlen Proteus mirabilis Genom. További öt Proteus Genom van jelen a Kraken-GB adatbázisában, lehetővé téve a Kraken-GB számára, hogy jobban osztályozza az olvasásokat ebből a nemzetségből. Ezenkívül a MiSeq metagenom öt genomot tartalmaz az Enterobacteriaceae családból (Citrobacter, Enterobacter, Klebsiella, Proteus és Salmonella). Az ebben a családban található nemzetségek közötti nagy szekvenciahasonlóság megnehezítheti a nemzetségek megkülönböztetését bármely osztályozó számára.

a Szimba-5 metagenomot a refseq-ben található teljes bakteriális és archeális genomokból származó leolvasások szimulálásával hozták létre. Azokból a genomokból származó replikonokat akkor használtuk, ha olyan taxonnal társultak, amelynek bejegyzése a nemzetséghez volt társítva rang, ami 607 nemzetségből származó replikonok halmazát eredményezte. Ezután a Mason read szimulátort Illumina modelljével 10 millió 100 bp-os leolvasást állítottunk elő ezekből a genomokból. Először szimulált genomokat hoztunk létre minden fajhoz, 0,1% SNP és 0,1% indel Arány felhasználásával (mindkét alapértelmezett paraméter), amelyekből létrehoztuk az olvasásokat. A szimulált olvasások esetében megszoroztuk az alapértelmezett eltérést és az indel arányokat ötvel, ami átlagosan 2% – os eltérést eredményezett (az olvasás elején 1% – tól a végén 6% – ig terjed), és 1% – os indel arányt (0,5% beillesztési valószínűség és 0,5% deléciós valószínűség). A simBA – 5 metagenome esetében a 10 000 olvasási készletet a 10 millió olvasási készlet véletlenszerű mintájából generálták.

pontosság és sebesség értékelése

úgy döntöttünk, hogy a pontosságot elsősorban a nemzetség szintjén mérjük, amely a legalacsonyabb szint volt, amelyre könnyen meghatározhattuk a phymmbl és az NBC előrejelzéseinek taxonómiai adatait automatizált módon. (Ez annak köszönhető, hogy a phymmbl és az NBC beszámol az eredményekről). Mivel néhány genomnak nincs taxonómiai bejegyzése mind a hét rangban (faj, nemzetség, család, rend, osztály, törzs és Királyság), a nemzetségi szintű érzékenységet a/B-ként határoztuk meg, ahol A a hozzárendelt nemhez tartozó olvasások száma, amelyeket helyesen osztályoztak ebben a rangban, és B a hozzárendelt nemzetséghez tartozó olvasások teljes száma. Hasonlóan határoztuk meg az érzékenységet más taxonómiai rangok esetében is.

mivel a Kraken osztályozhat egy leolvasást a faj feletti szinteken, annak pontosságának mérése megköveteli, hogy meghatározzuk a helyes nemzetség (például) hozzárendelésének pontosságára gyakorolt hatását, miközben egyáltalán nem rendelünk hozzá egy fajt. Ezért a rang-szintű pontosságot C/(D + E) – ként határoztuk meg, ahol C a mért rangnál a helyes taxonon vagy az alatt címkézett olvasások száma, D a mért rangnál vagy az alatt címkézett olvasások száma, E pedig a mért rang felett helytelenül címkézett olvasások száma. Például, ha egy olvasási R – t Escherichia coli-nak kell címkézni, az R E. coli, E. fergusonii vagy Escherichia jelölése javítaná a nemzetség szintű pontosságát. Az Enterobacteriaceae (helyes család) vagy a Proteobaktériumok (helyes törzs) címkéje nem lenne hatással a nemzetség szintű pontosságára. A Bacillus (helytelen nemzetség) vagy Firmicutes (helytelen törzs) r címkéje csökkentené a nemzetség szintű pontosságát.

a PhymmBL pontosságának értékelésekor a fejlesztők tanácsát követve kiválasztottuk a nemzetség megbízhatósági küszöbét összehasonlításainkhoz. A szimulált közepes komplexitású (simmc) adatkészletből 3333 olvasatot választottunk ki, amelyek 31 különböző nemzetséget fednek le. A Simmc készlet Sanger szekvencia adataiból származó rövid olvasások szimulálásához az egyes olvasmányok közül az utolsó 100 bp-t választottuk ki. Ezután lefuttattuk a PhymmBL-t a 100 bp-es leolvasásokkal, és kiértékeltük a phymmbl előrejelzéseinek nemzetségi szintű érzékenységét és pontosságát 0-tól 1-ig terjedő nemzetségi konfidencia küszöbértékekkel, 0,05-ös lépésekben. Megállapítottuk, hogy a 0,65-ös küszöb adta a legmagasabb F-pontszámot (az érzékenység és pontosság harmonikus átlaga), 0,60 és 0,70-nél pedig az F-pontszám is 0-on belül van.A maximális érték 5 százalékpontja (1.kiegészítő fájl: S2 táblázat). Ezért összehasonlításunkban a 0,65 nemzetségi konfidencia küszöböt használtuk. Bár a küszöb kiválasztása a felhasználó egyéni igényeitől függ, és így bizonyos mértékig önkényes, az így kiválasztott küszöb megfelelőbb összehasonlítást nyújt egy olyan szelektív osztályozóval, mint a Kraken, mint egyáltalán nincs küszöb.

a Megablast osztályozóként történő alkalmazásakor kapott idő-és pontossági eredményeket a PhymmBL által előállított naplóadatokból nyertük, mivel a PhymmBL a Megablast-ot használja az igazítási lépéshez. Amikor hozzárendel egy taxonómiai címkét egy olvasási Megablast, használtuk a taxon társított első jelentett igazítás. A Megablast alapértelmezett beállításokkal futott.

a sebességet az egyes programok egyszálú működésével értékeltük (az NBC kivételével). A PhymmBL-t úgy változtatták meg, hogy a blastn programhoz való hívása kettő helyett egy szálat használt. Az NBC-t 36 egyidejű folyamattal futtatták a genomkönyvtárban lévő diszjunkt genomkészleteken, és az osztályozó teljes idejét az egyes genomok dekompressziós és pontozási idejének összegzésével határozták meg. A falióra idejét minden osztályozó számára rögzítették. A Kraken és a többi osztályozó összehasonlításakor a BLAST+ 2.2.27, a PhymmBL 4.0, az NBC 1.1 és a Metaflan 1.7.6-ot használtuk. Az osztályozók mind ugyanazon a számítógépen futottak, 48 AMD Opteron 6172 2,1 GHz-es CPU-val és 252 GB RAM-mal, Red Hat Enterprise Linux 5-tel. A sebességértékeléshez használt adatkészletek mindegyike 10 000 olvasást tartalmazott a Kraken (és annak változatai) és a Metaflan kivételével, amelyek 10 000 000 olvasási adatkészletet használtak. Magasabb olvasási számokat használtak ezekkel a gyorsabb programokkal, hogy minimalizálják a programok végrehajtása során végbemenő kezdeti és végső műveletek hatását.

bár a Kraken az egyetlen az általunk vizsgált programok közül, amely kifejezetten olyan műveleteket hajt végre, amelyek biztosítják, hogy az adatok fizikai memóriában legyenek az osztályozás előtt, biztosak akartunk lenni abban, hogy minden programot hasonló módon értékeltek. A sebesség értékelésekor Minden program esetében elolvassuk az összes adatbázisfájlt (pl. IMM fájlok és BLAST adatbázisok PhymmBL, k-mer frekvencia listák NBC és a Bowtie index Metaflan) a memóriába háromszor futtatása előtt a program, annak érdekében, hogy az adatbázis tartalmát az operációs rendszer cache (amely tárolja a fizikai memóriában).

csökkentett adatbázisméretek

a MiniKraken eredmények 4 GB-os adatbázisának létrehozásához eltávolítottuk a standard Kraken adatbázis minden blokkjának 18 rekordját. A 19-es zsugorodási tényezőt választották, mivel ez volt a legkisebb egész tényező, amely 4 GB-nál kisebbre csökkenti a méretet, ez a méret könnyen beilleszthető sok általános személyi számítógép memóriájába. Azoknak a felhasználóknak, akiknek több RAM áll rendelkezésre, a Kraken lehetővé teszi egy kisebb zsugorodási tényező használatát, ami fokozott érzékenységet biztosít.

vázlat genomok használata

a Kraken-GB adatbázis felépítésekor észrevettük, hogy több olyan kapcsolat volt, amelyek végén ismert adapterszekvenciák voltak. A későbbi tesztekben, azt is megállapítottuk, hogy a nagy mennyiségű emberi szekvenciát tartalmazó minták egyes szekvenciáit ez az adatbázis következetesen tévesen osztályozta, arra enged következtetni, hogy a genomok tervezetében valószínűleg szennyeződés volt jelen. Annak érdekében, hogy ellensúlyozzuk ezt a szennyeződést, eltávolítottuk az adatbázisból azokat a k-Mer-eket az ismert adaptersorozatokból, valamint az első és az utolsó 20 k-Mer-t minden egyes tervezetből. Bár ez javította az osztályozást, nem szüntette meg a téves osztályozás problémáját. Ezért úgy gondoljuk, hogy ha a Kraken adatbázisban vázlatos genomokat használnak, nagyon szigorú intézkedéseket kell alkalmazni a szennyező szekvenciák eltávolítására a genomi könyvtárból.

Kládkizárási kísérletek

a kládkizárási kísérleteink simBA-5 adatkészletének újraelemzésekor néhány leolvasást nem használtak bizonyos mért és kizárt rangpárok esetében. Ha az olvasás eredetének nem volt taxonómiai bejegyzése a mért vagy kizárt rangok egyikében sem, akkor azt nem használták az adott kísérlethez.

ezenkívül egy kísérletben nem használtunk leolvasást, kivéve, ha az adatbázisunkban képviselt legalább két másik taxon (a kizárt kládon kívül) a kizárt rangban megosztotta a származási klád taxonját a mért rangon. Például a G nemzetségből származó olvasmány nem használható olyan kísérletben, amely pontosságot mér az osztály rangján, kivéve a nemzetség rangját, kivéve, ha G otthoni osztályának legalább két másik nemzetsége volt genomokkal Kraken genomikai könyvtárában. E szűrési lépés nélkül, ha egy nemzetséget kizártak volna, amikor ez volt az egyetlen nemzetség az osztályában, Kraken nem tudta megnevezni a megfelelő osztályt, mivel az adatbázis összes bejegyzése ebből az osztályból is kizárásra kerülne. Ez ugyanaz a megközelítés, amelyet hasonló kísérletekben alkalmaztak a PhymmBL értékelésére .

humán mikrobiom adatok osztályozása

a humán mikrobiom projekt adatait egy Kraken adatbázis segítségével osztályoztuk, amely teljes RefSeq bakteriális, archeális és vírusgenomokból, valamint a GRCh37 humán genomból készült. A három csatlakozás (SRS019120, SRS014468 és SRS015055) szekvenciáit az NCBI Sequence Read archívumából szereztük be, és minden csatlakozásnak két futása volt. Minden olvasást levágtak, hogy eltávolítsák az alacsony minőségű bázisokat és az adapter szekvenciákat. A Krona-t az összes taxonómiai eloszlási terület előállításához használták.

mivel a szekvenciák mind párosított olvasmányok voltak, úgy kötöttük össze az olvasásokat, hogy összefűztük a társakat egy ‘NNNNN’ szekvenciával közöttük. Kraken figyelmen kívül hagyja a kétértelmű nukleotidokkal rendelkező k-Mer-eket, így az ezeket az ‘N’ karaktereket átfogó k-Mer-ek nem befolyásolják az osztályozást. Ez a művelet lehetővé tette Kraken számára, hogy egy olvasópárt egyetlen egységként osztályozzon, ahelyett, hogy külön kellene osztályoznia a társakat.

szoftverek és adatok elérhetősége