Kraken: ultrafast metagenomisk sekvensklassificering ved hjælp af nøjagtige justeringer

Sekvensklassificeringsalgoritme

for at klassificere en DNA-sekvens S samler vi alle k-mers inden for denne sekvens i et sæt, betegnet som K(s). Vi kortlægger derefter hver k-mer i K(S) ved hjælp af algoritmen beskrevet nedenfor til LCA-taksonen for alle genomer, der indeholder den k-mer. Disse LCA-takser og deres forfædre i taksonomitræet danner det, vi kalder klassificeringstræet, et beskåret undertræ, der bruges til at klassificere S. hver knude i klassificeringstræet vægtes med antallet af k-mers i K(S), der er kortlagt til taksonet, der er knyttet til den knude. Derefter scorer hver rod-til-blad (RTL) sti i klassificeringstræet ved at beregne summen af alle knudevægte langs stien. Den maksimale scorings-RTL-sti i klassificeringstræet er klassificeringsstien, og S tildeles etiketten svarende til dens blad (hvis der er flere maksimalt scoringsstier, vælges LCA for alle disse Stiers blade). Denne algoritme, illustreret i Figur 1, giver Kraken mulighed for at betragte hver k-mer inden for en sekvens som et separat bevis og derefter forsøge at løse eventuelle modstridende beviser, hvis det er nødvendigt. Bemærk, at for et passende valg af k, de fleste k-mers vil kortlægge entydigt til en enkelt art, i høj grad forenkle klassificeringsprocessen. Sekvenser, for hvilke ingen af k-merne i K(S) findes i noget genom, efterlades uklassificeret af denne algoritme.

brugen af RTL-stiscoring i klassificeringstræet er nødvendig i lyset af de uundgåelige forskelle mellem sekvenserne, der skal klassificeres, og sekvenserne, der findes i ethvert bibliotek med genomer. Sådanne forskelle kan, selv for store værdier af k, resultere i en k-mer, der er til stede i biblioteket, men forbundet med en art langt væk fra den sande kildeart. Ved at score de forskellige RTL-stier i klassificeringstræet kan vi kompensere for disse forskelle og korrekt klassificere sekvenser, selv når et lille mindretal af k-mers i en sekvens indikerer, at sekvensen skal tildeles en forkert taksonomisk etiket.

oprettelse af Database

effektiv implementering af Krakens klassificeringsalgoritme kræver, at kortlægningen af k-mers til taksa udføres ved at forespørge en forudberegnet database. Kraken opretter denne database gennem en flertrinsproces, der begynder med udvælgelsen af et bibliotek med genomiske sekvenser. Kraken inkluderer et standardbibliotek baseret på udfyldte mikrobielle genomer i National Center for Biotechnology Information ‘ s (NCBI) referencedatabase, men biblioteket kan tilpasses efter behov af individuelle brugere .

når biblioteket er valgt, bruger vi Jellyfish multithreaded k-mer-tælleren til at oprette en database, der indeholder hver særskilt 31-mer i biblioteket. Når databasen er færdig, bruges de 4-byte mellemrum vandmænd, der bruges til at gemme k-mer-tællingerne i databasefilen, i stedet af Kraken til at gemme de taksonomiske ID-numre for k-mers’ LCA-værdier. Når databasen er oprettet af vandmænd, behandles de genomiske sekvenser i biblioteket en ad gangen. For hver sekvens bruges den tilknyttede takson til at indstille de lagrede LCA-værdier for alle k-mers i sekvensen. Når sekvenser behandles, hvis en k-mer fra en sekvens har haft sin LCA-værdi tidligere indstillet, beregnes LCA for den lagrede værdi og den aktuelle sekvenss takson, og at LCA gemmes for k-mer. Oplysninger om taksonomi fås fra NCBI-taksonomidatabasen.

databasestruktur og søgealgoritme

fordi Kraken meget ofte bruger en k-mer som en databaseforespørgsel umiddelbart efter at have spurgt en tilstødende k-mer, og fordi tilstødende k-mers deler en betydelig mængde sekvens, bruger vi minimeringskonceptet til at gruppere lignende k-mers sammen. For at forklare vores anvendelse af dette koncept, vi definerer her den kanoniske repræsentation af en DNA-sekvens s som den leksikografisk mindre af S og det omvendte komplement af S. For at bestemme en k-Mers minimering af længde M, betragter vi den kanoniske repræsentation af alle M-mers i k-mer og vælger den leksikografisk mindste af disse M-mers som k-Mers minimering. I praksis vil tilstødende k-mers ofte have den samme minimering.

i Krakens database gemmes alle k-mers med samme minimering fortløbende og sorteres i leksikografisk rækkefølge efter deres kanoniske repræsentationer. En forespørgsel om en k-mer R kan derefter behandles ved at kigge op i et indeks positionerne i databasen, hvor k-mers med R ‘ S minimering ville blive gemt, og derefter udføre en binær søgning inden for denne region (figur 5). Da tilstødende k-mers ofte har den samme minimering, er søgeområdet ofte det samme mellem to på hinanden følgende forespørgsler, og søgningen i den første forespørgsel kan ofte bringe data ind i CPU-cachen, der vil blive brugt i den anden forespørgsel. Ved at tillade hukommelsesadgang i efterfølgende forespørgsler at få adgang til data i CPU-cachen i stedet for RAM, gør denne strategi efterfølgende forespørgsler meget hurtigere, end de ellers ville være. Indekset, der indeholder forskydningerne for hver gruppe af k-mers i databasen, kræver 8 milliliter 4m bytes. Som standard bruger Kraken 15-bp minimisatorer, men brugeren kan ændre denne værdi; for eksempel ved oprettelse af MiniKraken brugte vi 13-bp minimisatorer for at sikre, at den samlede databasestørrelse forblev under 4 GB.

Kraken database struktur. Hver k-mer, der skal forespørges mod databasen, har en specifik substring, der er dens minimering. For at søge efter en k-mer i databasen undersøges positionerne i databasen, der indeholder k-mers med samme minimering. Disse positioner findes hurtigt ved at undersøge minimeringsforskydningsarrayet for startpositionerne for poster med k-Mers minimering (orange) og den næste mulige minimering (blå). Inden for en række poster, der er knyttet til en given minimering, sorteres poster efter leksikografisk rækkefølge af deres k-mers, så en forespørgsel kan udfyldes ved hjælp af en binær søgning over dette interval.

ved implementeringen af Kraken foretog vi yderligere optimeringer af strukturen og søgealgoritmen beskrevet ovenfor. Først, som bemærket af Roberts et al. , en simpel leksikografisk rækkefølge af M-mers kan resultere i en skæv fordeling af minimeringsmidler, der overrepræsenterer M-mers med lav kompleksitet. I Kraken ville en sådan bias skabe mange store søgeområder, hvilket ville kræve mere tid til at søge. For at skabe en mere jævn fordeling af minimisatorer (og dermed fremskynde søgninger) bruger vi den eksklusive-eller-handling til at skifte halvdelen af bitene i hver M-Mers kanoniske repræsentation, inden vi sammenligner M-mers med hinanden ved hjælp af leksikografisk rækkefølge. Denne operation krypterer effektivt standardbestillingen og forhindrer den store bias mod minimering af lav kompleksitet.

vi drager også fordel af det faktum, at søgeområdet ofte er det samme mellem forespørgsler for at gøre Krakens forespørgsler hurtigere. I stedet for at beregne minimeringen hver gang vi udfører en forespørgsel, søger vi først i det forrige interval. Hvis vores forespurgte k-mer findes i dette interval, kan forespørgslen straks vende tilbage. Hvis k-mer ikke findes, beregnes minimisatoren; hvis k-mer ‘s minimisator er den samme som de sidste forespurgte k-mer’ er, mislykkes forespørgslen, da minimeringens søgeplads har vist sig ikke at have k-mer. Kun hvis minimeringen er ændret, skal Kraken justere søgeområdet og søge igen efter k-mer.

konstruktion af simulerede metagenomer

Hisek og Misek metagenomer blev bygget ved hjælp af 20 sæt bakterielle helgenom haglgevær læser. Disse læsninger blev fundet enten som en del af Gage-B-projektet eller i NCBI-Sekvenslæsearkivet. Hvert metagenom indeholder sekvenser fra ti genomer (yderligere fil 1: tabel S1). For både 10.000 og 10 millioner læste prøver af hvert af disse metagenomer blev 10% af deres sekvenser valgt fra hvert af de ti komponentgenomdatasæt (dvs.hvert genom havde samme sekvensmængde). Alle sekvenser blev trimmet for at fjerne baser og adaptersekvenser af lav kvalitet.

sammensætningen af disse to metagenomer udgør visse udfordringer for vores klassifikatorer. For eksempel kan Pelosinus fermentans, der findes i vores hisek metagenom, ikke identificeres korrekt på slægtsniveau af Kraken (eller nogen af de andre tidligere beskrevne klassifikatorer), fordi der ikke er nogen Pelosinus genomer i den komplette genomdatabase; der er dog syv sådanne genomer i Kraken-GB ‘ s database, herunder seks stammer af P. fermentans. Tilsvarende klassificeres Proteus vulgaris i vores Misek metagenom ofte forkert på slægtsniveau, fordi det eneste Proteusgenom i Krakens database er et enkelt Proteus mirabilis genom. Fem flere Proteus genomer er til stede i Kraken-GBS database, så Kraken-GB kan klassificere læser bedre fra den slægt. Derudover indeholder Metagenomet fem genomer fra Enterobacteriaceae-familien (Citrobacter, Enterobacter, Klebsiella, Proteus og Salmonella). Den høje sekvenslighed mellem slægterne i denne familie kan gøre det vanskeligt at skelne mellem slægter for enhver klassifikator.

Simba-5-metagenomet blev skabt ved at simulere læsninger fra sættet af komplette bakterielle og arkaeale genomer i Refseks. Replikoner fra disse genomer blev brugt, hvis de var forbundet med en takson, der havde en post forbundet med slægtsrangeringen, hvilket resulterede i et sæt replikoner fra 607 slægter. Vi brugte derefter Mason read simulator med sin Illumina model til at producere 10 millioner 100-bp læser fra disse genomer. Først oprettede vi simulerede genomer for hver art ved hjælp af en SNP-hastighed på 0,1% og en indel-hastighed på 0,1% (begge standardparametre), hvorfra vi genererede læsningerne. For de simulerede aflæsninger multiplicerede vi standard mismatch-og indel-satserne med fem, hvilket resulterede i en gennemsnitlig mismatch-sats på 2% (fra 1% i begyndelsen af læsninger til 6% i enderne) og en indel-sats på 1% (0,5% indsættelsessandsynlighed og 0,5% sletningssandsynlighed). For simBA – 5 metagenome blev 10.000 læsesæt genereret fra en tilfældig prøve af 10 millioner læsesæt.

evaluering af nøjagtighed og hastighed

vi valgte at måle nøjagtighed primært på slægtsniveau, hvilket var det laveste niveau, som vi let kunne bestemme taksonomioplysningerne for PhymmBL og NBC ‘ s forudsigelser på en automatiseret måde. (Dette skyldes den måde, hvorpå PhymmBL og NBC rapporterer deres resultater). Fordi nogle genomer ikke har taksonomiske poster i alle syv rækker (art, slægt, familie, orden, klasse, phylum og Kongerige), definerede vi genus-niveau følsomhed som A/B, hvor A er antallet af læsninger med en tildelt slægt, der blev korrekt klassificeret på den rang, og B er det samlede antal læsninger med enhver tildelt Slægt. Vi definerede følsomhed på samme måde for andre taksonomiske rækker.

fordi Kraken kan klassificere en læsning på niveauer over arten, kræver måling af dens præcision, at vi definerer effekten på præcisionen ved at tildele den korrekte Slægt (for eksempel) uden at tildele en art overhovedet. Af denne grund definerede vi præcision på rangniveau som C/(D + E), hvor C er antallet af læsninger mærket ved eller under den korrekte takson ved den målte rang, D er antallet af læsninger mærket ved eller under den målte rang, og E er antallet af læsninger forkert mærket over den målte rang. For eksempel, givet en læse R, der skal mærkes som Escherichia coli, en mærkning af R som E. coli, E. fergusonii eller Escherichia ville forbedre genus-niveau præcision. En etiket af Enterobacteriaceae (korrekt familie) eller Proteobacteria (korrekt phylum) ville ikke have nogen effekt på genus-niveau præcision. En etiket til r af Bacillus (forkert slægt) eller Firmicutes (forkert phylum) ville reducere genus-niveau præcision.

ved evaluering af Phymmbls nøjagtighed valgte vi efter udviklernes råd en genus confidence threshold for vores sammenligninger. Vi valgte 3.333 læser fra simuleret medium kompleksitet (simMC) datasæt, der dækker 31 forskellige slægter. For at simulere korte læsninger fra Sanger-sekvensdataene i simMC-Sættet valgte vi de sidste 100 bp fra hver af læsningerne. Vi løb derefter PhymmBL mod disse 100-bp læser og vurderede genus-niveau følsomhed og præcision af Phymmbls forudsigelser med genus tillidstærskler fra 0 til 1 i trin på 0,05. Vi fandt ud af, at en tærskel på 0,65 gav den højeste F-score (det harmoniske gennemsnit af følsomhed og præcision), hvor 0,60 og 0,70 også havde F-score inden for 0.5 procentpoint af maksimum (yderligere fil 1: tabel S2). Vi brugte derfor 0.65 genus confidence threshold i vores sammenligninger. Selvom valget af en tærskel afhænger af en brugers individuelle behov, og det er også til en vis grad vilkårligt, giver en tærskel valgt på denne måde en mere korrekt sammenligning med en selektiv klassifikator som f.eks Kraken end slet ingen tærskel.

tids-og nøjagtighedsresultaterne ved brug af Megablast som klassifikator blev opnået fra logdata produceret af PhymmBL, da PhymmBL bruger Megablast til dets justeringstrin. Ved tildeling af en taksonomisk etiket til en læsning med Megablast, vi brugte taksonen i forbindelse med den første rapporterede justering. Megablast blev kørt med standardindstillinger.

hastighed blev evalueret ved hjælp af enkelttrådet drift af hvert program (undtagen NBC). PhymmBL blev ændret, så dets opkald til blastn-programmet brugte en tråd i stedet for to. NBC blev kørt med 36 samtidige processer, der opererer på uensartede sæt genomer i dets genomiske bibliotek, og den samlede tid for klassifikatoren blev bestemt ved at opsummere dekompressions-og scoringstiderne for hvert genom. Vægurtider blev registreret for alle klassifikatorer. Ved sammenligning af Kraken med de andre klassifikatorer brugte vi BLAST+ 2.2.27, PhymmBL 4.0, NBC 1.1 og Metaflan 1.7.6. Klassifikatorer blev alle kørt på den samme computer med 48 AMD Opteron 6172 2.1 CPU ‘ er og 252 GB RAM, der kørte Red Hat Enterprise 5. Datasættene, der blev brugt til hastighedsevaluering, havde 10.000 læsninger hver for alle andre programmer end Kraken (og dens varianter) og Metaflan, der brugte 10.000.000 læsedatasæt. Højere læsetal blev brugt med disse hurtigere programmer for at minimere effekten af de indledende og endelige operationer, der finder sted under programmernes udførelse.

selvom Kraken er det eneste af de programmer, vi undersøgte, der eksplicit udfører operationer for at sikre, at dens data er i fysisk hukommelse før klassificering, ønskede vi at være sikre på, at alle programmer blev evalueret på en lignende måde. Ved evaluering hastighed, for hvert program, vi læser alle databasefiler (f. eks. IMM filer og BLAST databaser for PhymmBL, k-mer frekvens lister for NBC og Butterfly indeks for Metaflan) i hukommelsen tre gange, før du kører programmet, for at placere databaseindholdet i operativsystemet cache (som er gemt i fysisk hukommelse).

reducerede databasestørrelser

for at generere 4 GB-databasen til vores minikraken-resultater fjernede vi de første 18 af hver blok med 19 poster i standard Kraken-databasen. En krympende faktor på 19 blev valgt, da det var den mindste heltalsfaktor, der ville reducere størrelsen til mindre end 4 GB, En størrelse, der let kan passe ind i hukommelsen på mange almindelige personlige computere. For brugere, der har mere RAM til rådighed, tillader Kraken, at der anvendes en mindre krympefaktor, hvilket vil give øget følsomhed.

brug af udkast til genomer

ved konstruktion af Kraken-GB-databasen bemærkede vi, at der var flere contigs med kendte adaptersekvenser i enderne. I efterfølgende tests fandt vi også, at nogle sekvenser i prøver med store mængder menneskelig sekvens konsekvent blev klassificeret forkert af denne database, hvilket førte til, at vi konkluderede, at forurening sandsynligvis var til stede i udkastet til genomer. I et forsøg på at modvirke denne forurening, vi fjernede disse k-mers fra databasen fra kendte adaptersekvenser, såvel som den første og sidste 20 k-mers fra hver af udkastet til contigs. Selvom dette forbedrede klassificeringen, eliminerede det ikke fejlklassificeringsproblemet. Af denne grund mener vi, at hvis udkast til genomer anvendes i en Kraken-database, bør der anvendes meget strenge foranstaltninger til at fjerne forurenende sekvenser fra det genomiske bibliotek.

Clade-udelukkelseseksperimenter

ved genanalyse af Simba-5-datasættet til vores clade-udelukkelseseksperimenter blev nogle læsninger ikke brugt til visse par målte og udelukkede rækker. Hvis en læses Oprindelse manglede en taksonomisk post i en af de målte eller udelukkede rækker, blev den ikke brugt til det pågældende eksperiment.

derudover blev en læsning ikke brugt i et eksperiment, medmindre mindst to andre takser repræsenteret i vores database (bortset fra den ekskluderede klade) på den ekskluderede rang delte oprindelseskladens takson på den målte rang. For eksempel, en læsning fra Slægt G ville ikke blive brugt i et eksperiment, der måler nøjagtighed ved klassens rang og ekskluderer slægtsrangeringen, medmindre GS hjemmeklasse havde mindst to andre slægter med genomer i Krakens genomiske bibliotek. Uden dette filtreringstrin, var en slægt udelukket, da den var den eneste slægt i sin klasse, Kraken kunne umuligt navngive den rigtige klasse, da alle poster i databasen fra denne klasse også ville blive udelukket. Dette er den samme tilgang taget i lignende eksperimenter, der blev brugt til at evaluere PhymmBL .

klassificering af humane mikrobiomdata

vi klassificerede de humane Mikrobiomprojektdata ved hjælp af en Kraken-database lavet af komplette bakterielle, arkaeale og virale genomer sammen med det menneskelige genom GRCh37. Vi hentede sekvenserne af tre tiltrædelser (SRS019120, SRS014468 og SRS015055) fra NCBI-Sekvenslæsearkivet, og hver tiltrædelse havde to kørsler indsendt. Alle læsninger blev trimmet for at fjerne baser og adaptersekvenser af lav kvalitet. Krona blev brugt til at generere alle taksonomiske distributionsplotter.

fordi sekvenserne alle var parrede læser, sluttede vi os sammen ved at sammenkæde kammeraterne med en sekvens af ‘NNNNN’ mellem dem. Kraken ignorerer k-mers med tvetydige nukleotider, så k-mers, der spænder over disse ‘N’ tegn, påvirker ikke klassificeringen. Denne operation gjorde det muligt for Kraken at klassificere et par læsninger som en enkelt enhed i stedet for at skulle klassificere kammeraterne separat.

tilgængelighed af programmel og data