Kraken: ultrafast metagenomic sequence classification using exact alignments

Sequence classification algorithm

för att klassificera en DNA-sekvens S samlar vi alla k-mers inom den sekvensen i en uppsättning, betecknad som K(S). Vi kartlägger sedan varje k – mer i K(S), med hjälp av algoritmen som beskrivs nedan, till LCA taxon av alla genom som innehåller den k-meren. Dessa LCA taxa och deras förfäder i taxonomiträdet bildar vad vi kallar klassificeringsträdet, ett beskuren underträd som används för att klassificera S. varje nod i klassificeringsträdet viktas med antalet k-mers i K(S) som mappas till taxonet associerat med den noden. Därefter görs varje rot-till-blad (RTL) – bana i klassificeringsträdet genom att beräkna summan av alla nodvikter längs banan. Den maximala poäng RTL-sökvägen i klassificeringsträdet är klassificeringsvägen och S tilldelas etiketten som motsvarar dess blad (om det finns flera maximalt poängvägar väljs LCA för alla dessa banors blad). Denna algoritm, illustrerad i Figur 1, tillåter Kraken att betrakta varje k-mer inom en sekvens som ett separat bevis och sedan försöka lösa eventuella motstridiga bevis om det behövs. Observera att för ett lämpligt val av k kommer de flesta k-mers att kartlägga unikt till en enda art, vilket förenklar klassificeringsprocessen. Sekvenser för vilka ingen av k-mers i K(S) finns i något genom lämnas oklassificerade av denna algoritm.

användningen av RTL-sökvägspoäng i klassificeringsträdet är nödvändigt mot bakgrund av de oundvikliga skillnaderna mellan sekvenserna som ska klassificeras och sekvenserna som finns i alla genombibliotek. Sådana skillnader kan, även för stora värden på k, resultera i en k-mer som finns i biblioteket men associerad med en art långt borta från den sanna källarten. Genom att poängsätta de olika RTL-banorna i klassificeringsträdet kan vi kompensera för dessa skillnader och korrekt klassificera sekvenser även när en liten minoritet av k-mers i en sekvens indikerar att sekvensen ska tilldelas en felaktig taxonomisk etikett.

Databasskapande

effektiv implementering av Krakens klassificeringsalgoritm kräver att kartläggningen av k-mers till taxa utförs genom att fråga en förberäknad databas. Kraken skapar denna databas genom en flerstegsprocess, som börjar med valet av ett bibliotek med genomiska sekvenser. Kraken innehåller ett standardbibliotek, baserat på slutförda mikrobiella genom i National Center for Biotechnology Information ’ s (NCBI) RefSeq-databas, men biblioteket kan anpassas efter behov av enskilda användare .

när biblioteket har valts använder vi Manet multithreaded k-Mer-räknaren för att skapa en databas som innehåller varje distinkt 31-mer i biblioteket. När databasen är klar, de 4-byte utrymmen maneter används för att lagra k-Mer räknas i databasfilen används istället av Kraken för att lagra taxonomiska ID-nummer för k-mers’ LCA-värden. Efter att databasen har skapats av maneter bearbetas de genomiska sekvenserna i biblioteket en i taget. För varje sekvens används taxonet som är associerat med det för att ställa in de lagrade LCA-värdena för alla k-mers i sekvensen. När sekvenser bearbetas, om en k-Mer från en sekvens har haft sitt LCA-värde tidigare inställt, beräknas LCA för det lagrade värdet och den aktuella sekvensens taxon och att LCA lagras för k-meren. Taxon information erhålls från NCBI taxonomy database.

databasstruktur och sökalgoritm

eftersom Kraken mycket ofta använder en k-Mer som en databasfråga omedelbart efter att ha frågat en intilliggande k-mer, och eftersom intilliggande k-mers delar en betydande mängd sekvens, använder vi minimizer-konceptet för att gruppera liknande k-mers tillsammans. För att förklara vår tillämpning av detta koncept definierar vi här den kanoniska representationen av en DNA-sekvens S som den lexikografiskt mindre av S och det omvända komplementet av S. För att bestämma en k-Mers minimizer av Längd M, betraktar vi den kanoniska representationen av alla M-mers i k-meren och väljer den lexikografiskt minsta av dessa M-mers som k-Mers minimizer. I praktiken har intilliggande k-mers ofta samma minimizer.

i Krakens databas lagras alla k-mers med samma minimizer i följd och sorteras i lexikografisk ordning av deras kanoniska representationer. En fråga för en k-mer R kan sedan behandlas genom att leta upp i ett index positionerna i databasen där k-mers med R: S minimizer skulle lagras och sedan utföra en binär sökning inom den regionen (Figur 5). Eftersom intilliggande k-mers ofta har samma minimizer är sökområdet ofta detsamma mellan två på varandra följande frågor, och sökningen i den första frågan kan ofta ta med data till CPU-cachen som kommer att användas i den andra frågan. Genom att tillåta minnesåtkomst i efterföljande frågor för att komma åt data i CPU-cachen istället för RAM, gör denna strategi efterföljande frågor mycket snabbare än de annars skulle vara. Indexet som innehåller förskjutningarna för varje grupp av k-mers i databasen kräver 8 4m-byte av 4 miljoner. Som standard använder Kraken 15-BP minimizers, men användaren kan ändra detta värde; till exempel, för att skapa MiniKraken, använde vi 13-bp minimizers för att säkerställa att den totala databasstorleken stannade under 4 GB.

Kraken databasstruktur. Varje k-mer som ska frågas mot databasen har en specifik delsträng som är dess minimizer. För att söka efter en k-mer i databasen undersöks positionerna i databasen som innehåller k-mers med samma minimizer. Dessa positioner hittas snabbt genom att undersöka minimizer offset array för startpositionerna för poster med k-mer ’ s minimizer (orange) och nästa möjliga minimizer (blå). Inom ett intervall av poster som är associerade med en given minimizer sorteras poster efter lexikografisk ordning av deras k-mers, vilket gör att en fråga kan slutföras med hjälp av en binär sökning över detta intervall.

vid implementeringen av Kraken gjorde vi ytterligare optimeringar av strukturen och sökalgoritmen som beskrivs ovan. Först, som noterat av Roberts et al. , en enkel lexikografisk beställning av M – mers kan resultera i en skev fördelning av minimisatorer som överrepresenterar M-Mers med låg komplexitet. I Kraken skulle en sådan bias skapa många stora sökområden, vilket skulle kräva mer tid att söka. För att skapa en jämnare fördelning av minimizers (och därmed påskynda sökningar) använder vi den exklusiva-eller (XOR) operationen för att växla hälften av bitarna i varje M-Mers kanoniska representation innan man jämför M-mers med varandra med hjälp av lexikografisk beställning. Denna XOR-operation förvränger effektivt standardbeställningen och förhindrar den stora förspänningen mot minimeringsmedel med låg komplexitet.

vi utnyttjar också det faktum att sökområdet ofta är detsamma mellan frågor för att göra Krakens frågor snabbare. I stället för att beräkna minimizer varje gång vi utför en fråga söker vi först i föregående intervall. Om vår frågade k-mer finns i detta intervall kan frågan returneras omedelbart. Om k-meren inte hittas, beräknas minimizer; om k-merens minimizer är densamma som den senast frågade k-mer, misslyckas frågan, eftersom minimizerens sökutrymme har visat sig inte ha k-meren. Endast om minimizer har ändrats måste Kraken justera sökområdet och söka igen efter k-meren.

konstruera simulerade metagenomer

hiseq-och MiSeq-metagenomerna byggdes med 20 uppsättningar av bakteriell helgenom Hagelgevär läser. Dessa läsningar hittades antingen som en del av gage-B-projektet eller i NCBI Sequence Read Archive. Varje metagenom innehåller sekvenser från tio genom (ytterligare fil 1: Tabell S1). För både 10 000 och 10 miljoner lästa prover av var och en av dessa metagenomer valdes 10% av deras sekvenser från var och en av de tio komponentgenomdatauppsättningarna (dvs. varje Genom hade lika sekvensöverflöd). Alla sekvenser trimmades för att ta bort baser och adaptersekvenser av låg kvalitet.

sammansättningen av dessa två metagenomer utgör vissa utmaningar för våra klassificerare. Till exempel kan Pelosinus fermentans, som finns i vårt HiSeq-metagenom, inte identifieras korrekt på genusnivå av Kraken (eller någon av de andra tidigare beskrivna klassificerarna), eftersom det inte finns några Pelosinus-genom i RefSeq complete genomes database; det finns dock sju sådana Genom i Kraken-GB: s databas, inklusive sex stammar av P. fermentans. På samma sätt klassificeras Proteus vulgaris i vårt MiSeq-metagenom ofta felaktigt på genusnivå eftersom det enda Proteusgenomet i Krakens databas är ett enda Proteus mirabilis-genom. Ytterligare fem Proteus-genom finns i Kraken-GB: s databas, vilket gör att Kraken-GB kan klassificera läser bättre från det släktet. Dessutom innehåller MiSeq-metagenomet fem Genom från familjen Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus och Salmonella). Den höga sekvenslikheten mellan släktena i denna familj kan göra det svårt att skilja mellan släkten för alla klassificerare.

simBA – 5-metagenomet skapades genom att simulera läsningar från uppsättningen kompletta bakteriella och archaeala genom i RefSeq. Replikoner från dessa genom användes om de var associerade med en taxon som hade en post associerad med släktet rank, vilket resulterade i en uppsättning replikoner från 607 släkten. Vi använde sedan Mason read-simulatorn med sin Illumina-modell för att producera 10 miljoner 100-bp-läsningar från dessa genom. Först skapade vi simulerade genom för varje art, med hjälp av en SNP-hastighet på 0,1% och en indelhastighet på 0,1% (båda standardparametrarna), från vilka vi genererade läsningarna. För de simulerade läsningarna multiplicerade vi standardmatchnings-och indelhastigheterna med fem, vilket resulterade i en genomsnittlig felmatchningsgrad på 2% (från 1% i början av läsningar till 6% i ändarna) och en indelhastighet på 1% (0,5% införingssannolikhet och 0,5% borttagningssannolikhet). För Simba-5 metagenom genererades 10 000 läsuppsättningar från ett slumpmässigt urval av 10 miljoner läsuppsättningar.

utvärdering av noggrannhet och hastighet

vi valde att mäta noggrannhet främst på genusnivå, vilket var den lägsta nivån för vilken vi lätt kunde bestämma taxonomiinformationen för PhymmBL och NBC: s förutsägelser på ett automatiserat sätt. (Detta beror på hur PhymmBL och NBC rapporterar sina resultat). Eftersom vissa Genom inte har taxonomiska poster i alla sju LED (Art, släkt, familj, ordning, klass, fylum och rike) definierade vi genusnivåkänslighet som A/B, där A är antalet läsningar med ett tilldelat släkt som korrekt klassificerades vid den rangordningen, och B är det totala antalet läsningar med något tilldelat släkt. Vi definierade känslighet på samma sätt för andra taxonomiska LED.

eftersom Kraken kan klassificera en läsning på nivåer över arten, kräver mätning av dess precision att vi definierar effekten på precisionen att tilldela rätt släkt (till exempel) medan vi inte tilldelar en art alls. Av denna anledning definierade vi rangnivåprecision som C / (D + E), där C är antalet läsningar märkta vid eller under rätt taxon vid den uppmätta rangen, D är antalet läsningar märkta vid eller under den uppmätta rangen, och E är antalet läsningar felaktigt märkta ovanför den uppmätta rangen. Till exempel, med tanke på en läs R som bör märkas som Escherichia coli, skulle en märkning av R som E. coli, E. fergusonii eller Escherichia förbättra precisionen på genusnivå. En etikett av Enterobacteriaceae (korrekt familj) eller Proteobakterier (korrekt phylum) skulle inte ha någon effekt på genusnivå precision. En etikett för R av Bacillus (felaktigt släkt) eller Firmicutes (felaktigt fylum) skulle minska släktnivåprecisionen.

när vi utvärderade Phymmbls noggrannhet , efter Utvecklarens råd, valde vi en genusförtroende tröskel för våra jämförelser. Vi valde 3,333 läser från simulerad medium komplexitet (simMC) datamängd, som täcker 31 olika släkten. För att simulera korta läsningar från Sanger-sekvensdata i simMC-uppsättningen valde vi de sista 100 bp från var och en av läsningarna. Vi sprang sedan PhymmBL mot de 100-bp-läsningar och utvärderade genusnivåkänsligheten och precisionen hos Phymmbls förutsägelser med genusförtroende tröskelvärden från 0 till 1, i steg om 0.05. Vi fann att en tröskel på 0.65 gav den högsta F-poängen (det harmoniska medelvärdet av känslighet och precision), med 0.60 och 0.70 som också hade F-poäng inom 0.5 procentenheter av det maximala (tilläggsfil 1: Tabell S2). Vi använde därför 0.65 genus confidence threshold i våra jämförelser. Även om valet av ett tröskelvärde beror på användarens individuella behov, och det är till viss del godtyckligt, ger en tröskel som väljs på detta sätt en mer korrekt jämförelse med en selektiv klassificerare som Kraken än ingen tröskel alls.

tid och noggrannhet resultat vid användning Megablast som en klassificerare erhölls från loggdata som produceras av PhymmBL, som PhymmBL använder Megablast för dess inriktningssteg. När du tilldelar en taxonomisk etikett till en läsning med Megablast, vi använde taxon i samband med den första rapporterade anpassningen. Megablast kördes med standardalternativ.

hastighet utvärderades med hjälp av den enkelgängade operationen för varje program (förutom NBC). PhymmBL ändrades så att dess samtal till blastn-programmet använde en tråd istället för två. NBC kördes med 36 samtidiga processer som arbetar på separata uppsättningar genom i sitt genomiska bibliotek, och den totala tiden för klassificeraren bestämdes genom att summera dekompressions-och poängstiderna för varje genom. Väggklocktider registrerades för alla klassificerare. Vid jämförelse av Kraken med de andra klassificerarna använde vi BLAST+ 2.2.27, PhymmBL 4.0, NBC 1.1 och MetaPhlAn 1.7.6. Klassificerare kördes alla på samma dator, med 48 AMD Opteron 6172 2.1 GHz-processorer och 252 GB RAM, kör Red Hat Enterprise Linux 5. Datamängderna som användes för hastighetsutvärdering hade 10 000 läsningar vardera för alla andra program än Kraken (och dess varianter) och MetaPhlAn, som använde 10 000 000 lästa datamängder. Högre läsnummer användes med dessa snabbare program för att minimera effekten av de initiala och slutliga operationerna som äger rum under programmens körning.

även om Kraken är det enda av de program vi undersökte som uttryckligen utför operationer för att säkerställa att dess data finns i fysiskt minne före klassificering, ville vi vara säkra på att alla program utvärderades på ett liknande sätt. Vid utvärdering av hastighet, för varje program, läser vi alla databasfiler (t. ex. IMM-filer och BLAST-databaser för PhymmBL, k-Mer-frekvenslistor för NBC och Bowtie index för MetaPhlAn) i minnet tre gånger innan programmet körs för att placera databasinnehållet i operativsystemets cache (som lagras i fysiskt minne).

reducerade databasstorlekar

för att generera 4-GB-databasen för våra MiniKraken-resultat tog vi bort de första 18 i varje block med 19 poster i standard Kraken-databasen. En krympande faktor på 19 valdes eftersom det var den minsta heltalsfaktorn som skulle minska storleken till mindre än 4 GB, en storlek som lätt kan passa in i minnet på många vanliga persondatorer. För användare som har mer RAM tillgängligt tillåter Kraken en mindre krympningsfaktor att användas, vilket ger ökad känslighet.

användning av utkast till genom

vid konstruktion av Kraken-GB-databasen märkte vi att det fanns flera contigs med kända adaptersekvenser i ändarna. I efterföljande tester fann vi också att vissa sekvenser i prover med stora mängder mänsklig sekvens konsekvent felklassificerades av denna databas, vilket ledde oss att dra slutsatsen att kontaminering sannolikt var närvarande i utkastet till genom. I ett försök att motverka denna förorening, vi bort från databasen dessa k-mers från kända adaptersekvenser, liksom de första och sista 20 k-mers från var och en av utkastet contigs. Även om detta förbättrade klassificeringen eliminerade det inte felklassificeringsproblemet. Av denna anledning tror vi att om utkast till genom används i en Kraken-databas, bör mycket stränga åtgärder användas för att ta bort föroreningssekvenser från det genomiska biblioteket.

Clade-uteslutningsexperiment

vid omanalys av simBA-5-datamängden för våra clade-uteslutningsexperiment användes vissa läsningar inte för vissa par av uppmätta och uteslutna LED. Om en läses ursprung saknade en taxonomisk post i någon av de uppmätta eller uteslutna leden, användes den inte för just det experimentet.

dessutom användes inte en läsning i ett experiment om inte minst två andra taxa representerade i vår databas (bortsett från den uteslutna kladen) vid den uteslutna rangen delade ursprungsbladets taxon vid den uppmätta rankningen. Till exempel skulle en läsning från genus G inte användas i ett experiment som mäter noggrannhet vid klassrankningen och utesluter genusrankningen om inte G: s hemklass hade minst två andra släkten med genom i Krakens genomiska bibliotek. Utan detta filtreringssteg, var ett släkte uteslutet när det var det enda släktet i sin klass, kunde Kraken omöjligt namnge rätt klass, eftersom alla poster i databasen från den klassen också skulle uteslutas. Detta är samma tillvägagångssätt som tagits i liknande experiment som användes för att utvärdera PhymmBL .

Human microbiome Data classification

vi klassificerade Human Microbiome Project-data med hjälp av en Kraken-databas gjord av kompletta RefSeq-bakteriella, archaeala och virala genom, tillsammans med grch37-mänskliga genomet. Vi hämtade sekvenserna av tre anslutningar (SRS019120, SRS014468 och SRS015055) från NCBI Sequence Read Archive, och varje anslutning hade två inlämnade körningar. Alla läsningar trimmades för att ta bort baser och adaptersekvenser av låg kvalitet. Krona användes för att generera alla taxonomiska fördelningsområden.

eftersom sekvenserna var alla parade läsningar, gick vi ihop läsarna genom att sammanfoga kompisarna med en sekvens av ’NNNNN’ mellan dem. Kraken ignorerar k-mers med tvetydiga nukleotider, så k-mers som spänner över dessa ’N’ – tecken påverkar inte klassificeringen. Denna operation gjorde det möjligt för Kraken att klassificera ett par läsningar som en enda enhet snarare än att behöva klassificera kompisarna separat.

programvara och datatillgänglighet