Kraken: clasificarea secvenței metagenomice ultrarapide folosind aliniamentele exacte

algoritmul de clasificare a secvenței

pentru a clasifica o secvență ADN S, colectăm toți k-mers din acea secvență într-un set, notat ca K(S). Apoi mapăm fiecare k-mer în K(S), folosind algoritmul descris mai jos, la taxonul LCA al tuturor genomilor care conțin acel k-mer. Acești taxoni LCA și strămoșii lor din arborele taxonomiei formează ceea ce numim arborele de clasificare, un subarbor tăiat care este folosit pentru a clasifica S. fiecare nod din arborele de clasificare este ponderat cu numărul de K-mers în K(S) care mapat la taxonul asociat cu acel nod. Apoi, fiecare cale rădăcină-frunză (RTL) din arborele de clasificare este marcată prin calcularea sumei tuturor greutăților nodului de-a lungul căii. Calea maximă de notare RTL în arborele de clasificare este calea de clasificare, iar lui S i se atribuie eticheta corespunzătoare frunzei sale (dacă există mai multe căi de notare maximă, este selectată LCA a tuturor frunzelor acestor căi). Acest algoritm, ilustrat în Figura 1, îi permite lui Kraken să ia în considerare fiecare k-mer dintr-o secvență ca o dovadă separată și apoi să încerce să rezolve orice dovadă conflictuală, dacă este necesar. Rețineți că, pentru o alegere adecvată a k, majoritatea k-mers vor mapa în mod unic la o singură specie, simplificând foarte mult procesul de clasificare. Secvențele pentru care niciunul dintre k-mers în K(S) nu se găsește în niciun genom sunt lăsate neclasificate de acest algoritm.

utilizarea notării căii RTL în arborele de clasificare este necesară în lumina diferențelor inevitabile dintre secvențele care trebuie clasificate și secvențele prezente în orice bibliotecă de genomi. Astfel de diferențe pot, chiar și pentru valori mari de k, duce la un k-mer care este prezent în bibliotecă, dar asociat cu o specie departe de adevărata specie sursă. Prin notarea diferitelor căi RTL din arborele de clasificare, putem compensa aceste diferențe și clasifica corect secvențele chiar și atunci când o mică minoritate de k-mers într-o secvență indică faptul că secvenței ar trebui să i se atribuie o etichetă taxonomică incorectă.

Crearea bazei de date

implementarea eficientă a algoritmului de clasificare Kraken necesită ca maparea k-mers la taxoni să fie efectuată prin interogarea unei baze de date precalculate. Kraken creează această bază de date printr-un proces în mai multe etape, începând cu selectarea unei biblioteci de secvențe genomice. Kraken include o bibliotecă implicită, bazată pe genomii microbieni finalizați în baza de date RefSeq a Centrului Național pentru Informații Biotehnologice (NCBI), dar biblioteca poate fi personalizată după cum este necesar de către utilizatorii individuali .

odată ce biblioteca este aleasă, folosim contorul K-mer meduze multithreaded pentru a crea o bază de date care conține fiecare 31-mer distincte în bibliotecă. Odată ce baza de date este completă, meduzele de spații de 4 octeți utilizate pentru stocarea numărului k-mer în fișierul bazei de date sunt utilizate în schimb de Kraken pentru a stoca numerele de identificare taxonomice ale valorilor LCA ale k-mers. După ce baza de date a fost creată de meduze, secvențele genomice din bibliotecă sunt procesate una câte una. Pentru fiecare secvență, taxonul asociat cu acesta este utilizat pentru a seta valorile LCA stocate ale tuturor k-mers din secvență. Pe măsură ce secvențele sunt procesate, dacă un k-mer dintr-o secvență a avut valoarea LCA setată anterior, atunci LCA a valorii stocate și taxonul secvenței curente este calculat și că LCA este stocat pentru k-mer. Informațiile despre Taxon sunt obținute din Baza de date NCBI taxonomy.

structura bazei de date și algoritmul de căutare

deoarece Kraken folosește foarte frecvent un k-mer ca interogare a bazei de date imediat după interogarea unui k-mer adiacent și pentru că K-mers adiacent împărtășesc o cantitate substanțială de secvență, Folosim conceptul minimizer pentru a grupa K-mers similare împreună. Pentru a explica aplicarea acestui concept, definim aici reprezentarea canonică a unei secvențe ADN S ca lexicografic mai mic de S și complementul invers al lui S. Pentru a determina un minimizator de lungime M al lui k-mer, luăm în considerare reprezentarea canonică a tuturor M-mers în k-mer și selectăm lexicografic cel mai mic dintre acei m-mers ca minimizator al lui k-mer. În practică, k-mers adiacente vor avea adesea același minimizator.

în baza de date Kraken, toate K-mers cu același minimizator sunt stocate consecutiv și sunt sortate în ordinea lexicografică a reprezentărilor lor canonice. O interogare pentru un k-mer R poate fi apoi procesată căutând într-un index pozițiile din Baza de date în care ar fi stocate k-mers cu minimizatorul lui R și apoi efectuând o căutare binară în acea regiune (Figura 5). Deoarece k-mers-urile adiacente au adesea același minimizator, intervalul de căutare este adesea același între două interogări consecutive, iar căutarea din prima interogare poate aduce adesea date în memoria cache a procesorului care va fi utilizată în a doua interogare. Permițând accesul la memorie în interogările ulterioare pentru a accesa datele din memoria cache a procesorului în loc de RAM, această strategie face interogările ulterioare mult mai rapide decât ar fi altfel. Indexul care conține compensările fiecărui grup de K-mers din Baza de date necesită 8 octeți de 4 m XT. În mod implicit, Kraken folosește minimizatoare de 15 bp, dar utilizatorul poate modifica această valoare; de exemplu, în crearea MiniKraken, am folosit minimizatoare de 13 bp pentru a ne asigura că dimensiunea totală a bazei de date a rămas sub 4 GB.

structura bazei de date Kraken. Fiecare k-mer care trebuie interogat împotriva bazei de date are un subșir specific care este minimizatorul său. Pentru a căuta un k-mer în baza de date, sunt examinate pozițiile din Baza de date care conțin k-mers cu același minimizator. Aceste poziții sunt găsite rapid prin examinarea matricei de compensare a minimizatorului pentru pozițiile de pornire ale înregistrărilor cu minimizatorul k-mer (portocaliu) și următorul minimizator posibil (albastru). Într-un interval de înregistrări asociate cu un minimizator dat, înregistrările sunt sortate după ordonarea lexicografică A k-mers-urilor lor, permițând completarea unei interogări utilizând o căutare binară în acest interval.

în implementarea Kraken, am făcut optimizări suplimentare pentru structura și algoritmul de căutare descris mai sus. În primul rând, după cum a remarcat Roberts și colab. , o simplă ordonare lexicografică A M-mers poate duce la o distribuție înclinată a minimizatorilor care reprezintă în exces m-mers cu complexitate redusă. În Kraken, o astfel de părtinire ar crea multe Game mari de căutare, ceea ce ar necesita mai mult timp pentru a căuta. Pentru a crea o distribuție mai uniformă a minimizatorilor (și astfel a accelera căutările), folosim operația exclusivă-sau (XOR) pentru a comuta jumătate din biții reprezentării canonice a fiecărui M-mer înainte de a compara M-mers între ei folosind ordonarea lexicografică. Această operație XOR amestecă în mod eficient comanda standard și previne părtinirea mare față de minimizatoarele cu complexitate redusă.

de asemenea, profităm de faptul că intervalul de căutare este adesea același între interogări pentru a face interogările Kraken mai rapide. În loc să calculăm minimizatorul de fiecare dată când efectuăm o interogare, căutăm mai întâi intervalul anterior. Dacă k-mer-ul nostru interogat se găsește în acest interval, interogarea poate reveni imediat. Dacă k-mer nu este găsit, atunci minimizatorul este calculat; dacă minimizatorul k-mer este același cu ultimul k-mer interogat, atunci interogarea eșuează, deoarece spațiul de căutare al minimizatorului s-a dovedit că nu are k-mer. Numai dacă minimizatorul s-a schimbat, Kraken trebuie să ajusteze intervalul de căutare și să caute din nou k-mer.

construirea metagenomilor simulați

metagenomii HiSeq și MiSeq au fost construiți folosind 20 de seturi de citiri de pușcă cu genom întreg bacterian. Aceste citiri au fost găsite fie ca parte a proiectului GAGE-B, fie în arhiva de citire a secvenței NCBI. Fiecare metagenom conține secvențe din zece genomi (fișier suplimentar 1: Tabelul S1). Atât pentru cele 10.000, cât și pentru 10 milioane de eșantioane citite ale fiecăruia dintre aceste metagenomi, 10% din secvențele lor au fost selectate din fiecare dintre cele zece seturi de date ale genomului component (adică, fiecare genom avea o abundență de secvențe egale). Toate secvențele au fost tăiate pentru a elimina bazele de calitate scăzută și secvențele adaptorului.

compoziția acestor două metagenomi ridică anumite provocări pentru clasificatorii noștri. De exemplu, pelosinus fermentans, găsit în metagenomul nostru HiSeq, nu poate fi identificat corect la nivel de gen de către Kraken (sau oricare dintre ceilalți clasificatori descriși anterior), deoarece nu există genomi Pelosinus în baza de date completă a genomilor RefSeq; cu toate acestea, există șapte astfel de genomi în baza de date Kraken-GB, inclusiv șase tulpini de P. fermentans. În mod similar, în metagenomul nostru MiSeq, Proteus vulgaris este adesea clasificat incorect la nivel de gen, deoarece singurul genom Proteus din Baza de date Kraken este un singur genom Proteus mirabilis. Încă cinci genomi Proteus sunt prezenți în baza de date Kraken-GB, permițând Kraken-GB să clasifice citirile mai bine din acel gen. În plus, Metagenomul MiSeq conține cinci genomi din familia Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus și Salmonella). Similitudinea mare a secvenței dintre genurile din această familie poate face dificilă distincția între genuri pentru orice clasificator.

metagenomul simBA-5 a fost creat prin simularea citirilor din setul de genomi bacterieni și arheali compleți din RefSeq. Repliconii din acele genomi au fost folosiți dacă au fost asociați cu un taxon care avea o intrare asociată cu rangul genului, rezultând un set de repliconi din 607 genuri. Apoi am folosit simulatorul Mason read cu modelul său Illumina pentru a produce 10 milioane de citiri de 100-bp din aceste genomi. Mai întâi am creat genomi simulați pentru fiecare specie, folosind o rată SNP de 0,1% și o rată indel de 0,1% (ambii parametri impliciți), din care am generat citirile. Pentru citirile simulate, am înmulțit ratele implicite de nepotrivire și indel cu cinci, rezultând o rată medie de nepotrivire de 2% (variind de la 1% la începutul citirilor la 6% la capete) și o rată de indel de 1% (0,5% probabilitate de inserție și 0,5% probabilitate de ștergere). Pentru metagenomul simBA-5, setul de 10.000 de citire a fost generat dintr-un eșantion aleatoriu din setul de 10 milioane de citire.

evaluarea acurateței și vitezei

am ales să măsurăm precizia în primul rând la nivelul genului, care a fost cel mai scăzut nivel pentru care am putut determina cu ușurință informațiile taxonomice pentru predicțiile PhymmBL și NBC într-un mod automat. (Acest lucru se datorează modului în care PhymmBL și NBC își raportează rezultatele). Deoarece unele genomi nu au intrări taxonomice la toate cele șapte rânduri (specie, gen, familie, ordine, clasă, filum și regat), am definit sensibilitatea la nivel de gen ca A/B, unde A este numărul de citiri cu un gen atribuit care au fost clasificate corect la acel rang și B este numărul total de citiri cu orice gen atribuit. Am definit sensibilitatea în mod similar pentru alte ranguri taxonomice.

deoarece Kraken poate clasifica o citire la niveluri deasupra speciei, măsurarea preciziei sale ne cere să definim efectul asupra preciziei atribuirii genului corect (de exemplu), în timp ce nu atribuim deloc o specie. Din acest motiv, am definit precizia nivelului de rang ca C/(D + E), unde C este numărul de citiri etichetate la sau sub taxonul corect la rangul măsurat, D este numărul de citiri etichetate la sau sub rangul măsurat și E este numărul de citiri etichetate incorect deasupra rangului măsurat. De exemplu, având în vedere o citire R care ar trebui etichetată ca Escherichia coli, o etichetare a lui R ca E. coli, E. fergusonii sau Escherichia ar îmbunătăți precizia la nivel de gen. O etichetă de Enterobacteriaceae (familie corectă) sau proteobacterii (filum corect) nu ar avea niciun efect asupra preciziei la nivel de gen. O etichetă pentru R de bacil (gen incorect) sau Firmicutes (filum incorect) ar scădea precizia la nivel de gen.

când evaluăm acuratețea PhymmBL, în urma sfaturilor dezvoltatorilor săi , am selectat un prag de încredere al genului pentru comparațiile noastre. Am selectat 3.333 de citiri din setul de date simulat de complexitate medie (simMC), acoperind 31 de genuri diferite. Pentru a simula citirile scurte din datele secvenței Sanger din setul simMC, am selectat ultimele 100 bp din fiecare citire. Apoi am rulat PhymmBL împotriva acelor citiri de 100-bp și am evaluat sensibilitatea la nivel de gen și precizia predicțiilor PhymmBL cu praguri de încredere a genului de la 0 la 1, în trepte de 0,05. Am constatat că un prag de 0,65 a dat cel mai mare scor F (media armonică a sensibilității și preciziei), 0,60 și 0,70 având, de asemenea, scoruri F în 0.5 puncte procentuale din valoarea maximă (fișier suplimentar 1: Tabelul S2). Prin urmare, am folosit pragul de încredere al genului 0.65 în comparațiile noastre. Deși selectarea unui prag depinde de nevoile individuale ale unui utilizator și, prin urmare, este într-o oarecare măsură arbitrară, un prag selectat în acest mod oferă o comparație mai adecvată cu un clasificator selectiv, cum ar fi Kraken, decât niciun prag.

rezultatele de timp și precizie atunci când se utilizează Megablast ca clasificator au fost obținute din datele de jurnal produse de PhymmBL, deoarece PhymmBL folosește Megablast pentru pasul său de aliniere. La atribuirea unei etichete taxonomice la o citire cu Megablast, am folosit taxonul asociat cu prima aliniere raportată. Megablast a fost rulat cu opțiuni implicite.

viteza a fost evaluată utilizând funcționarea cu un singur fir a fiecărui program (cu excepția NBC). PhymmBL a fost modificat astfel încât apelul său la programul blastn a folosit un fir în loc de două. NBC a fost rulat cu 36 de procese concurente care funcționează pe seturi disjuncte de genomi în biblioteca sa genomică, iar timpul total pentru clasificator a fost determinat prin însumarea timpilor de decompresie și notare pentru fiecare genom. Orele de ceas de perete au fost înregistrate pentru toți clasificatorii. În compararea Kraken cu ceilalți clasificatori, am folosit BLAST + 2.2.27, PhymmBL 4.0, NBC 1.1 și MetaPhlAn 1.7.6. Clasificatorii au fost executați pe același computer, cu 48 de procesoare AMD Opteron 6172 2.1 GHz și 252 GB RAM, care rulează Red Hat Enterprise Linux 5. Seturile de date utilizate pentru evaluarea vitezei au avut 10.000 de citiri fiecare pentru toate programele, altele decât Kraken (și variantele sale) și MetaPhlAn, care au folosit 10.000.000 de seturi de date citite. Au fost utilizate numere de citire mai mari cu aceste programe mai rapide pentru a minimiza efectul operațiilor inițiale și finale care au loc în timpul execuției programelor.

deși Kraken este singurul dintre programele pe care le-am examinat care efectuează în mod explicit operațiuni pentru a se asigura că datele sale sunt în memoria fizică înainte de clasificare, am vrut să fim siguri că toate programele au fost evaluate într-un mod similar. La evaluarea vitezei, pentru fiecare program, citim toate fișierele bazei de date (de ex. Fișiere IMM și baze de date BLAST pentru PhymmBL, liste de frecvențe k-mer pentru NBC și indicele Bowtie pentru MetaPhlAn) în memorie de trei ori înainte de a rula programul, pentru a plasa conținutul bazei de date în memoria cache a sistemului de operare (care este stocată în memoria fizică).

dimensiuni reduse ale bazei de date

pentru a genera baza de date de 4 GB pentru rezultatele noastre MiniKraken, am eliminat primele 18 din fiecare bloc de 19 înregistrări din Baza de date standard Kraken. A fost selectat un factor de scădere de 19, deoarece a fost cel mai mic factor întreg care ar reduce dimensiunea la mai puțin de 4 GB, o dimensiune care se poate încadra cu ușurință în memoria multor computere personale obișnuite. Pentru utilizatorii care au mai multă memorie RAM disponibilă, Kraken permite utilizarea unui factor de scădere mai mic, ceea ce va oferi o sensibilitate crescută.

utilizarea genomurilor draft

la construirea bazei de date Kraken-GB, am observat că au existat mai multe contiguri cu secvențe de adaptor cunoscute la capete. În testele ulterioare, am constatat, de asemenea, că unele secvențe din eșantioane cu cantități mari de secvență umană au fost în mod constant clasificate greșit de această bază de date, ceea ce ne-a condus la concluzia că contaminarea a fost probabil prezentă în proiectul genomului. În încercarea de a contracara această contaminare, am eliminat din Baza de date acele k-mers din secvențe de adaptor cunoscute, precum și primul și ultimul 20 K-mers din fiecare dintre contigurile de proiect. Deși acest lucru a îmbunătățit clasificarea, nu a eliminat problema clasificării greșite. Din acest motiv, credem că, dacă genomii de proiect sunt utilizați într-o bază de date Kraken, ar trebui utilizate măsuri foarte stricte pentru a elimina secvențele contaminante din biblioteca genomică.

experimente de excludere Clade

la reanalizarea setului de date simBA-5 pentru experimentele noastre de excludere clade, unele citiri nu au fost utilizate pentru anumite perechi de ranguri măsurate și excluse. Dacă originea unei citiri nu avea o intrare taxonomică la oricare dintre rândurile măsurate sau excluse, aceasta nu a fost utilizată pentru acel experiment special.

în plus, o citire nu a fost utilizată într-un experiment decât dacă cel puțin alți doi taxoni reprezentați în Baza noastră de date (în afară de clada exclusă) la rangul exclus au împărtășit taxonul cladei de origine la rangul măsurat. De exemplu, o citire din genul G nu ar fi utilizată într-un experiment care măsoară precizia la rangul de clasă și exclude rangul genului, cu excepția cazului în care clasa de origine a lui G avea cel puțin alte două genuri cu genomi în biblioteca genomică a lui Kraken. Fără această etapă de filtrare, dacă un gen ar fi exclus atunci când ar fi singurul gen din clasa sa, Kraken nu ar putea numi clasa corectă, deoarece toate intrările din Baza de date din acea clasă ar fi excluse și ele. Aceasta este aceeași abordare luată în experimente similare care au fost utilizate pentru a evalua PhymmBL .

Clasificarea datelor microbiomului uman

am clasificat datele proiectului microbiomului uman folosind o bază de date Kraken realizată din genomi bacterieni, arheali și virali complet RefSeq, împreună cu genomul uman GRCh37. Am recuperat secvențele a trei aderări (SRS019120, SRS014468 și SRS015055) din arhiva de citire a secvenței NCBI și fiecare aderare a avut două runde depuse. Toate citirile au fost tăiate pentru a elimina bazele de calitate scăzută și secvențele adaptorului. Coroana a fost utilizată pentru a genera toate parcelele de distribuție taxonomică.

pentru că secvențele erau toate citiri împerecheate, am unit citirile împreună prin concatenarea perechilor cu o secvență de ‘NNNNN’ între ele. Kraken ignoră k-mers cu nucleotide ambigue, astfel încât k-mers care acoperă aceste ‘N’ caractere nu afectează clasificarea. Această operație i-a permis lui Kraken să clasifice o pereche de lecturi ca o singură unitate, mai degrabă decât să trebuiască să clasifice colegii separat.

disponibilitatea Software-ului și a datelor