Kraken: classificazione delle sequenze metagenomiche ultraveloci con allineamenti esatti

Algoritmo di classificazione delle sequenze

Per classificare una sequenza di DNA S, raccogliamo tutti i k-mer all’interno di quella sequenza in un insieme, indicato come K(S). Quindi mappiamo ogni k-mer in K (S), usando l’algoritmo descritto di seguito, al taxon LCA di tutti i genomi che contengono quel k-mer. Questi taxa LCA e i loro antenati nell’albero della tassonomia formano ciò che chiamiamo albero di classificazione, un sottoalbero potato che viene utilizzato per classificare S. Ogni nodo nell’albero di classificazione è ponderato con il numero di k-mers in K(S) mappato al taxon associato a quel nodo. Quindi, ogni percorso RTL (root-to-leaf) nell’albero di classificazione viene valutato calcolando la somma di tutti i pesi dei nodi lungo il percorso. Il percorso RTL di punteggio massimo nell’albero di classificazione è il percorso di classificazione e a S viene assegnata l’etichetta corrispondente alla sua foglia (se ci sono più percorsi di punteggio massimo, viene selezionata l’LCA di tutte le foglie di quei percorsi). Questo algoritmo, illustrato in Figura 1, permette Kraken di considerare ogni k-mer all’interno di una sequenza come un pezzo separato di prova, e quindi tentare di risolvere eventuali prove in conflitto, se necessario. Si noti che per una scelta appropriata di k, la maggior parte dei k-mer mapperà in modo univoco a una singola specie, semplificando notevolmente il processo di classificazione. Sequenze per le quali nessuno dei k-mers in K(S) si trovano in qualsiasi genoma sono lasciati non classificati da questo algoritmo.

L’uso di RTL path scoring nell’albero di classificazione è necessario alla luce delle inevitabili differenze tra le sequenze da classificare e le sequenze presenti in qualsiasi libreria di genomi. Tali differenze possono, anche per grandi valori di k, portare a un k-mer che è presente nella libreria ma associato a una specie lontana dalla vera specie di origine. Segnando i vari percorsi RTL nell’albero di classificazione, possiamo compensare queste differenze e classificare correttamente le sequenze anche quando una piccola minoranza di k-mer in una sequenza indica che alla sequenza dovrebbe essere assegnata un’etichetta tassonomica errata.

Creazione di database

L’implementazione efficiente dell’algoritmo di classificazione di Kraken richiede che la mappatura di k-mers a taxa venga eseguita interrogando un database pre-calcolato. Kraken crea questo database attraverso un processo multi-step, iniziando con la selezione di una libreria di sequenze genomiche. Kraken include una libreria predefinita, basata su genomi microbici completati nel database RefSeq del National Center for Biotechnology Information (NCBI), ma la libreria può essere personalizzata in base alle esigenze dei singoli utenti .

Una volta scelta la libreria, usiamo il contatore k-mer multithread Jellyfish per creare un database contenente ogni 31-mer distinto nella libreria. Una volta completato il database, gli spazi di 4 byte Jellyfish utilizzati per memorizzare i conteggi k-mer nel file di database vengono invece utilizzati da Kraken per memorizzare i numeri ID tassonomici dei valori LCA dei k-mer. Dopo che il database è stato creato da Jellyfish, le sequenze genomiche nella libreria vengono elaborate una alla volta. Per ogni sequenza, il taxon ad esso associato viene utilizzato per impostare i valori LCA memorizzati di tutti i k-mer nella sequenza. Mentre le sequenze vengono elaborate, se un k-mer da una sequenza ha avuto il suo valore LCA precedentemente impostato, allora l’LCA del valore memorizzato e il taxon della sequenza corrente viene calcolato e che LCA è memorizzato per il k-mer. Informazioni Taxon è ottenuto dal database NCBI tassonomia.

Struttura del database e algoritmo di ricerca

Poiché Kraken utilizza molto frequentemente un k-mer come query di database immediatamente dopo aver interrogato un k-mer adiacente e poiché i k-mer adiacenti condividono una notevole quantità di sequenza, utilizziamo il concetto di minimizer per raggruppare k-mer simili insieme. Per spiegare la nostra applicazione di questo concetto, qui definiamo la rappresentazione canonica di una sequenza di DNA S come lessicograficamente più piccola di S e il complemento inverso di S. Per determinare il minimizzatore di lunghezza M di un k-mer, consideriamo la rappresentazione canonica di tutte le M-mer nel k-mer e selezioniamo il lessicograficamente più piccolo di quelle M-mer come minimizzatore del k-mer. In pratica, i k-mer adiacenti avranno spesso lo stesso minimizzatore.

Nel database di Kraken, tutti i k-mer con lo stesso minimizzatore sono memorizzati consecutivamente e sono ordinati in ordine lessicografico delle loro rappresentazioni canoniche. Una query per un k-mer R può quindi essere elaborata cercando in un indice le posizioni nel database in cui verranno memorizzati i k-mer con il minimizzatore di R e quindi eseguendo una ricerca binaria all’interno di quella regione (Figura 5). Poiché i k-mer adiacenti hanno spesso lo stesso minimizzatore, l’intervallo di ricerca è spesso lo stesso tra due query consecutive e la ricerca nella prima query può spesso portare dati nella cache della CPU che verrà utilizzata nella seconda query. Consentendo agli accessi alla memoria nelle query successive di accedere ai dati nella cache della CPU anziché nella RAM, questa strategia rende le query successive molto più veloci di quanto non sarebbero altrimenti. L’indice contenente gli offset di ciascun gruppo di k-mer nel database richiede 8 × 4 M byte. Per impostazione predefinita Kraken utilizza 15-bp minimizzatori, ma l’utente può modificare questo valore; per esempio, nella creazione di MiniKraken, abbiamo usato 13-bp minimizzatori per garantire la dimensione totale del database rimasto sotto 4 GB.

Struttura del database Kraken. Ogni k-mer da interrogare sul database ha una sottostringa specifica che è il suo minimizzatore. Per cercare un k-mer nel database, vengono esaminate le posizioni nel database che contengono k-mer con lo stesso minimizzatore. Queste posizioni vengono trovate rapidamente esaminando l’array di offset minimizzatore per le posizioni iniziali dei record con il minimizzatore di k-mer (arancione) e il successivo minimizzatore possibile (blu). All’interno di un intervallo di record associati a un dato minimizzatore, i record sono ordinati in base all’ordine lessicografico dei loro k-mer, consentendo di completare una query utilizzando una ricerca binaria su questo intervallo.

Nell’implementazione di Kraken, abbiamo apportato ulteriori ottimizzazioni alla struttura e all’algoritmo di ricerca descritti sopra. In primo luogo, come notato da Roberts et al. , un semplice ordinamento lessicografico di M-mers può risultare in una distribuzione distorta di minimizzatori che sovrarappresenta M-mers a bassa complessità. In Kraken, un tale pregiudizio creerebbe molti ampi intervalli di ricerca, che richiederebbero più tempo per la ricerca. Per creare una distribuzione più uniforme dei minimizzatori (e quindi accelerare le ricerche), usiamo l’operazione exclusive-or (XOR) per alternare metà dei bit della rappresentazione canonica di ogni M-mer prima di confrontare le M-mer tra loro usando l’ordinamento lessicografico. Questa operazione XOR rimescola efficacemente l’ordinamento standard e impedisce la grande polarizzazione verso minimizzatori a bassa complessità.

Sfruttiamo anche il fatto che l’intervallo di ricerca è spesso lo stesso tra le query per rendere più veloci le query di Kraken. Invece di calcolare il minimizzatore ogni volta che eseguiamo una query, cerchiamo prima l’intervallo precedente. Se il nostro k-mer interrogato si trova in questo intervallo, la query può tornare immediatamente. Se il k-mer non viene trovato, viene calcolato il minimizzatore; se il minimizzatore del k-mer è lo stesso dell’ultimo k-mer interrogato, la query fallisce, poiché lo spazio di ricerca del minimizzatore ha dimostrato di non avere il k-mer. Solo se il minimizzatore è cambiato Kraken deve regolare l’intervallo di ricerca e cercare nuovamente il k-mer.

Costruzione di metagenomi simulati

I metagenomi HiSeq e MiSeq sono stati costruiti utilizzando 20 set di shotgun batterici a genoma intero. Queste letture sono state trovate come parte del progetto GAGE-B o nell’archivio di lettura sequenza NCBI. Ogni metagenoma contiene sequenze da dieci genomi (file aggiuntivo 1: Tabella S1). Per entrambi i 10.000 e 10 milioni di campioni letti di ciascuno di questi metagenomi, il 10% delle loro sequenze è stato selezionato da ciascuno dei dieci set di dati del genoma componente (cioè, ogni genoma aveva uguale abbondanza di sequenza). Tutte le sequenze sono state tagliate per rimuovere basi di bassa qualità e sequenze di adattatori.

La composizione di questi due metagenomi pone alcune sfide ai nostri classificatori. Ad esempio, Pelosinus fermentans, trovato nel nostro metagenoma HiSeq, non può essere correttamente identificato a livello di genere da Kraken (o da uno qualsiasi degli altri classificatori precedentemente descritti), perché non ci sono genomi Pelosinus nel database Refseq complete genomes; tuttavia, ci sono sette genomi di questo tipo nel database di Kraken-GB, inclusi sei ceppi di P. fermentans. Allo stesso modo, nel nostro metagenoma MiSeq, Proteus vulgaris è spesso classificato erroneamente a livello di genere perché l’unico genoma Proteus nel database di Kraken è un singolo genoma Proteus mirabilis. Altri cinque genomi Proteus sono presenti nel database di Kraken-GB, consentendo a Kraken-GB di classificare meglio le letture da quel genere. Inoltre, il metagenoma MiSeq contiene cinque genomi della famiglia delle Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus e Salmonella). L’elevata somiglianza di sequenza tra i generi di questa famiglia può rendere difficile la distinzione tra generi per qualsiasi classificatore.

Il metagenoma simBA-5 è stato creato simulando le letture dall’insieme di genomi batterici e archaeal completi in RefSeq. I repliconi di questi genomi sono stati utilizzati se sono stati associati a un taxon che aveva una voce associata al rango del genere, risultando in un insieme di repliconi da 607 generi. Abbiamo quindi utilizzato il Mason read simulator con il suo modello Illumina per produrre 10 milioni di letture da 100 bp da questi genomi. Per prima cosa abbiamo creato genomi simulati per ogni specie, utilizzando un tasso SNP dello 0,1% e un tasso indel dello 0,1% (entrambi i parametri predefiniti), da cui abbiamo generato le letture. Per le letture simulate, abbiamo moltiplicato il disallineamento predefinito e i tassi indel per cinque, risultando in un tasso medio di disallineamento del 2% (che va dall ‘1% all’inizio delle letture al 6% alle estremità) e un tasso indel dell’ 1% (0,5% di probabilità di inserimento e 0,5% di probabilità di cancellazione). Per il metagenoma simBA-5, il set di 10.000 letture è stato generato da un campione casuale del set di 10 milioni di letture.

Valutazione dell’accuratezza e della velocità

Abbiamo scelto di misurare l’accuratezza principalmente a livello di genere, che era il livello più basso per il quale potevamo facilmente determinare le informazioni sulla tassonomia per le previsioni di PhymmBL e NBC in modo automatizzato. (Ciò è dovuto al modo in cui PhymmBL e NBC riportano i loro risultati). Poiché alcuni genomi non hanno voci tassonomiche in tutti e sette i ranghi (specie, genere, famiglia, ordine, classe, phylum e regno), abbiamo definito la sensibilità a livello di genere come A/B, dove A è il numero di letture con un genere assegnato che sono state correttamente classificate in quel rango, e B è il numero totale di letture con qualsiasi genere assegnato. Abbiamo definito la sensibilità in modo simile per altri ranghi tassonomici.

Poiché Kraken può classificare una lettura a livelli superiori alla specie, misurare la sua precisione ci richiede di definire l’effetto sulla precisione dell’assegnazione del genere corretto (ad esempio) senza assegnare affatto una specie. Per questo motivo, abbiamo definito la precisione a livello di rango come C / (D + E), dove C è il numero di letture etichettate in corrispondenza o al di sotto del taxon corretto al rango misurato, D è il numero di letture etichettate in corrispondenza o al di sotto del rango misurato, ed E è il numero di letture etichettate erroneamente al di sopra del rango misurato. Ad esempio, data una lettura R che dovrebbe essere etichettata come Escherichia coli, un’etichettatura di R come E. coli, E. fergusonii o Escherichia migliorerebbe la precisione a livello di genere. Un’etichetta di Enterobacteriaceae (famiglia corretta) o Proteobacteria (phylum corretto) non avrebbe alcun effetto sulla precisione a livello di genere. Un’etichetta per R di Bacillus (genere errato) o Firmicutes (phylum errato) diminuirebbe la precisione a livello di genere.

Nel valutare l’accuratezza di PhymmBL, seguendo il consiglio dei suoi sviluppatori , abbiamo selezionato una soglia di confidenza del genere per i nostri confronti. Abbiamo selezionato 3.333 letture dal set di dati simulato media complessità (simMC), che copre 31 generi diversi. Per simulare brevi letture dai dati della sequenza Sanger nel set simMC, abbiamo selezionato gli ultimi 100 bp da ciascuna delle letture. Abbiamo quindi eseguito PhymmBL contro quelle letture da 100 bp e valutato la sensibilità a livello di genere e la precisione delle previsioni di PhymmBL con soglie di confidenza del genere da 0 a 1, con incrementi di 0.05. Abbiamo scoperto che una soglia di 0.65 ha prodotto il punteggio F più alto (la media armonica di sensibilità e precisione), con 0.60 e 0.70 che hanno anche punteggi F entro 0.5 punti percentuali del massimo (file aggiuntivo 1: Tabella S2). Abbiamo quindi utilizzato la soglia di confidenza del genere 0.65 nei nostri confronti. Sebbene la selezione di una soglia dipenda dalle esigenze individuali di un utente, e quindi sia in una certa misura arbitraria, una soglia selezionata in questo modo fornisce un confronto più corretto con un classificatore selettivo come Kraken rispetto a nessuna soglia.

I risultati di tempo e precisione quando si utilizza Megablast come classificatore sono stati ottenuti dai dati di log prodotti da PhymmBL, poiché PhymmBL utilizza Megablast per la sua fase di allineamento. Quando si assegna un’etichetta tassonomica a una lettura con Megablast, abbiamo usato il taxon associato al primo allineamento riportato. Megablast è stato eseguito con opzioni predefinite.

La velocità è stata valutata utilizzando l’operazione a thread singolo di ciascun programma (ad eccezione di NBC). PhymmBL è stato modificato in modo che la sua chiamata al programma blastn utilizzasse un thread invece di due. NBC è stato eseguito con 36 processi simultanei che operano su insiemi disgiunti di genomi nella sua libreria genomica, e il tempo totale per il classificatore è stato determinato sommando i tempi di decompressione e punteggio per ogni genoma. I tempi dell’orologio da parete sono stati registrati per tutti i classificatori. Nel confrontare Kraken con gli altri classificatori, abbiamo usato BLAST+ 2.2.27, PhymmBL 4.0, NBC 1.1 e MetaPhlAn 1.7.6. I classificatori erano tutti eseguiti sullo stesso computer, con 48 CPU AMD Opteron 6172 2.1 GHz e 252 GB di RAM, con Red Hat Enterprise Linux 5. I set di dati utilizzati per la valutazione della velocità avevano 10.000 letture ciascuno per tutti i programmi diversi da Kraken (e le sue varianti) e MetaPhlAn, che utilizzavano 10.000.000 di set di dati letti. Numeri di lettura più elevati sono stati utilizzati con questi programmi più veloci per ridurre al minimo l’effetto delle operazioni iniziali e finali che avvengono durante l’esecuzione dei programmi.

Sebbene Kraken sia l’unico dei programmi che abbiamo esaminato che esegue esplicitamente operazioni per garantire che i suoi dati siano nella memoria fisica prima della classificazione, volevamo essere sicuri che tutti i programmi fossero valutati in modo simile. Nel valutare la velocità, per ogni programma, leggiamo tutti i file di database (ad es. File IMM e database BLAST per PhymmBL, liste di frequenza k-mer per NBC e l’indice Bowtie per MetaPhlAn) in memoria tre volte prima di eseguire il programma, al fine di inserire il contenuto del database nella cache del sistema operativo (che viene memorizzato nella memoria fisica).

Dimensioni ridotte del database

Per generare il database da 4 GB per i nostri risultati MiniKraken, abbiamo rimosso i primi 18 di ogni blocco di 19 record nel database Kraken standard. Un fattore di contrazione di 19 è stato selezionato come era il più piccolo fattore intero che ridurrebbe la dimensione a meno di 4 GB, una dimensione che può facilmente inserirsi nella memoria di molti personal computer comuni. Per gli utenti che hanno più RAM disponibile, Kraken consente di utilizzare un fattore di restringimento più piccolo, che darà una maggiore sensibilità.

Uso di progetti di genomi

Durante la costruzione del database Kraken-GB, abbiamo notato che c’erano diversi contig con sequenze di adattatori note alle estremità. Nei test successivi, abbiamo anche scoperto che alcune sequenze in campioni con grandi quantità di sequenza umana sono state costantemente classificate erroneamente da questo database, portandoci a concludere che la contaminazione era probabilmente presente nel progetto di genomi. Nel tentativo di contrastare questa contaminazione, abbiamo rimosso dal database quei k-mer da sequenze di adattatori note, così come il primo e l’ultimo 20 k-mer da ciascuno dei draft contig. Sebbene ciò abbia migliorato la classificazione, non ha eliminato il problema di errata classificazione. Per questo motivo, riteniamo che se i progetti di genomi vengono utilizzati in un database Kraken, dovrebbero essere utilizzate misure molto severe per rimuovere le sequenze di contaminanti dalla libreria genomica.

Esperimenti di esclusione del Clade

Quando si analizza nuovamente il set di dati simBA-5 per i nostri esperimenti di esclusione del clade, alcune letture non sono state utilizzate per alcune coppie di ranghi misurati ed esclusi. Se l’origine di una lettura mancava di una voce tassonomica in uno dei ranghi misurati o esclusi, non è stata utilizzata per quel particolare esperimento.

Inoltre, una lettura non è stata utilizzata in un esperimento a meno che almeno altri due taxa rappresentati nel nostro database (a parte il clade escluso) al rango escluso condividessero il taxon del clade di origine al rango misurato. Ad esempio, una lettura dal genere G non sarebbe stata utilizzata in un esperimento che misurava la precisione al rango di classe ed escludeva il rango di genere a meno che la classe di origine di G non avesse almeno altri due generi con genomi nella libreria genomica di Kraken. Senza questo passaggio di filtraggio, se un genere fosse escluso quando era l’unico genere della sua classe, Kraken non poteva nominare la classe corretta, poiché anche tutte le voci nel database di quella classe sarebbero state escluse. Questo è lo stesso approccio adottato in esperimenti simili che sono stati utilizzati per valutare PhymmBL .

Classificazione dei dati del microbioma umano

Abbiamo classificato i dati del progetto del microbioma umano utilizzando un database Kraken composto da genomi batterici, archaeali e virali RefSeq completi, insieme al genoma umano GRCh37. Abbiamo recuperato le sequenze di tre adesioni (SRS019120, SRS014468 e SRS015055) dall’archivio di lettura sequenza NCBI e ogni adesione ha presentato due esecuzioni. Tutte le letture sono state tagliate per rimuovere basi di bassa qualità e sequenze di adattatori. Krona è stato utilizzato per generare tutte le trame di distribuzione tassonomica.

Poiché le sequenze erano tutte letture accoppiate, abbiamo unito le letture insieme concatenando i compagni con una sequenza di ‘NNNNN’ tra di loro. Kraken ignora k-mer con nucleotidi ambigui, quindi i k-mer che coprono questi caratteri ‘N’ non influenzano la classificazione. Questa operazione ha permesso a Kraken di classificare una coppia di letture come una singola unità piuttosto che dover classificare i compagni separatamente.

Disponibilità di software e dati