Kraken: classification de séquence métagénomique ultrarapide utilisant des alignements exacts

Algorithme de classification de séquence

Pour classer une séquence d’ADN S, nous collectons tous les k-mers de cette séquence en un ensemble, noté K(S). Nous mappons ensuite chaque k-mer en K(S), en utilisant l’algorithme décrit ci-dessous, au taxon LCA de tous les génomes qui contiennent ce k-mer. Ces taxons de LCA et leurs ancêtres dans l’arbre de taxonomie forment ce que nous appelons l’arbre de classification, un sous-arbre taillé qui est utilisé pour classer S. Chaque nœud de l’arbre de classification est pondéré avec le nombre de k-mers en K(S) qui ont été mappés au taxon associé à ce nœud. Ensuite, chaque chemin de la racine à la feuille (RTL) dans l’arbre de classification est noté en calculant la somme de tous les poids de nœuds le long du chemin. Le chemin RTL de notation maximale dans l’arbre de classification est le chemin de classification, et S se voit attribuer l’étiquette correspondant à sa feuille (s’il y a plusieurs chemins de notation maximale, l’ACV de toutes les feuilles de ces chemins est sélectionnée). Cet algorithme, illustré à la figure 1, permet à Kraken de considérer chaque k-mer dans une séquence comme une preuve distincte, puis de tenter de résoudre toute preuve contradictoire si nécessaire. Notez que pour un choix approprié de k, la plupart des k-mers seront cartographiées uniquement à une seule espèce, ce qui simplifiera grandement le processus de classification. Les séquences pour lesquelles aucun des k-mers de K(S) ne se trouve dans un génome quelconque ne sont pas classées par cet algorithme.

L’utilisation de la notation de chemin RTL dans l’arbre de classification est nécessaire au vu des différences inévitables entre les séquences à classer et les séquences présentes dans toute bibliothèque de génomes. De telles différences peuvent, même pour de grandes valeurs de k, aboutir à un k-mer présent dans la bibliothèque mais associé à une espèce très éloignée de l’espèce source véritable. En marquant les différents chemins RTL dans l’arbre de classification, nous pouvons compenser ces différences et classer correctement les séquences même lorsqu’une petite minorité de k-mers dans une séquence indique qu’une étiquette taxonomique incorrecte devrait être attribuée à la séquence.

Création de base de données

La mise en œuvre efficace de l’algorithme de classification de Kraken nécessite que le mappage des k-mers aux taxons soit effectué en interrogeant une base de données pré-calculée. Kraken crée cette base de données à travers un processus en plusieurs étapes, en commençant par la sélection d’une bibliothèque de séquences génomiques. Kraken comprend une bibliothèque par défaut, basée sur des génomes microbiens complétés dans la base de données RefSeq du National Center for Biotechnology Information (NCBI), mais la bibliothèque peut être personnalisée selon les besoins des utilisateurs individuels.

Une fois la bibliothèque choisie, nous utilisons le compteur k-mer multithread de Jellyfish pour créer une base de données contenant chaque 31-mer distinct de la bibliothèque. Une fois la base de données terminée, les espaces de 4 octets utilisés par les méduses pour stocker les nombres de k-mer dans le fichier de base de données sont à la place utilisés par Kraken pour stocker les numéros d’identification taxonomiques des valeurs de l’ACV des k-mers. Une fois la base de données créée par les méduses, les séquences génomiques de la bibliothèque sont traitées une à la fois. Pour chaque séquence, le taxon qui lui est associé est utilisé pour définir les valeurs de l’ACV stockées de tous les k-mers de la séquence. Au fur et à mesure que les séquences sont traitées, si une k-mer d’une séquence a eu sa valeur d’ACV précédemment définie, alors l’ACV de la valeur stockée et du taxon de la séquence actuelle est calculée et cette ACV est stockée pour le k-mer. Les informations sur les taxons sont obtenues à partir de la base de données de taxonomie NCBI.

Structure de base de données et algorithme de recherche

Parce que Kraken utilise très fréquemment un k-mer comme requête de base de données immédiatement après avoir interrogé un k-mer adjacent, et parce que les k-mers adjacentes partagent une quantité substantielle de séquence, nous utilisons le concept de minimiseur pour regrouper des k-mers similaires. Pour expliquer notre application de ce concept, nous définissons ici la représentation canonique d’une séquence d’ADN S comme le plus petit lexicographiquement de S et le complément inverse de S. Pour déterminer le minimiseur de longueur M d’un k-mer, nous considérons la représentation canonique de tous les M-mers dans le k-mer, et sélectionnons le plus petit lexicographiquement de ces M-mers comme minimiseur du k-mer. En pratique, les k-mers adjacentes auront souvent le même minimiseur.

Dans la base de données de Kraken, tous les k-mers avec le même minimiseur sont stockés consécutivement et sont triés dans l’ordre lexicographique de leurs représentations canoniques. Une requête pour un k-mer R peut ensuite être traitée en recherchant dans un index les positions dans la base de données où les k-mers avec le minimiseur de R seraient stockés, puis en effectuant une recherche binaire dans cette région (Figure 5). Étant donné que les k-mers adjacentes ont souvent le même minimiseur, la plage de recherche est souvent la même entre deux requêtes consécutives, et la recherche dans la première requête peut souvent apporter des données dans le cache du processeur qui sera utilisé dans la deuxième requête. En permettant aux accès à la mémoire dans les requêtes ultérieures d’accéder aux données du cache du processeur au lieu de la RAM, cette stratégie rend les requêtes ultérieures beaucoup plus rapides qu’elles ne le seraient autrement. L’index contenant les décalages de chaque groupe de k-mers dans la base de données nécessite 8 × 4 millions d’octets. Par défaut, Kraken utilise des minimiseurs de 15 pb, mais l’utilisateur peut modifier cette valeur ; par exemple, lors de la création de MiniKraken, nous avons utilisé des minimiseurs de 13 pb pour garantir que la taille totale de la base de données reste inférieure à 4 Go.

Structure de la base de données Kraken. Chaque k-mer à interroger par rapport à la base de données a une sous-chaîne spécifique qui est son minimiseur. Pour rechercher un k-mer dans la base de données, les positions dans la base de données qui contiennent des k-mers avec le même minimiseur sont examinées. Ces positions sont rapidement trouvées en examinant le tableau de décalage du minimiseur pour les positions de début des enregistrements avec le minimiseur du k-mer (orange) et le minimiseur suivant possible (bleu). Dans une plage d’enregistrements associés à un minimiseur donné, les enregistrements sont triés par ordre lexicographique de leurs k-mers, ce qui permet de compléter une requête en utilisant une recherche binaire sur cette plage.

Lors de la mise en œuvre de Kraken, nous avons effectué d’autres optimisations de la structure et de l’algorithme de recherche décrits ci-dessus. Premièrement, comme l’ont noté Roberts et coll. , un simple ordre lexicographique des M-mers peut entraîner une distribution biaisée des minimiseurs qui sur-représente les M-mers de faible complexité. Dans Kraken, un tel biais créerait de nombreuses grandes plages de recherche, ce qui nécessiterait plus de temps pour rechercher. Pour créer une distribution plus uniforme des minimiseurs (et ainsi accélérer les recherches), nous utilisons l’opération exclusive-or (XOR) pour basculer la moitié des bits de la représentation canonique de chaque M-mer avant de comparer les M-mers les uns aux autres en utilisant l’ordre lexicographique. Cette opération XOR brouille efficacement la commande standard et empêche le biais important vers des minimiseurs à faible complexité.

Nous profitons également du fait que la plage de recherche est souvent la même entre les requêtes pour accélérer les requêtes de Kraken. Plutôt que de calculer le minimiseur chaque fois que nous effectuons une requête, nous recherchons d’abord la plage précédente. Si notre k-mer interrogé se trouve dans cette plage, la requête peut revenir immédiatement. Si le k-mer n’est pas trouvé, le minimiseur est calculé; si le minimiseur du k-mer est le même que le dernier k-mer interrogé, la requête échoue, car l’espace de recherche du minimiseur n’a pas le k-mer. Ce n’est que si le minimiseur a changé que Kraken doit ajuster la plage de recherche et rechercher à nouveau le k-mer.

Construction de métagénomes simulés

Les métagénomes HiSeq et MiSeq ont été construits à l’aide de 20 séries de lectures bactériennes à génome entier. Ces lectures ont été trouvées soit dans le cadre du projet GAGE-B, soit dans l’archive de lecture de séquences NCBI. Chaque métagénome contient des séquences de dix génomes (fichier supplémentaire 1 : Tableau S1). Pour les 10 000 et 10 millions d’échantillons lus de chacun de ces métagénomes, 10% de leurs séquences ont été sélectionnées dans chacun des ensembles de données génomiques à dix composants (c’est-à-dire que chaque génome avait une abondance de séquences égale). Toutes les séquences ont été coupées pour éliminer les bases et les séquences d’adaptateur de mauvaise qualité.

La composition de ces deux métagénomes pose certains défis à nos classificateurs. Par exemple, Pelosinus fermentans, trouvé dans notre métagénome HiSeq, ne peut pas être correctement identifié au niveau du genre par Kraken (ou l’un des autres classificateurs précédemment décrits), car il n’y a pas de génomes de Pelosinus dans la base de données complète sur les génomes de RefSeq; cependant, il y a sept génomes de ce type dans la base de données de Kraken-GB, y compris six souches de P. fermentans. De même, dans notre métagénome MiSeq, Proteus vulgaris est souvent classé incorrectement au niveau du genre car le seul génome de Proteus dans la base de données de Kraken est un génome unique de Proteus mirabilis. Cinq autres génomes de Proteus sont présents dans la base de données de Kraken-GB, ce qui permet à Kraken-GB de mieux classer les lectures de ce genre. De plus, le métagénome MiSeq contient cinq génomes de la famille des Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus et Salmonella). La forte similitude de séquence entre les genres de cette famille peut rendre difficile la distinction entre les genres pour tout classificateur.

Le métagénome simBA-5 a été créé en simulant des lectures de l’ensemble des génomes bactériens et archaïques complets dans RefSeq. Les réplicons de ces génomes ont été utilisés s’ils étaient associés à un taxon dont l’entrée était associée au rang de genre, ce qui a donné un ensemble de réplicons de 607 genres. Nous avons ensuite utilisé le simulateur de lecture Mason avec son modèle Illumina pour produire 10 millions de lectures de 100 pb à partir de ces génomes. Nous avons d’abord créé des génomes simulés pour chaque espèce, en utilisant un taux de SNP de 0,1% et un taux d’indel de 0,1% (les deux paramètres par défaut), à partir desquels nous avons généré les lectures. Pour les lectures simulées, nous avons multiplié par cinq les taux d’inadéquation par défaut et d’indel, ce qui donne un taux d’inadéquation moyen de 2% (allant de 1% au début des lectures à 6% aux extrémités) et un taux d’indel de 1% (probabilité d’insertion de 0,5% et probabilité de suppression de 0,5%). Pour le métagénome simBA-5, l’ensemble de 10 000 lectures a été généré à partir d’un échantillon aléatoire de l’ensemble de 10 millions de lectures.

Évaluation de la précision et de la vitesse

Nous avons choisi de mesurer la précision principalement au niveau du genre, qui était le niveau le plus bas pour lequel nous pouvions facilement déterminer les informations de taxonomie pour les prédictions de PhymmBL et NBC de manière automatisée. (Cela est dû à la manière dont PhymmBL et NBC rapportent leurs résultats). Parce que certains génomes n’ont pas d’entrées taxonomiques aux sept rangs (espèce, genre, famille, ordre, classe, phylum et royaume), nous avons défini la sensibilité au niveau du genre comme A / B, où A est le nombre de lectures avec un genre attribué qui ont été correctement classées à ce rang, et B est le nombre total de lectures avec un genre attribué. Nous avons défini la sensibilité de la même manière pour d’autres rangs taxonomiques.

Comme Kraken peut classer une lecture à des niveaux supérieurs à l’espèce, la mesure de sa précision nous oblige à définir l’effet sur la précision de l’attribution du genre correct (par exemple) tout en n’attribuant aucune espèce. Pour cette raison, nous avons défini la précision au niveau du rang comme C / (D + E), où C est le nombre de lectures étiquetées à ou en dessous du taxon correct au rang mesuré, D est le nombre de lectures étiquetées à ou en dessous du rang mesuré, et E est le nombre de lectures incorrectement étiquetées au-dessus du rang mesuré. Par exemple, étant donné une lecture de R qui devrait être étiquetée comme Escherichia coli, un étiquetage de R comme E. coli, E. fergusonii ou Escherichia améliorerait la précision au niveau du genre. Une étiquette d’Enterobacteriaceae (famille correcte) ou de Protéobactéries (phylum correct) n’aurait aucun effet sur la précision au niveau du genre. Une étiquette pour R de Bacillus (genre incorrect) ou de Firmicutes (phylum incorrect) diminuerait la précision au niveau du genre.

Lors de l’évaluation de la précision de PhymmBL, suivant les conseils de ses développeurs, nous avons sélectionné un seuil de confiance en genre pour nos comparaisons. Nous avons sélectionné 3 333 lectures de l’ensemble de données simMC (simmed medium complexity), couvrant 31 genres différents. Pour simuler des lectures courtes à partir des données de séquence de Sanger dans l’ensemble simMC, nous avons sélectionné les 100 derniers pb de chacune des lectures. Nous avons ensuite comparé PhymmBL à ces lectures de 100 pb et évalué la sensibilité au niveau du genre et la précision des prédictions de PhymmBL avec des seuils de confiance du genre de 0 à 1, par incréments de 0,05. Nous avons constaté qu’un seuil de 0,65 donnait le score F le plus élevé (la moyenne harmonique de sensibilité et de précision), 0,60 et 0,70 ayant également des scores F à l’intérieur de 0.5 points de pourcentage du maximum (fichier supplémentaire 1 : Tableau S2). Nous avons donc utilisé le seuil de confiance de 0,65 pour le genre dans nos comparaisons. Bien que la sélection d’un seuil dépende des besoins individuels d’un utilisateur et soit donc dans une certaine mesure arbitraire, un seuil ainsi sélectionné fournit une comparaison plus appropriée avec un classificateur sélectif tel que Kraken qu’aucun seuil du tout.

Les résultats de temps et de précision lors de l’utilisation de Megablast comme classificateur ont été obtenus à partir des données de journal produites par PhymmBL, car PhymmBL utilise Megablast pour son étape d’alignement. Lors de l’attribution d’une étiquette taxonomique à une lecture avec Megablast, nous avons utilisé le taxon associé au premier alignement signalé. Megablast a été exécuté avec des options par défaut.La vitesse

a été évaluée en utilisant le fonctionnement mono-thread de chaque programme (sauf pour NBC). PhymmBL a été modifié de sorte que son appel au programme blastn utilisait un thread au lieu de deux. NBC a été exécuté avec 36 processus simultanés opérant sur des ensembles disjoints de génomes dans sa bibliothèque génomique, et le temps total pour le classificateur a été déterminé en additionnant les temps de décompression et de notation pour chaque génome. Les heures d’horloge murale ont été enregistrées pour tous les classificateurs. En comparant Kraken aux autres classificateurs, nous avons utilisé BLAST+2.2.27, PhymmBL 4.0, NBC 1.1 et MetaPhlAn 1.7.6. Les classificateurs étaient tous exécutés sur le même ordinateur, avec 48 processeurs AMD Opteron 6172 2,1 GHz et 252 Go de RAM, exécutant Red Hat Enterprise Linux 5. Les ensembles de données utilisés pour l’évaluation de la vitesse avaient chacun 10 000 lectures pour tous les programmes autres que Kraken (et ses variantes) et MetaPhlAn, qui utilisaient 10 000 000 ensembles de données en lecture. Des nombres de lecture plus élevés ont été utilisés avec ces programmes plus rapides pour minimiser l’effet des opérations initiales et finales qui ont lieu pendant l’exécution des programmes.

Bien que Kraken soit le seul des programmes que nous avons examinés à effectuer explicitement des opérations pour s’assurer que ses données sont en mémoire physique avant la classification, nous voulions être sûrs que tous les programmes étaient évalués de la même manière. Lors de l’évaluation de la vitesse, pour chaque programme, nous lisons tous les fichiers de base de données (par ex. Fichiers IMM et bases de données BLAST pour PhymmBL, listes de fréquences k-mer pour NBC et index Bowtie pour MetaPhlAn) en mémoire trois fois avant d’exécuter le programme, afin de placer le contenu de la base de données dans le cache du système d’exploitation (qui est stocké dans la mémoire physique).

Tailles de base de données réduites

Pour générer la base de données de 4 Go pour nos résultats MiniKraken, nous avons supprimé les 18 premiers enregistrements de chaque bloc de 19 dans la base de données Kraken standard. Un facteur de réduction de 19 a été sélectionné car il s’agissait du plus petit facteur entier qui réduirait la taille à moins de 4 Go, une taille qui peut facilement s’intégrer dans la mémoire de nombreux ordinateurs personnels courants. Pour les utilisateurs qui ont plus de RAM disponible, Kraken permet d’utiliser un facteur de rétrécissement plus petit, ce qui donnera une sensibilité accrue.

Utilisation de génomes d’ébauche

Lors de la construction de la base de données Kraken-GB, nous avons remarqué qu’il y avait plusieurs contigs avec des séquences d’adaptateurs connues aux extrémités. Lors de tests ultérieurs, nous avons également constaté que certaines séquences dans des échantillons contenant de grandes quantités de séquences humaines étaient systématiquement mal classées par cette base de données, ce qui nous a amenés à conclure que la contamination était probablement présente dans les génomes provisoires. Pour tenter de contrer cette contamination, nous avons retiré de la base de données les k-mers des séquences d’adaptateurs connues, ainsi que les 20 premiers et derniers k-mers de chacun des contigs de brouillon. Bien que cela ait amélioré la classification, cela n’a pas éliminé le problème de la mauvaise classification. Pour cette raison, nous pensons que si des génomes brouillons sont utilisés dans une base de données Kraken, des mesures très strictes devraient être utilisées pour éliminer les séquences de contaminants de la banque génomique.

Expériences d’exclusion de clade

Lors de la réanalyse de l’ensemble de données simBA-5 pour nos expériences d’exclusion de clade, certaines lectures n’ont pas été utilisées pour certaines paires de rangs mesurés et exclus. Si l’origine d’une lecture n’avait pas d’entrée taxonomique à l’un des rangs mesurés ou exclus, elle n’a pas été utilisée pour cette expérience particulière.

De plus, une lecture n’a pas été utilisée dans une expérience à moins qu’au moins deux autres taxons représentés dans notre base de données (à l’exception du clade exclu) au rang exclu ne partagent le taxon du clade d’origine au rang mesuré. Par exemple, une lecture du genre G ne serait pas utilisée dans une expérience mesurant la précision au rang de classe et excluant le rang de genre à moins que la classe d’origine de G ait au moins deux autres genres avec des génomes dans la bibliothèque génomique de Kraken. Sans cette étape de filtrage, si un genre était exclu alors qu’il était le seul genre de sa classe, Kraken ne pouvait pas nommer la classe correcte, car toutes les entrées de la base de données de cette classe seraient également exclues. C’est la même approche adoptée dans des expériences similaires qui ont été utilisées pour évaluer PhymmBL.

Classification des données sur le microbiome humain

Nous avons classé les données du Projet sur le microbiome humain à l’aide d’une base de données Kraken composée de génomes bactériens, archaïques et viraux RefSeq complets, ainsi que du génome humain GRCh37. Nous avons récupéré les séquences de trois accessions (SRS019120, SRS014468 et SRS015055) de l’archive de lecture de séquences NCBI, et chaque adhésion avait deux séries soumises. Toutes les lectures ont été coupées pour supprimer les bases et les séquences d’adaptateur de mauvaise qualité. Krona a été utilisé pour générer toutes les parcelles de répartition taxonomique.

Parce que les séquences étaient toutes des lectures appariées, nous avons joint les lectures ensemble en concaténant les compagnons avec une séquence de ‘NNNNN’ entre elles. Kraken ignore les k-mers avec des nucléotides ambigus, de sorte que les k-mers qui couvrent ces caractères « N » n’affectent pas la classification. Cette opération a permis à Kraken de classer une paire de lectures en une seule unité plutôt que d’avoir à classer les compagnons séparément.

Disponibilité des logiciels et des données