Kraken: ultrafast metagenomic sequence classification using exact alignments

Sequence classification algorithm

Um eine DNA-Sequenz zu klassifizieren, sammeln wir alle k-mers innerhalb dieser Sequenz in einer Menge, die als K(S) bezeichnet wird. Wir ordnen dann jedes k-Mer in K (S) unter Verwendung des unten beschriebenen Algorithmus dem LCA-Taxon aller Genome zu, die dieses k-Mer enthalten. Diese Ökobilanztaxa und ihre Vorfahren im Taxonomiebaum bilden den sogenannten Klassifikationsbaum, einen beschnittenen Teilbaum, der zur Klassifizierung von S verwendet wird. Anschließend wird jeder RTL-Pfad (Root-to-Leaf) im Klassifikationsbaum bewertet, indem die Summe aller Knotengewichte entlang des Pfads berechnet wird. Der maximale Bewertungspfad im Klassifikationsbaum ist der Klassifikationspfad, und S wird die Bezeichnung zugewiesen, die seinem Blatt entspricht (wenn mehrere Pfade mit maximaler Bewertung vorhanden sind, wird die Ökobilanz aller Blätter dieser Pfade ausgewählt). Dieser in Abbildung 1 dargestellte Algorithmus ermöglicht es Kraken, jedes k-Mer innerhalb einer Sequenz als separates Beweisstück zu betrachten und dann bei Bedarf zu versuchen, widersprüchliche Beweise aufzulösen. Beachten Sie, dass für eine geeignete Wahl von k, Die meisten k-mers werden eindeutig einer einzelnen Art zugeordnet, Dies vereinfacht den Klassifizierungsprozess erheblich. Sequenzen, für die keines der k-Mere in K (S) in einem Genom gefunden wird, werden von diesem Algorithmus nicht klassifiziert.

Die Verwendung einer Pfadbewertung im Klassifikationsbaum ist angesichts der unvermeidlichen Unterschiede zwischen den zu klassifizierenden Sequenzen und den in einer beliebigen Genombibliothek vorhandenen Sequenzen erforderlich. Solche Unterschiede können selbst bei großen Werten von k zu einem k-mer führen, der in der Bibliothek vorhanden ist, aber einer Spezies zugeordnet ist, die weit von der wahren Quellspezies entfernt ist. Durch die Bewertung der verschiedenen RTL-Pfade im Klassifikationsbaum können wir diese Unterschiede ausgleichen und Sequenzen korrekt klassifizieren, selbst wenn eine kleine Minderheit von k-Mern in einer Sequenz angibt, dass der Sequenz eine falsche taxonomische Bezeichnung zugewiesen werden sollte.

Datenbankerstellung

Die effiziente Implementierung des Kraken-Klassifizierungsalgorithmus erfordert, dass die Zuordnung von k-mers zu Taxa durch Abfragen einer vorberechneten Datenbank erfolgt. Kraken erstellt diese Datenbank in einem mehrstufigen Prozess, beginnend mit der Auswahl einer Bibliothek genomischer Sequenzen. Kraken enthält eine Standardbibliothek, die auf abgeschlossenen mikrobiellen Genomen in der RefSeq-Datenbank des National Center for Biotechnology Information (NCBI) basiert, aber die Bibliothek kann von einzelnen Benutzern nach Bedarf angepasst werden .

Sobald die Bibliothek ausgewählt ist, verwenden wir den Quallen-Multithread-k-Mer-Zähler, um eine Datenbank zu erstellen, die alle eindeutigen 31-Mer in der Bibliothek enthält. Sobald die Datenbank vollständig ist, werden die 4-Byte-Leerzeichen, die Quallen zum Speichern der k-Mer-Zählungen in der Datenbankdatei verwenden, stattdessen von Kraken zum Speichern der taxonomischen ID-Nummern der Ökobilanzwerte der k-mers verwendet. Nachdem die Datenbank von Jellyfish erstellt wurde, werden die genomischen Sequenzen in der Bibliothek einzeln verarbeitet. Für jede Sequenz wird das ihr zugeordnete Taxon verwendet, um die gespeicherten Ökobilanzwerte aller k-Mere in der Sequenz festzulegen. Wenn bei der Verarbeitung von Sequenzen für ein k-Mer aus einer Sequenz der LCA-Wert zuvor festgelegt wurde, wird die LCA des gespeicherten Werts und des Taxons der aktuellen Sequenz berechnet und diese LCA für das k-Mer gespeichert. Taxon Informationen werden von der NCBI Taxonomie-Datenbank erhalten.

Datenbankstruktur und Suchalgorithmus

Da Kraken sehr häufig einen K-Mer als Datenbankabfrage unmittelbar nach der Abfrage eines benachbarten k-Mers verwendet und benachbarte k-Mer eine erhebliche Menge an Sequenz gemeinsam haben, verwenden wir das Minimizer-Konzept, um ähnliche k-Mer zu gruppieren. Um unsere Anwendung dieses Konzepts zu erklären, definieren wir hier die kanonische Darstellung einer DNA-Sequenz S als die lexikographisch kleinere von S und das umgekehrte Komplement von S. Um einen k-Mer-Minimizer der Länge M zu bestimmen, betrachten wir die kanonische Darstellung aller M-mer im k-mer und wählen den lexikographisch kleinsten dieser M-mer als Minimizer des k-mer aus. In der Praxis haben benachbarte k-mer oft den gleichen Minimierer.

In Krakens Datenbank werden alle K-mer mit demselben Minimierer nacheinander gespeichert und in lexikographischer Reihenfolge ihrer kanonischen Darstellungen sortiert. Eine Abfrage nach einem k-mer R kann dann verarbeitet werden, indem in einem Index die Positionen in der Datenbank nachgeschlagen werden, an denen die k-mer mit dem Minimizer von R gespeichert würden, und dann eine binäre Suche innerhalb dieser Region durchgeführt wird (Abbildung 5). Da benachbarte k-Mer häufig denselben Minimierer haben, ist der Suchbereich zwischen zwei aufeinanderfolgenden Abfragen häufig derselbe, und die Suche in der ersten Abfrage kann häufig Daten in den CPU-Cache bringen, die in der zweiten Abfrage verwendet werden. Indem Speicherzugriffe in nachfolgenden Abfragen auf Daten im CPU-Cache anstelle von RAM zugreifen können, werden nachfolgende Abfragen mit dieser Strategie viel schneller ausgeführt, als dies sonst der Fall wäre. Der Index, der die Offsets jeder Gruppe von k-mers in der Datenbank enthält, benötigt 8 × 4M Bytes. Standardmäßig verwendet Kraken 15-bp-Minimierer, aber der Benutzer kann diesen Wert ändern; Bei der Erstellung von MiniKraken haben wir beispielsweise 13-bp-Minimierer verwendet, um sicherzustellen, dass die Gesamtdatenbankgröße unter 4 GB bleibt.

Kraken Datenbankstruktur. Jeder k-Mer, der gegen die Datenbank abgefragt werden soll, hat einen bestimmten Teilstring, der sein Minimierer ist. Um nach einem k-Mer in der Datenbank zu suchen, werden die Positionen in der Datenbank untersucht, die k-Mer mit demselben Minimierer enthalten. Diese Positionen werden schnell gefunden, indem das Minimizer-Offset-Array auf die Startpositionen von Datensätzen mit dem Minimizer des k-mer (orange) und dem nächstmöglichen Minimizer (blau) untersucht wird. Innerhalb eines Bereichs von Datensätzen, die einem bestimmten Minimierer zugeordnet sind, werden Datensätze nach lexikographischer Reihenfolge ihrer k-mers sortiert, sodass eine Abfrage mithilfe einer binären Suche über diesen Bereich abgeschlossen werden kann.

Bei der Implementierung von Kraken haben wir die oben beschriebene Struktur und den Suchalgorithmus weiter optimiert. Erstens, wie von Roberts et al. , eine einfache lexikographische Ordnung von M-mers kann zu einer schiefen Verteilung von Minimierern führen, die M-Mers mit geringer Komplexität überrepräsentieren. In Kraken würde eine solche Verzerrung viele große Suchbereiche erzeugen, was mehr Zeit für die Suche erfordern würde. Um eine gleichmäßigere Verteilung der Minimierer zu erstellen (und somit die Suche zu beschleunigen), verwenden wir die exklusive-or (XOR) -Operation, um die Hälfte der Bits der kanonischen Darstellung jedes M-Mer umzuschalten, bevor die M-mer mithilfe der lexikografischen Reihenfolge miteinander verglichen werden. Diese XOR-Operation verschlüsselt effektiv die Standardreihenfolge und verhindert die große Tendenz zu Minimierern mit geringer Komplexität.

Wir nutzen auch die Tatsache, dass der Suchbereich zwischen Abfragen oft gleich ist, um Krakens Abfragen schneller zu machen. Anstatt den Minimierer bei jeder Abfrage zu berechnen, durchsuchen wir zuerst den vorherigen Bereich. Wenn unser abgefragter k-mer in diesem Bereich gefunden wird, kann die Abfrage sofort zurückkehren. Wenn der Minimizer des k-Mer mit dem zuletzt abgefragten k-mer übereinstimmt, schlägt die Abfrage fehl, da gezeigt wurde, dass der Suchraum des Minimizers keinen k-Mer aufweist. Nur wenn sich der Minimizer geändert hat, muss Kraken den Suchbereich anpassen und erneut nach dem k-Mer suchen.

Konstruktion simulierter Metagenome

Die HiSeq- und MiSeq-Metagenome wurden unter Verwendung von 20 Sätzen bakterieller Ganzgenom-Shotgun-Lesevorgänge erstellt. Diese Lesevorgänge wurden entweder im Rahmen des GAGE-B-Projekts oder im NCBI Sequence Read Archive gefunden. Jedes Metagenom enthält Sequenzen aus zehn Genomen (Zusatzdatei 1: Tabelle S1). Sowohl für die 10.000 als auch für die 10 Millionen Leseproben jedes dieser Metagenome wurden 10% ihrer Sequenzen aus jedem der zehn Komponentengenomdatensätze ausgewählt (d. h. Jedes Genom hatte die gleiche Sequenzhäufigkeit). Alle Sequenzen wurden getrimmt, um minderwertige Basen und Adaptersequenzen zu entfernen.

Die Zusammensetzung dieser beiden Metagenome stellt unsere Klassifikatoren vor gewisse Herausforderungen. Zum Beispiel kann Pelosinus fermentans, der in unserem HiSeq-Metagenom gefunden wird, auf Gattungsebene von Kraken (oder einem der anderen zuvor beschriebenen Klassifikatoren) nicht korrekt identifiziert werden, da es in der vollständigen RefSeq-Genomdatenbank keine Pelosinus-Genome gibt; Es gibt jedoch sieben solcher Genome in der Datenbank von Kraken-GB, darunter sechs Stämme von P. fermentans. In ähnlicher Weise wird Proteus vulgaris in unserem MiSeq-Metagenom auf Gattungsebene häufig falsch klassifiziert, da das einzige Proteus-Genom in Krakens Datenbank ein einzelnes Proteus mirabilis-Genom ist. Fünf weitere Proteus-Genome sind in der Datenbank von Kraken-GB vorhanden, so dass Kraken-GB Tiere dieser Gattung besser klassifizieren kann. Darüber hinaus enthält das MiSeq-Metagenom fünf Genome aus der Familie der Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus und Salmonella). Die hohe Sequenzähnlichkeit zwischen den Gattungen in dieser Familie kann die Unterscheidung zwischen Gattungen für jeden Klassifikator schwierig machen.

Das simBA-5-Metagenom wurde durch Simulation von Lesevorgängen aus dem Satz vollständiger bakterieller und archaealer Genome in RefSeq erstellt. Replikone aus diesen Genomen wurden verwendet, wenn sie mit einem Taxon assoziiert waren, das einen Eintrag mit dem Gattungsrang hatte, was zu einem Satz von Replikonen aus 607 Gattungen führte. Wir haben dann den Mason Read Simulator mit seinem Illumina-Modell verwendet, um 10 Millionen 100-bp-Lesevorgänge aus diesen Genomen zu erzeugen. Zuerst erstellten wir simulierte Genome für jede Spezies mit einer SNP-Rate von 0,1% und einer Indel-Rate von 0,1% (beide Standardparameter), aus denen wir die Lesevorgänge generierten. Für die simulierten Lesevorgänge multiplizierten wir die Standard-Mismatch- und Indel-Raten mit fünf, was zu einer durchschnittlichen Mismatch-Rate von 2% (von 1% zu Beginn der Lesevorgänge bis zu 6% am Ende) und einer Indel-Rate von 1% (0,5% Einfügewahrscheinlichkeit und 0,5% Löschwahrscheinlichkeit) führte. Für das simBA-5-Metagenom wurde der 10.000-Lesesatz aus einer Zufallsstichprobe des 10-Millionen-Lesesatzes generiert.

Bewertung der Genauigkeit und Geschwindigkeit

Wir haben uns dafür entschieden, die Genauigkeit hauptsächlich auf der Gattungsebene zu messen, der niedrigsten Ebene, für die wir die Taxonomieinformationen für die Vorhersagen von PhymmBL und NBC auf automatisierte Weise leicht bestimmen konnten. (Dies liegt an der Art und Weise, in der PhymmBL und NBC ihre Ergebnisse melden). Da einige Genome keine taxonomischen Einträge in allen sieben Rängen (Spezies, Gattung, Familie, Ordnung, Klasse, Stamm und Königreich) haben, haben wir die Empfindlichkeit auf Gattungsebene als A / B definiert, wobei A die Anzahl der Lesevorgänge mit einer zugewiesenen Gattung ist, die korrekt klassifiziert wurden In diesem Rang und B ist die Gesamtzahl der Lesevorgänge mit einer zugewiesenen Gattung. Wir haben die Empfindlichkeit für andere taxonomische Ränge ähnlich definiert.

Da Kraken eine Gattung auf Ebenen oberhalb der Art klassifizieren kann, müssen wir bei der Messung ihrer Genauigkeit den Effekt auf die Genauigkeit definieren, wenn wir beispielsweise die richtige Gattung zuweisen, während wir überhaupt keine Art zuweisen. Aus diesem Grund haben wir die Genauigkeit auf Rangebene als C / (D + E) definiert, wobei C die Anzahl der Lesevorgänge ist, die am oder unter dem richtigen Taxon am gemessenen Rang markiert sind, D die Anzahl der Lesevorgänge ist, die am oder unter dem gemessenen Rang markiert sind, und E ist die Anzahl der Lesevorgänge, die falsch über dem gemessenen Rang markiert sind. Zum Beispiel würde eine Markierung von R als E. coli, E. fergusonii oder Escherichia die Genauigkeit auf Gattungsebene verbessern. Eine Markierung von Enterobacteriaceae (korrekte Familie) oder Proteobacteria (korrektes Phylum) hätte keinen Einfluss auf die Genauigkeit auf Gattungsebene. Eine Markierung für R von Bacillus (falsche Gattung) oder Firmicutes (falsches Phylum) würde die Genauigkeit auf Gattungsebene verringern.

Bei der Bewertung der Genauigkeit von PhymmBL haben wir nach dem Rat der Entwickler eine bestimmte Vertrauensschwelle für unsere Vergleiche ausgewählt. Wir wählten 3.333 Lesevorgänge aus dem Datensatz simulierte mittlere Komplexität (simMC) aus, der 31 verschiedene Gattungen abdeckt. Um kurze Lesevorgänge aus den Sanger-Sequenzdaten im simMC-Set zu simulieren, haben wir die letzten 100 bp aus jedem der Lesevorgänge ausgewählt. Wir liefen dann PhymmBL gegen diese 100-bp liest, und bewertete die Genus-Ebene Empfindlichkeit und Präzision der PhymmBL Vorhersagen mit Genus Vertrauen Schwellenwerte von 0 bis 1, in Schritten von 0,05. Wir fanden heraus, dass ein Schwellenwert von 0,65 den höchsten F-Score (das harmonische Mittel von Empfindlichkeit und Präzision) ergab, wobei 0,60 und 0,70 ebenfalls F-Scores innerhalb von 0 aufwiesen.5 Prozentpunkte des Maximums (Zusatzdatei 1: Tabelle S2). Wir haben daher in unseren Vergleichen die Konfidenzschwelle von 0,65 % verwendet. Obwohl die Auswahl eines Schwellenwerts von den individuellen Bedürfnissen eines Benutzers abhängt und daher in gewissem Maße willkürlich ist, bietet ein auf diese Weise ausgewählter Schwellenwert einen angemesseneren Vergleich mit einem selektiven Klassifikator wie Kraken als überhaupt kein Schwellenwert.

Die Zeit- und Genauigkeitsergebnisse bei Verwendung von Megablast als Klassifikator wurden aus den von PhymmBL erzeugten Protokolldaten erhalten, da PhymmBL Megablast für seinen Ausrichtungsschritt verwendet. Beim Zuweisen einer taxonomischen Bezeichnung zu einem Lesevorgang mit Megablast haben wir das Taxon verwendet, das der ersten gemeldeten Ausrichtung zugeordnet ist. Megablast wurde mit Standardoptionen ausgeführt.

Die Geschwindigkeit wurde anhand der Single-Threaded-Operation jedes Programms (mit Ausnahme von NBC) bewertet. PhymmBL wurde so geändert, dass sein Aufruf des blastn-Programms einen Thread anstelle von zwei verwendete. NBC wurde mit 36 gleichzeitigen Prozessen ausgeführt, die auf disjunkten Genomsätzen in seiner genomischen Bibliothek betrieben wurden, und die Gesamtzeit für den Klassifikator wurde durch Summieren der Dekompressions- und Bewertungszeiten für jedes Genom bestimmt. Wanduhrzeiten wurden für alle Klassifikatoren aufgezeichnet. Beim Vergleich von Kraken mit den anderen Klassifikatoren haben wir BLAST + 2.2.27, PhymmBL 4.0, NBC 1.1 und MetaPhlAn 1.7.6 verwendet. Klassifikatoren wurden alle auf demselben Computer mit 48 AMD Opteron 6172 2,1 GHz CPUs und 252 GB RAM mit Red Hat Enterprise Linux 5 ausgeführt. Die für die Geschwindigkeitsbewertung verwendeten Datensätze hatten jeweils 10.000 Lesevorgänge für alle anderen Programme als Kraken (und seine Varianten) und VLAN, die 10.000.000 Lesedatensätze verwendeten. Höhere Lesenummern wurden mit diesen schnelleren Programmen verwendet, um die Auswirkungen der Anfangs- und Endoperationen zu minimieren, die während der Ausführung der Programme stattfinden.

Obwohl Kraken das einzige der von uns untersuchten Programme ist, das explizit Operationen ausführt, um sicherzustellen, dass sich seine Daten vor der Klassifizierung im physischen Speicher befinden, wollten wir sicherstellen, dass alle Programme auf ähnliche Weise ausgewertet wurden. Bei der Auswertung der Geschwindigkeit lesen wir für jedes Programm alle Datenbankdateien (z. B. IMM-Dateien und BLAST-Datenbanken für PhymmBL, k-Mer-Frequenzlisten für NBC und den Bowtie-Index für VLAN) vor dem Ausführen des Programms dreimal in den Speicher, um den Datenbankinhalt im Cache des Betriebssystems (der im physischen Speicher gespeichert ist) zu speichern.

Reduzierte Datenbankgrößen

Um die 4-GB-Datenbank für unsere MiniKraken-Ergebnisse zu generieren, haben wir die ersten 18 jedes Blocks von 19 Datensätzen in der Standard-Kraken-Datenbank entfernt. Ein Schrumpffaktor von 19 wurde ausgewählt, da es der kleinste ganzzahlige Faktor war, der die Größe auf weniger als 4 GB reduzieren würde, eine Größe, die leicht in den Speicher vieler gängiger PCS passen kann. Für Benutzer, die mehr RAM zur Verfügung haben, ermöglicht Kraken die Verwendung eines kleineren Schrumpfungsfaktors, wodurch die Empfindlichkeit erhöht wird.

Verwendung von Entwurfsgenomen

Beim Aufbau der Kraken-GB-Datenbank stellten wir fest, dass an den Enden mehrere Contigs mit bekannten Adaptersequenzen vorhanden waren. In nachfolgenden Tests stellten wir auch fest, dass einige Sequenzen in Proben mit großen Mengen menschlicher Sequenz von dieser Datenbank durchweg falsch klassifiziert wurden, was zu dem Schluss führte, dass eine Kontamination wahrscheinlich in den Entwurfsgenomen vorhanden war. Um dieser Kontamination entgegenzuwirken, haben wir diese k-mer aus bekannten Adaptersequenzen sowie die ersten und letzten 20 k-mer aus jedem der Entwurfskontingente aus der Datenbank entfernt. Dies verbesserte zwar die Klassifizierung, beseitigte jedoch nicht das Problem der Fehlklassifizierung. Aus diesem Grund glauben wir, dass, wenn Entwurfsgenome in einer Kraken-Datenbank verwendet werden, sehr strenge Maßnahmen ergriffen werden sollten, um kontaminante Sequenzen aus der genomischen Bibliothek zu entfernen.

Clade-Ausschlussexperimente

Bei der erneuten Analyse des simBA-5-Datensatzes für unsere Clade-Ausschlussexperimente wurden einige Lesevorgänge für bestimmte Paare von gemessenen und ausgeschlossenen Rängen nicht verwendet. Wenn dem Ursprung eines Lesers ein taxonomischer Eintrag in einer der gemessenen oder ausgeschlossenen Ränge fehlte, Es wurde für dieses bestimmte Experiment nicht verwendet.

Darüber hinaus wurde in einem Experiment kein Read verwendet, es sei denn, mindestens zwei andere in unserer Datenbank vertretene Taxa (abgesehen von der ausgeschlossenen Klasse) auf dem ausgeschlossenen Rang teilten sich das Taxon der Ursprungsklasse auf dem gemessenen Rang. Zum Beispiel würde ein Lesewert aus der Gattung G nicht in einem Experiment verwendet, das die Genauigkeit am Klassenrang misst und den Gattungsrang ausschließt, es sei denn, die Heimatklasse von G hatte mindestens zwei andere Gattungen mit Genomen in Krakens genomischer Bibliothek. Ohne diesen Filterschritt würde eine Gattung ausgeschlossen, wenn sie die einzige Gattung in ihrer Klasse wäre, Kraken könnte unmöglich die richtige Klasse benennen, da alle Einträge in der Datenbank aus dieser Klasse ebenfalls ausgeschlossen würden. Dies ist der gleiche Ansatz in ähnlichen Experimenten, die zur Bewertung von PhymmBL verwendet wurden .

Human microbiome data classification

Wir klassifizierten die Daten des Human Microbiome Project mithilfe einer Kraken-Datenbank, die aus vollständigen RefSeq-Genomen von Bakterien, Archaeen und Viren zusammen mit dem Humangenom GRCh37 erstellt wurde. Wir haben die Sequenzen von drei Beitritten (SRS019120, SRS014468 und SRS015055) aus dem NCBI-Sequenzlesearchiv abgerufen, und für jeden Beitritt wurden zwei Läufe eingereicht. Alle Lesevorgänge wurden getrimmt, um Basen und Adaptersequenzen von geringer Qualität zu entfernen. Krona wurde verwendet, um alle taxonomischen Verteilungsdiagramme zu generieren.

Da die Sequenzen alle gepaarten Lesevorgänge waren, haben wir die Lesevorgänge zusammengefügt, indem wir die Übereinstimmungen mit einer Sequenz von ‚NNNNN‘ zwischen ihnen verkettet haben. Kraken ignoriert K-mers mit mehrdeutigen Nukleotiden, so dass die K-Mers, die diese ‚N‘ -Zeichen umfassen, die Klassifizierung nicht beeinflussen. Diese Operation ermöglichte es Kraken, ein Paar Lesevorgänge als eine Einheit zu klassifizieren, anstatt die Lesevorgänge separat klassifizieren zu müssen.

Software- und Datenverfügbarkeit