Kraken: ultrasnelle metagenomic sequence classification using exact alignments

Sequence classification algorithm

om een DNA sequentie S te classificeren, verzamelen we alle k-mers binnen die sequentie in een verzameling, aangeduid als K(S). Vervolgens brengen we elke k-mer in K(S) in kaart, met behulp van het algoritme dat hieronder wordt beschreven, aan de LCA-taxon van alle genomen die die k-mer bevatten. Deze LCA-taxa en hun voorouders in de taxonomieboom vormen wat we de classificatieboom noemen, een gesnoeide subboom die wordt gebruikt om S. te classificeren. elk knooppunt in de classificatieboom wordt gewogen met het aantal k-mers in K(S) dat is toegewezen aan het taxon geassocieerd met dat knooppunt. Vervolgens wordt elk root-to-leaf (RTL) pad in de classificatieboom gescoord door de som van alle knooppuntgewichten langs het pad te berekenen. Het maximale scorebord RTL-pad in de classificatiestructuur is het classificatiepad, en S krijgt het label toegewezen dat overeenkomt met het blad (als er meerdere maximale scoreborden zijn, wordt de LCA van de bladeren van al die paden geselecteerd). Dit algoritme, geïllustreerd in Figuur 1, staat Kraken toe om elke k-mer binnen een sequentie te beschouwen als een afzonderlijk bewijsstuk, en vervolgens te proberen om eventueel tegenstrijdig bewijs op te lossen. Merk op dat Voor een geschikte keuze van k, de meeste k-mers uniek in kaart zullen brengen op een enkele soort, waardoor het classificatieproces aanzienlijk wordt vereenvoudigd. De opeenvolgingen waarvoor geen van de k-mers in K(S) in om het even welk genoom worden gevonden worden niet geclassificeerd door dit algoritme verlaten.

het gebruik van RTL-padscore in de classificatieboom is noodzakelijk in het licht van de onvermijdelijke verschillen tussen de te classificeren sequenties en de sequenties die aanwezig zijn in een bibliotheek van genomen. Dergelijke verschillen kunnen, zelfs voor grote waarden van k, resulteren in een k-mer die in de bibliotheek aanwezig is, maar geassocieerd wordt met een soort die ver verwijderd is van de echte bronsoort. Door de verschillende RTL-paden in de classificatieboom te noteren, kunnen we deze verschillen compenseren en sequenties correct classificeren, zelfs wanneer een kleine minderheid van k-mers in een sequentie aangeeft dat de sequentie een onjuist taxonomisch label moet krijgen.

databasecreatie

efficiënte implementatie van Kraken ’s classificatiealgoritme vereist dat de mapping van k-MER’ s naar taxa wordt uitgevoerd door een vooraf berekende database te bevragen. Kraken creëert deze database door middel van een multi-stap proces, te beginnen met de selectie van een bibliotheek van genomische sequenties. Kraken bevat een standaardbibliotheek, gebaseerd op voltooide microbiële genomen in de RefSeq database van het National Center for Biotechnology Information (NCBI), maar de bibliotheek kan worden aangepast als dat nodig is door individuele gebruikers .

zodra de bibliotheek is gekozen, gebruiken we de kwal multithreaded K-Mer-teller om een database aan te maken die elke afzonderlijke 31-Mer in de bibliotheek bevat. Zodra de database is voltooid, worden de 4-byte spaties kwallen gebruikt om de k-mer tellingen in het database bestand op te slaan in plaats daarvan gebruikt door Kraken om de taxonomische ID nummers van de K-mers’ LCA waarden op te slaan. Nadat de database door kwallen is gemaakt, worden de genomische opeenvolgingen in de bibliotheek één voor één verwerkt. Voor elke reeks wordt het daaraan gekoppelde taxon gebruikt om de opgeslagen LCA-waarden van alle k-mers in de reeks in te stellen. Als een k-mer van een reeks zijn LCA-waarde eerder heeft ingesteld, wordt de LCA van de opgeslagen waarde en het taxon van de huidige reeks berekend en wordt die LCA opgeslagen voor de k-Mer. Taxon informatie wordt verkregen uit de NCBI taxonomy database.

databasestructuur en zoekalgoritme

omdat Kraken zeer vaak een k-mer gebruikt als een database query onmiddellijk na het bevragen van een aangrenzende k-mer, en omdat aangrenzende k-mers een aanzienlijke hoeveelheid sequentie delen, gebruiken we het minimizer concept om soortgelijke k-mers samen te groeperen. Om onze toepassing van dit concept te verklaren, definiëren we hier de canonieke representatie van een DNA-sequentie S als de lexicografisch kleinere van S en de omgekeerde complement van S. Om een K-mer ’s minimizer van lengte M te bepalen, beschouwen we de canonieke representatie van alle M-mers in de k-mer, en selecteren we de lexicografisch kleinste van die m-mers als de k-MER’ s minimizer. In de praktijk zullen aangrenzende k-mers vaak dezelfde minimizer hebben.

in de database van Kraken worden alle k-mers met dezelfde minimizer opeenvolgend opgeslagen en gesorteerd in lexicografische volgorde van hun canonieke representaties. Een query voor een k-mer R kan dan worden verwerkt door op te zoeken in een index de posities in de database waar de k-mers met R ‘ S minimizer zou worden opgeslagen, en dan het uitvoeren van een binaire zoekopdracht binnen die regio (Figuur 5). Omdat aangrenzende k-mers vaak dezelfde minimizer hebben, is het zoekbereik vaak hetzelfde tussen twee opeenvolgende query ‘ s, en de zoekopdracht in de eerste query kan vaak gegevens in de CPU-cache brengen die in de tweede query zal worden gebruikt. Door het toestaan van geheugen toegangen in de volgende query ’s om toegang te krijgen tot gegevens in de CPU-cache in plaats van RAM, deze strategie maakt de volgende query’ s veel sneller dan ze anders zou zijn. De index met de offsets van elke groep k-mers in de database vereist 8 × 4M bytes. Kraken gebruikt standaard 15-bp minimizers, maar de gebruiker kan deze waarde wijzigen; bijvoorbeeld, bij het maken van MiniKraken, gebruikten we 13-bp minimizers om ervoor te zorgen dat de totale databasegrootte onder de 4 GB bleef.

Kraken database structuur. Elke k-mer die wordt opgevraagd tegen de database heeft een specifieke substring die de minimizer is. Om te zoeken naar een k-mer in de database, worden de posities in de database die k-mers bevatten met dezelfde minimizer onderzocht. Deze posities worden snel gevonden door de minimizer offset array te onderzoeken voor de startposities van records met de k-mer ‘ s minimizer (oranje) en de volgende mogelijke minimizer (blauw). Binnen een reeks records die geassocieerd zijn met een bepaalde minimizer, worden records gesorteerd op lexicografische volgorde van hun k-mers, waardoor een query kan worden ingevuld met behulp van een binaire zoekopdracht over dit bereik.

bij de implementatie van Kraken hebben we verdere optimalisaties gemaakt voor de hierboven beschreven structuur-en zoekalgoritme. Ten eerste, zoals opgemerkt door Roberts et al. , een eenvoudige lexicografische ordening van M-mers kan resulteren in een scheef verdeelde verdeling van minimizers die over-vertegenwoordigt lage complexiteit M-mers. In Kraken zou zo ‘ n bias veel grote zoekbereiken creëren, wat meer tijd nodig zou hebben om te zoeken. Om een gelijkmatiger verdeling van minimizers te creëren (en zo zoekopdrachten te versnellen), gebruiken we de exclusieve-of (XOR) operatie om de helft van de bits van de canonieke representatie van elke M-mer om te schakelen voordat we de M-mers met elkaar vergelijken met behulp van lexicografische ordening. Deze XOR operatie scramble effectief de standaard bestelling, en voorkomt de grote bias naar lage complexiteit minimizers.

we maken ook gebruik van het feit dat het zoekbereik vaak hetzelfde is tussen query ’s om Kraken’ s query ‘ s sneller te maken. In plaats van het berekenen van de minimizer elke keer dat we een query uit te voeren, we eerst zoeken in het vorige bereik. Als onze queried k-mer in dit bereik wordt gevonden, kan de query onmiddellijk terugkeren. Als de k-mer niet wordt gevonden, dan wordt de minimizer berekend; als de minimizer van de k-mer hetzelfde is als de laatst opgevraagde k-mer ‘ s, dan mislukt de query, omdat de zoekruimte van de minimizer niet de k-Mer heeft. Alleen als de minimizer is veranderd, moet Kraken het zoekbereik aanpassen en opnieuw zoeken naar de k-mer.

construeren van gesimuleerde metagenomen

de HiSeq en MiSeq metagenomen werden gebouwd met behulp van 20 sets van bacteriële shotgun reads. Deze reads werden gevonden Als onderdeel van het GAGE-B project of in het NCBI Sequence Read Archive. Elk metagenoom bevat opeenvolgingen van tien genomen (aanvullend dossier 1: tabel S1). Voor zowel 10.000 als 10 miljoen gelezen steekproeven van elk van deze metagenomes, werden 10% van hun opeenvolgingen geselecteerd uit elk van de tien reeksen van de gegevens van het componentengenoom (d.w.z., had elk genoom gelijke opeenvolgingsovervloed). Alle opeenvolgingen werden bijgesneden om lage kwaliteitsbasissen en adapteropvolgingen te verwijderen.

de samenstelling van deze twee metagenomen stelt bepaalde uitdagingen voor onze classificeerders. Bijvoorbeeld, Pelosinus fermentans, gevonden in onze HiSeq Metagenome, kan niet correct worden geïdentificeerd op genus niveau door Kraken (of een van de andere eerder beschreven classifiers), omdat er geen pelosinus genomen in de RefSeq complete genomen database; Er zijn echter zeven van dergelijke genomen in Kraken-GB ‘ s database, waaronder zes stammen van P. fermentans. Op dezelfde manier wordt Proteus vulgaris in ons MiSeq-metagenoom vaak onjuist geclassificeerd op genus-niveau omdat het enige Proteus-genoom in de database van Kraken een enkel Proteus mirabilis-genoom is. Er zijn nog vijf Proteus genomen aanwezig in de database van Kraken-GB, waardoor Kraken-GB reads van dat geslacht beter kan classificeren. Daarnaast bevat het MiSeq Metagenoom vijf genomen uit de Enterobacteriaceae familie (Citrobacter, Enterobacter, Klebsiella, Proteus en Salmonella). De hoge sequentie gelijkenis tussen de geslachten in deze familie kan onderscheid maken tussen geslachten moeilijk voor elke classifier.

het simBA-5 metagenoom werd gecreëerd door het simuleren van reads uit de set van complete bacteriële en archaeale genomen in RefSeq. Replicons van die genomen werden gebruikt als ze werden geassocieerd met een taxon dat een vermelding geassocieerd met de genus rang had, resulterend in een set van replicons van 607 geslachten. Vervolgens gebruikten we de Mason read simulator met zijn Illumina model om 10 miljoen 100-bp reads te produceren uit deze genomen. Eerst creëerden we gesimuleerde genomen voor elke soort, met behulp van een SNP-tarief van 0,1% en een indel-tarief van 0,1% (beide standaardparameters), waaruit we de reads gegenereerd. Voor de gesimuleerde leest, hebben we de standaard mismatch en indel tarieven vermenigvuldigd met vijf, wat resulteert in een gemiddelde mismatch tarief van 2% (variërend van 1% Aan het begin van Leest tot 6% aan de uiteinden) en een indel tarief van 1% (0,5% insertie waarschijnlijkheid en 0,5% deletie waarschijnlijkheid). Voor simBA-5 metagenome, werd de 10.000 gelezen reeks gegenereerd uit een willekeurige steekproef van de 10 miljoen gelezen reeks.

evaluatie van nauwkeurigheid en snelheid

we kozen ervoor om nauwkeurigheid te meten op het genusniveau, het laagste niveau waarvoor we gemakkelijk de taxonomie-informatie voor PhymmBL en NBC ‘ s voorspellingen op geautomatiseerde wijze konden bepalen. (Dit komt door de manier waarop PhymmBL en NBC hun resultaten rapporteren). Omdat sommige genomen geen taxonomische gegevens hebben op alle zeven rangen (soorten, genus, familie, orde, klasse, phylum en Koninkrijk), definieerden we genus-niveau gevoeligheid als A/B, waar A het aantal reads is met een toegewezen genus die correct werden geclassificeerd op die rang, en B het totale aantal reads met een toegewezen genus is. We definieerden gevoeligheid op dezelfde manier voor andere taxonomische rangen.

omdat Kraken een read kan classificeren op niveaus boven de soort, vereist het meten van de precisie ervan dat we het effect op de precisie van het toewijzen van het juiste geslacht (bijvoorbeeld) definiëren, terwijl het helemaal niet toewijzen van een soort. Om deze reden hebben we rang-niveau precisie gedefinieerd als C/(D + E), waarbij C is het aantal reads gelabeld op of onder het juiste taxon bij de gemeten rang, D is het aantal reads gelabeld op of onder de gemeten rang, en E is het aantal reads onjuist gelabeld boven de gemeten rang. Bijvoorbeeld, gegeven een gelezen R die als Escherichia coli zou moeten worden geëtiketteerd, zou een etikettering van R als E. coli, E. fergusonii of Escherichia genus-niveau precisie verbeteren. Een label van Enterobacteriaceae (correcte familie) of Proteobacteria (correcte phylum) zou geen effect hebben op de nauwkeurigheid van het geslachtsniveau. Een label voor R van Bacillus (incorrect genus) of Firmicutes (incorrect phylum) zou de nauwkeurigheid van het genusniveau verminderen.

bij het evalueren van PhymmBL’ s nauwkeurigheid, op advies van de ontwikkelaars , hebben we een genus-betrouwbaarheidsdrempel gekozen voor onze vergelijkingen. We selecteerden 3.333 reads uit de simulated medium complexity (simMC) dataset, die 31 verschillende geslachten besloeg. Om korte leest van de sanger-opeenvolgingsgegevens in de simMC-reeks te simuleren, selecteerden wij laatste 100 bp van elk leest. Vervolgens hebben we PhymmBL vergeleken met die 100-bp reads, en de gevoeligheid en precisie van PhymmBL ‘ s voorspellingen geëvalueerd met genus betrouwbaarheidsdrempels van 0 tot 1, in stappen van 0,05. We vonden dat een drempel van 0,65 de hoogste F-score opleverde (het harmonische gemiddelde van gevoeligheid en precisie), met 0.60 en 0.70 ook F-scores binnen 0.5 procentpunten van het maximum (aanvullend dossier 1: tabel S2). Daarom hebben we de 0,65 genus betrouwbaarheidsdrempel gebruikt in onze vergelijkingen. Hoewel de keuze van een drempel afhankelijk is van de individuele behoeften van een gebruiker, en dit in zekere mate willekeurig is, biedt een op deze manier geselecteerde drempel een betere vergelijking met een selectieve classificeerder zoals Kraken dan helemaal geen drempel.

de tijd – en nauwkeurigheidsresultaten bij het gebruik van Megablast als classifier werden verkregen uit de loggegevens die door PhymmBL werden geproduceerd, omdat PhymmBL Megablast gebruikt voor zijn uitlijningsstap. Bij het toewijzen van een taxonomisch label aan een read met Megablast, gebruikten we het taxon geassocieerd met de eerste gerapporteerde uitlijning. Megablast werd uitgevoerd met standaard opties.

snelheid werd geëvalueerd met behulp van de single-threaded werking van elk programma (behalve voor NBC). PhymmBL werd gewijzigd zodat de oproep naar het blastn-programma gebruikt een thread in plaats van twee. NBC werd in werking gesteld met 36 gelijktijdige processen die op disjuncte reeksen genomen in zijn genomic bibliotheek werken, en de totale tijd voor de classifier werd bepaald door de decompressie en het noteren tijden voor elk genoom op te tellen. Wandkloktijden werden geregistreerd voor alle classifiers. Bij het vergelijken van Kraken met de andere classifiers, gebruikten we BLAST + 2.2.27, PhymmBL 4.0, NBC 1.1 en MetaPhlAn 1.7.6. Classifiers werden allemaal uitgevoerd op dezelfde computer, met 48 AMD Opteron 6172 2.1 GHz CPU ‘ s en 252 GB RAM, met Red Hat Enterprise Linux 5. De datasets die voor snelheidsevaluatie werden gebruikt, hadden elk 10.000 reads voor alle andere programma ‘ s dan Kraken (en zijn varianten) en MetaPhlAn, die 10.000.000 gelezen datasets gebruikten. Hogere leesnummers werden gebruikt met deze snellere programma’ s om het effect van de initiële en definitieve bewerkingen die plaatsvinden tijdens de uitvoering van de programma ‘ s te minimaliseren.

hoewel Kraken het enige programma is dat expliciet operaties uitvoert om ervoor te zorgen dat de gegevens in het fysieke geheugen zijn voordat ze worden geclassificeerd, wilden we er zeker van zijn dat alle programma ‘ s op dezelfde manier werden geëvalueerd. Bij het evalueren van snelheid, voor elk programma, lezen we alle databasebestanden (bijv. IMM-bestanden en BLAST-databases voor PhymmBL, K-Mer-frequentielijsten voor NBC en de Bowtie-index voor MetaPhlAn) drie keer in het geheugen voordat u het programma uitvoert, om de inhoud van de database in de cache van het besturingssysteem te plaatsen (die is opgeslagen in fysiek geheugen).

gereduceerde databasegrootte

om de 4-gb-database voor onze minikraken-resultaten te genereren, verwijderden we de eerste 18 van elk blok van 19 records in de standaard Kraken-database. Een krimpende factor van 19 werd gekozen omdat het de kleinste integer factor was die de grootte zou verminderen tot minder dan 4 GB, een grootte die gemakkelijk in het geheugen van veel gangbare personal computers past. Voor gebruikers die meer RAM beschikbaar hebben, Kraken kan een kleinere krimpfactor worden gebruikt, wat een verhoogde gevoeligheid zal geven.

gebruik van concept genomen

bij het construeren van de Kraken-GB database merkten we dat er verschillende contigs waren met bekende adaptersequenties aan de uiteinden. In latere tests, vonden we ook dat sommige sequenties in steekproeven met grote hoeveelheden menselijke sequentie consequent verkeerd werden geclassificeerd door deze database, wat ons leidde tot de conclusie dat besmetting waarschijnlijk aanwezig was in de concept genomen. In een poging om deze besmetting tegen te gaan, verwijderden we uit de database die k-mers van bekende adapter sequenties, evenals de eerste en laatste 20 k-mers van elk van de draft contigs. Hoewel dit de classificatie verbeterde, elimineerde het het probleem van de verkeerde indeling niet. Om deze reden, geloven wij dat als concept genomen in een Kraken database worden gebruikt, zeer strenge maatregelen zouden moeten worden gebruikt om de opeenvolgingen van de verontreiniging uit de genomic bibliotheek te verwijderen.

Clade-uitsluitingsexperimenten

bij het opnieuw analyseren van de Simba-5-dataset voor onze clade-uitsluitingsexperimenten werden sommige reads niet gebruikt voor bepaalde paren gemeten en uitgesloten rangen. Als de oorsprong van een read een taxonomische vermelding op een van de gemeten of uitgesloten rangen ontbrak, werd het niet gebruikt voor dat specifieke experiment.

bovendien werd een read niet gebruikt in een experiment, tenzij ten minste twee andere taxa in onze database (afgezien van de Uitgesloten clade) op de Uitgesloten rang deelden met het taxon van de oorsprong op de gemeten rang. Bijvoorbeeld, een read from genus G zou niet worden gebruikt in een experiment meten nauwkeurigheid op de klasse rang en exclusief de genus rang tenzij G ‘S thuis klasse had ten minste twee andere geslachten met genomen in Kraken’ s genomische bibliotheek. Zonder deze filterstap, zou een genus uitgesloten zijn toen het het enige genus in zijn klasse was, Kraken zou onmogelijk de juiste klasse kunnen benoemen, aangezien alle vermeldingen in de database van die klasse ook zouden worden uitgesloten. Dit is dezelfde aanpak genomen in soortgelijke experimenten die werden gebruikt om PhymmBL te evalueren .

Human microbiome data classification

we classificeerden de Human Microbiome Project data met behulp van een Kraken database gemaakt van complete RefSeq bacteriële, archaeale en virale genomen, samen met het grch37 menselijke genoom. We haalden de sequenties van drie toetredingen (SRS019120, SRS014468 en SRS015055) uit het NCBI Sequence Read Archive, en elke toetreding had twee runs ingediend. Alle reads werden bijgesneden om lage kwaliteit basissen en adapter sequenties te verwijderen. Krona werd gebruikt om alle taxonomische verdeelpunten te genereren.

omdat de sequenties allemaal gepaarde reads waren, voegden we de reads samen door de mates samen te voegen met een sequentie van ‘NNNNN’ ertussen. Kraken negeert k-mers met dubbelzinnige nucleotiden, dus de k-mers die deze ‘N’ karakters overspannen hebben geen invloed op de classificatie. Door deze operatie kon Kraken een paar reads als een enkele eenheid classificeren in plaats van de mates afzonderlijk te classificeren.

beschikbaarheid van Software en gegevens