Kraken: ultrarychlá klasifikace metagenomické sekvence pomocí exact alignments

sekvenční klasifikační algoritmus

pro klasifikaci sekvence DNA shromažďujeme všechny k-Mery v této sekvenci do sady označené jako K (S). Poté mapujeme každý k-mer v K(S) pomocí algoritmu popsaného níže na taxon LCA všech genomů, které obsahují tento K-mer. Tyto taxony LCA a jejich předkové ve stromu taxonomie tvoří to, co nazýváme klasifikační strom, prořezaný podstrom, který se používá ke klasifikaci s. každý uzel v klasifikačním stromu je vážen počtem k-mers v K (Y), které jsou mapovány na taxon spojený s tímto uzlem. Poté je každá cesta root-to-leaf (RTL) v klasifikačním stromu hodnocena výpočtem součtu všech hmotností uzlů podél cesty. Maximální bodovací RTL cesta v klasifikačním stromu je klasifikační cesta A S je přiřazen štítek odpovídající jeho listu (pokud existuje více maximálně bodovacích cest, je vybrána LCA listů všech těchto cest). Tento algoritmus, znázorněno na obrázku 1, umožňuje Kraken zvážit každý k-mer v sekvenci jako samostatný důkaz, a v případě potřeby se pokuste vyřešit jakékoli protichůdné důkazy. Všimněte si, že pro vhodnou volbu k, Většina k-mers bude mapovat jedinečně na jeden druh, což výrazně zjednoduší klasifikační proces. Sekvence, pro které se žádný z k-merů v K(Y) nenachází v žádném genomu, jsou tímto algoritmem neklasifikovány.

použití hodnocení cest RTL v klasifikačním stromu je nezbytné s ohledem na nevyhnutelné rozdíly mezi sekvencemi, které mají být klasifikovány, a sekvencemi přítomnými v jakékoli knihovně genomů. Takové rozdíly mohou, i pro velké hodnoty k, mít za následek k-mer, který je přítomen v knihovně, ale je spojen s druhem, který je daleko od skutečného zdrojového druhu. Bodováním různých cest RTL ve klasifikačním stromu, můžeme tyto rozdíly kompenzovat a správně klasifikovat sekvence, i když malá menšina k-mers v sekvenci naznačuje, že sekvenci by mělo být přiřazeno nesprávné taxonomické označení.

vytvoření databáze

efektivní implementace krakenova klasifikačního algoritmu vyžaduje, aby mapování k-mers na taxony bylo provedeno dotazováním na předem vypočítanou databázi. Kraken vytváří tuto databázi Pomocí vícestupňového procesu, počínaje výběrem knihovny genomických sekvencí. Kraken obsahuje výchozí knihovnu založenou na dokončených mikrobiálních genomech v databázi RefSeq Národního centra pro biotechnologické informace (NCBI), ale knihovnu lze přizpůsobit podle potřeby jednotlivých uživatelů .

jakmile je knihovna vybrána, použijeme čítač medúzy multithreaded k – mer k vytvoření databáze obsahující každý odlišný 31-mer v knihovně. Jakmile je databáze je kompletní, 4-byte mezery medúzy používané k ukládání k-mer počty v databázovém souboru jsou místo toho používá Kraken k ukládání taxonomických ID čísla hodnot LCA k-mers‘. Po vytvoření databáze medúzy se genomické sekvence v knihovně zpracovávají jeden po druhém. Pro každou sekvenci se taxon s ní spojený používá k nastavení uložených hodnot LCA všech k-mers v sekvenci. Při zpracování sekvencí, pokud má k-mer ze sekvence dříve nastavenou hodnotu LCA, vypočítá se LCA uložené hodnoty a taxonu aktuální sekvence a tato LCA je uložena pro K-mer. Informace o taxonu jsou získávány z databáze taxonomie NCBI.

struktura databáze a vyhledávací algoritmus

protože Kraken velmi často používá k-mer jako databázový dotaz bezprostředně po dotazování na sousední k-mer, a protože sousední k-Mer sdílejí značné množství sekvence, využíváme koncept minimizer ke seskupení podobných k-Mer dohromady. Vysvětlit naši aplikaci tohoto konceptu, zde definujeme kanonickou reprezentaci sekvence DNA s jako lexikograficky menší S a reverzní doplněk s. Abychom určili k-merův minimalizátor délky M, považujeme kanonickou reprezentaci všech M-Mer v k-mer, a vyberte lexikograficky nejmenší z těchto M-Mer jako minimalizátor k-mer. V praxi budou mít sousední k-mers často stejný minimalizátor.

v krakenově databázi jsou všechny k-Mery se stejným minimalizátorem uloženy postupně a jsou řazeny v lexikografickém pořadí jejich kanonických reprezentací. Dotaz na K-mer R pak může být zpracován vyhledáním v indexu pozic v databázi, kde by byly uloženy k-mers s minimalizátorem R, a poté provedením binárního vyhledávání v této oblasti (obrázek 5). Protože sousední k-Mery mají často stejný minimalizátor, rozsah vyhledávání je často stejný mezi dvěma po sobě jdoucími dotazy a vyhledávání v prvním dotazu může často přinést data do mezipaměti CPU, která budou použita ve druhém dotazu. Tím, že umožňuje přístup k paměti v následných dotazech pro přístup k datům v mezipaměti CPU místo paměti RAM, tato strategie umožňuje následné dotazy mnohem rychleji, než by jinak byly. Index obsahující offsety každé skupiny k-mers v databázi vyžaduje 8 × 4M bajtů. Ve výchozím nastavení Kraken používá 15-BP minimalizátory, ale uživatel může tuto hodnotu upravit; například při vytváření MiniKraken jsme použili 13-BP minimalizátory, abychom zajistili, že celková velikost databáze zůstane pod 4 GB.

struktura databáze Kraken. Každý k-mer, který má být dotazován proti databázi, má specifický podřetězec, který je jeho minimalizátorem. Chcete-li vyhledat k-mer v databázi, jsou zkoumány pozice v databázi, které obsahují k-mers se stejným minimalizátorem. Tyto pozice jsou rychle nalezeny zkoumáním pole minimizer offset pro počáteční pozice záznamů s minimalizátorem k-mer (oranžový) a dalším možným minimalizátorem (modrý). V rámci rozsahu záznamů spojených s daným minimalizátorem jsou záznamy seřazeny podle lexikografického uspořádání jejich k-mers, což umožňuje vyplnění dotazu pomocí binárního vyhledávání v tomto rozsahu.

při implementaci Krakenu jsme provedli další optimalizace struktury a vyhledávacího algoritmu popsaného výše. Za prvé, jak poznamenal Roberts et al. , jednoduché lexikografické uspořádání M-mers může mít za následek zkosené rozdělení minimalizátorů, které nadměrně představuje m-Mers s nízkou složitostí. V Krakenu by takové zkreslení vytvořilo mnoho velkých rozsahů vyhledávání, což by vyžadovalo více času na vyhledávání. Abychom vytvořili rovnoměrnější distribuci minimalizátorů (a tím urychlili vyhledávání), použijeme operaci exclusive-or (XOR) k přepínání poloviny bitů kanonické reprezentace každého M-mer před vzájemným porovnáním M-Mer pomocí lexikografického uspořádání. Tato operace XOR účinně zakóduje standardní uspořádání a zabraňuje velkému zkreslení směrem k minimalizátorům s nízkou složitostí.

využíváme také skutečnosti, že rozsah vyhledávání je mezi dotazy často stejný, aby se krakenovy dotazy zrychlily. Spíše než vypočítat minimalizátor pokaždé, když provedeme dotaz, nejprve prohledáme předchozí rozsah. Pokud se náš dotazovaný k-mer nachází v tomto rozsahu, může se dotaz okamžitě vrátit. Pokud k-mer není nalezen, pak se vypočítá minimalizátor; pokud je minimalizátor k-mer stejný jako poslední dotazovaný k-mer, pak dotaz selže, protože se ukázalo, že vyhledávací prostor minimizeru nemá k-mer. Pouze v případě, že se minimalizátor změnil, musí Kraken upravit rozsah vyhledávání a znovu vyhledat k-mer.

konstrukce simulovaných metagenomů

metagenomy HiSeq a MiSeq byly postaveny za použití 20 sad bakteriálních celogenomových brokovnic. Tato čtení byla nalezena buď jako součást projektu GAGE-B, nebo v archivu čtení sekvence NCBI. Každý metagenom obsahuje sekvence z deseti genomů (další soubor 1: tabulka S1). Pro 10 000 a 10 milionů čtených vzorků každého z těchto metagenomů bylo vybráno 10% jejich sekvencí z každé z deseti datových sad genomu (tj. Všechny sekvence byly oříznuty, aby se odstranily základny nízké kvality a sekvence adaptérů.

složení těchto dvou metagenomů představuje pro naše klasifikátory určité výzvy. Například pelosinus fermentans, nalezený v našem hiseq metagenomu, nemůže být správně identifikován na úrovni rodu Krakenem (nebo jiným dříve popsaným klasifikátorem), protože v databázi úplných genomů RefSeq nejsou žádné genomy Pelosinus; v databázi Kraken-GB je však sedm takových genomů, včetně šesti kmenů P.fermentans. Podobně, v našem MiSeq metagenomu, Proteus vulgaris je často klasifikován nesprávně na úrovni rodu, protože jediný Proteus genom v krakenově databázi je jediný Proteus mirabilis genom. V databázi Kraken-GB je přítomno dalších pět Proteusových genomů, což umožňuje Kraken-GB lépe klasifikovat čtení z tohoto rodu. Kromě toho Metagenom MiSeq obsahuje pět genomů z čeledi Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus a Salmonella). Vysoká sekvenční podobnost mezi rody v této rodině může ztížit rozlišení mezi rody pro každého klasifikátora.

metagenom simBA-5 byl vytvořen simulací čtení ze sady kompletních bakteriálních a archaeálních genomů v RefSeq. Replikony z těchto genomů byly použity, pokud byly spojeny s taxonem, který měl záznam spojený s rodovou hodností, což mělo za následek sadu replikonů z 607 rodů. Poté jsme použili Mason read simulator s jeho modelem Illumina k výrobě 10 milionů 100-bp čtení z těchto genomů. Nejprve jsme vytvořili simulované genomy pro každý druh, pomocí rychlosti SNP 0,1% a rychlosti indel 0,1% (oba výchozí parametry), ze kterých jsme generovali čtení. Pro simulované čtení jsme vynásobili výchozí míry nesouladu a indelu pěti, což mělo za následek průměrnou míru nesouladu 2% (v rozmezí od 1% na začátku čtení do 6% na koncích) a míru indelu 1% (0,5% pravděpodobnost Vložení a 0,5% pravděpodobnost vymazání). Pro metagenom simBA-5 byla sada 10 000 čtení generována z náhodného vzorku sady 10 milionů čtení.

vyhodnocení přesnosti a rychlosti

rozhodli jsme se měřit přesnost především na úrovni rodu, což byla nejnižší úroveň, pro kterou jsme mohli snadno automatizovaně určit informace o taxonomii pro předpovědi PhymmBL a NBC. (Je to kvůli způsobu, jakým PhymmBL a NBC hlásí své výsledky). Protože některé genomy nemají taxonomické záznamy na všech sedmi řadách (druh, rod, rodina ,řád, třída, kmen a království), definovali jsme citlivost na úrovni rodu jako A / B, kde A je počet čtení s přiřazeným rodem, které byly správně klasifikovány v této hodnosti, a B je celkový počet čtení s jakýmkoli přiřazeným rodem. Citlivost jsme definovali podobně pro ostatní taxonomické řady.

protože Kraken může klasifikovat čtení na úrovních nad druhem, měření jeho přesnosti vyžaduje, abychom definovali vliv na přesnost přiřazení správného rodu (například), zatímco vůbec nepřiřazujeme druh. Z tohoto důvodu jsme definovali přesnost na úrovni hodnosti jako C/(D + E), kde C je počet čtení označených na nebo pod správným taxonem na měřené hodnosti, D je počet čtení označených na nebo pod měřenou hodností a E je počet čtení nesprávně označených nad měřenou hodností. Například při čtení R, které by mělo být označeno jako Escherichia coli, by označení R jako E. coli, E. fergusonii nebo Escherichia zlepšilo přesnost na úrovni rodu. Označení Enterobacteriaceae (správná rodina) nebo Proteobacteria (správný kmen) by nemělo žádný vliv na přesnost na úrovni rodu. Označení pro R Bacillus (nesprávný rod) nebo Firmicutes (nesprávný kmen) by snížilo přesnost na úrovni rodu.

při hodnocení přesnosti PhymmBL jsme na základě doporučení vývojářů vybrali prahovou hodnotu spolehlivosti rodu pro naše Srovnání. Vybrali jsme 3,333 čtení z datové sady simulované střední složitosti (simMC), pokrývající 31 různých rodů. Pro simulaci krátkých čtení z dat sangerovy sekvence v sadě simMC jsme vybrali posledních 100 bp z každého čtení. Poté jsme provedli PhymmBL proti těmto 100-bp čtením a vyhodnotili citlivost na úrovni rodu a přesnost předpovědí PhymmBL s prahovými hodnotami spolehlivosti rodu od 0 do 1, v krocích po 0,05. Zjistili jsme, že práh 0,65 přinesl nejvyšší F-skóre (harmonický průměr citlivosti a přesnosti), s 0,60 a 0,70 také s F-skóre v rámci 0.5 procentních bodů maxima (další soubor 1: tabulka S2). V našem srovnání jsme proto použili práh spolehlivosti rodu 0, 65. Ačkoli výběr prahové hodnoty závisí na individuálních potřebách uživatele, a proto je do jisté míry libovolný, takto zvolená prahová hodnota poskytuje správnější srovnání se selektivním klasifikátorem, jako je Kraken, než vůbec žádný práh.

výsledky času a přesnosti při použití Megablastu jako klasifikátoru byly získány z dat protokolu produkovaných PhymmBL, protože PhymmBL používá Megablast pro svůj krok zarovnání. Při přiřazování taxonomického označení ke čtení pomocí Megablastu jsme použili taxon spojený s prvním hlášeným zarovnáním. Megablast byl spuštěn s výchozími možnostmi.

rychlost byla hodnocena pomocí jednovláknové operace každého programu (s výjimkou NBC). PhymmBL byl změněn tak, že jeho volání do programu blastn používalo jedno vlákno místo dvou. NBC byl spuštěn 36 souběžné procesy pracující na disjunktních sadách genomů ve své genomové knihovně, a celkový čas pro klasifikátor byl určen součtem dekompresních a bodovacích časů pro každý genom. Pro všechny klasifikátory byly zaznamenány časy nástěnných hodin. Při porovnání Krakenu s ostatními klasifikátory jsme použili BLAST + 2.2.27, PhymmBL 4.0, NBC 1.1 a MetaPhlAn 1.7.6. Klasifikátory byly všechny spuštěny na stejném počítači, s 48 procesory AMD Opteron 6172 2.1 GHz a 252 GB RAM, se systémem Red Hat Enterprise Linux 5. Datové sady používané pro vyhodnocení rychlosti měly 10 000 čtení pro všechny programy jiné než Kraken (a jeho varianty) a Metaflan, které používaly 10 000 000 čtení datových sad. Vyšší počty čtení byly použity s těmito rychlejšími programy, aby se minimalizoval účinek počátečních a konečných operací, ke kterým dochází během provádění programů.

ačkoli Kraken je jediný z programů, které jsme zkoumali a které explicitně provádějí operace, aby zajistily, že jeho data jsou ve fyzické paměti před klasifikací, chtěli jsme si být jisti, že všechny programy byly vyhodnoceny podobným způsobem. Při vyhodnocování rychlosti pro každý program čteme všechny databázové soubory (např. IMM soubory a BLAST databáze pro PhymmBL, k-Mer frekvenční seznamy pro NBC a Bowtie index pro Metaflan) do paměti třikrát před spuštěním programu, aby se obsah databáze v mezipaměti operačního systému (který je uložen ve fyzické paměti).

snížené velikosti databáze

abychom vygenerovali databázi 4 GB pro naše výsledky MiniKraken, odstranili jsme prvních 18 z každého bloku 19 záznamů ve standardní databázi Kraken. Byl vybrán zmenšující se faktor 19, protože to byl nejmenší celočíselný faktor, který by zmenšil velikost na méně než 4 GB, což je velikost, která se snadno vejde do paměti mnoha běžných osobních počítačů. Pro uživatele, kteří mají k dispozici více paměti RAM, Kraken umožňuje použití menšího smršťovacího faktoru, což zvýší citlivost.

použití návrhových genomů

při konstrukci databáze Kraken-GB jsme si všimli, že na koncích je několik spojů se známými adaptérovými sekvencemi. V následných testech, také jsme zjistili, že některé sekvence ve vzorcích s velkým množstvím lidské sekvence byly touto databází důsledně nesprávně klasifikovány, což nás vede k závěru, že kontaminace byla pravděpodobně přítomna v návrzích genomů. Ve snaze zabránit této kontaminaci, odstranili jsme z databáze tyto k-Mery ze známých sekvencí adaptérů, stejně jako první a poslední 20 k-mers z každého návrhu kontigů. I když to zlepšilo klasifikaci, neodstranilo to problém s nesprávným klasifikací. Z tohoto důvodu se domníváme, že pokud jsou v databázi Kraken použity návrhy genomů, měla by být použita velmi přísná opatření k odstranění kontaminačních sekvencí z genomické knihovny.

experimenty s vyloučením Clade

při opětovné analýze datové sady simBA-5 pro naše experimenty s vyloučením clade nebyla některá čtení použita pro určité páry měřených a vyloučených řad. Pokud původ čtení postrádal taxonomický záznam v jedné z měřených nebo vyloučených řad, nebyl použit pro tento konkrétní experiment.

kromě toho nebylo v experimentu použito čtení, pokud alespoň dva další taxony zastoupené v naší databázi (kromě vyloučeného clade) ve vyloučené hodnosti nesdílely taxon clade of origin v měřené hodnosti. Například, čtení z rodu G by nebylo použito v experimentu, který měří přesnost v hodnosti třídy a vylučuje hodnost rodu, pokud domácí třída G neměla v krakenově genomové knihovně alespoň dva další rody s genomy. Bez tohoto kroku filtrování, byl vyloučen Rod, když to byl jediný rod ve své třídě, Kraken nemohl pojmenovat správnou třídu, protože všechny položky v databázi z této třídy by byly také vyloučeny. Jedná se o stejný přístup přijatý v podobných experimentech, které byly použity k vyhodnocení PhymmBL .

klasifikace údajů o lidském mikrobiomu

klasifikovali jsme data Projektu lidského mikrobiomu pomocí databáze Kraken vytvořené z kompletních bakteriálních, archaeálních a virových genomů RefSeq spolu s lidským genomem GRCh37. Získali jsme sekvence tří přístupů (SRS019120, SRS014468 a SRS015055) z archivu čtení sekvence NCBI a každé přistoupení mělo předloženy dva běhy. Všechna čtení byla oříznuta, aby se odstranily základny nízké kvality a sekvence adaptérů. Krona byla použita k vytvoření všech taxonomických distribučních pozemků.

protože všechny sekvence byly spárovány, spojili jsme čtení dohromady spojením matů se sekvencí „NNNNN“ mezi nimi. Kraken ignoruje k-mers s dvojznačnými nukleotidy, takže k-mers, které pokrývají tyto‘ N ‚ znaky, neovlivňují klasifikaci. Tato operace umožnila Kraken klasifikovat dvojici čtení jako jednu jednotku, spíše než aby museli klasifikovat kamarády Samostatně.

dostupnost softwaru a dat