Kraken: clasificación de secuencias metagenómicas ultrarrápidas utilizando alineaciones exactas

Algoritmo de clasificación de secuencias

Para clasificar una secuencia de ADN S, recopilamos todos los k-mers dentro de esa secuencia en un conjunto, denotado como K(S). A continuación, mapeamos cada k-mer en K (S), utilizando el algoritmo descrito a continuación, al taxón LCA de todos los genomas que contienen ese k-mer. Estos taxones de ACV y sus ancestros en el árbol de taxonomía forman lo que denominamos árbol de clasificación, un subárbol podado que se utiliza para clasificar S. Cada nodo en el árbol de clasificación se pondera con el número de k-mers en K(S) que se asignan al taxón asociado con ese nodo. Luego, cada ruta de raíz a hoja (RTL) en el árbol de clasificación se puntúa calculando la suma de todos los pesos de nodo a lo largo de la ruta. La ruta de puntuación máxima de RTL en el árbol de clasificación es la ruta de clasificación, y a S se le asigna la etiqueta correspondiente a su hoja (si hay múltiples rutas de puntuación máxima, se selecciona el ACV de todas las hojas de esas rutas). Este algoritmo, ilustrado en la Figura 1, permite a Kraken considerar cada k-mer dentro de una secuencia como una pieza de evidencia separada, y luego intentar resolver cualquier evidencia conflictiva si es necesario. Tenga en cuenta que para una elección adecuada de k, la mayoría de los k-mers se asignarán exclusivamente a una sola especie, simplificando en gran medida el proceso de clasificación. Las secuencias para las que no se encuentra ninguno de los k-mers en K(S) en ningún genoma quedan sin clasificar por este algoritmo.

El uso de la puntuación de ruta RTL en el árbol de clasificación es necesario a la luz de las diferencias inevitables entre las secuencias a clasificar y las secuencias presentes en cualquier biblioteca de genomas. Tales diferencias pueden, incluso para valores grandes de k, resultar en un k-mer que está presente en la biblioteca pero asociado con una especie muy alejada de la verdadera especie de origen. Al anotar las diversas rutas de RTL en el árbol de clasificación, podemos compensar estas diferencias y clasificar correctamente las secuencias incluso cuando una pequeña minoría de k-mers en una secuencia indica que a la secuencia se le debe asignar una etiqueta taxonómica incorrecta.

Creación de bases de datos

La implementación eficiente del algoritmo de clasificación de Kraken requiere que la asignación de k-mers a taxones se realice consultando una base de datos pre-calculada. Kraken crea esta base de datos a través de un proceso de varios pasos, comenzando con la selección de una biblioteca de secuencias genómicas. Kraken incluye una biblioteca predeterminada, basada en genomas microbianos completos en la base de datos RefSeq del Centro Nacional de Información Biotecnológica (NCBI), pero la biblioteca se puede personalizar según lo necesiten los usuarios individuales .

Una vez elegida la biblioteca, utilizamos el contador k-mer de múltiples hilos de Medusas para crear una base de datos que contenga cada 31-mer distinto de la biblioteca. Una vez que la base de datos está completa, los espacios de medusa de 4 bytes utilizados para almacenar los recuentos de k-mer en el archivo de la base de datos son utilizados por Kraken para almacenar los números de identificación taxonómica de los valores de LCA de k-mers. Una vez que la base de datos ha sido creada por Medusas, las secuencias genómicas de la biblioteca se procesan una a la vez. Para cada secuencia, el taxón asociado a ella se usa para establecer los valores de ACV almacenados de todos los k-mers en la secuencia. A medida que se procesan las secuencias, si un k-mer de una secuencia ha tenido su valor de ACV previamente establecido, entonces se calcula el ACV del valor almacenado y el taxón de la secuencia actual y ese ACV se almacena para el k-mer. La información de taxones se obtiene de la base de datos taxonómica del NCBI.

Estructura de base de datos y algoritmo de búsqueda

Debido a que Kraken utiliza con mucha frecuencia un k-mer como consulta de base de datos inmediatamente después de consultar un k-mer adyacente, y debido a que los k-mers adyacentes comparten una cantidad sustancial de secuencia, utilizamos el concepto de minimizador para agrupar k-mers similares. Para explicar nuestra aplicación de este concepto, aquí definimos la representación canónica de una secuencia de ADN S como la lexicográficamente más pequeña de S y el complemento inverso de S. Para determinar el minimizador de longitud M de un k-mer, consideramos la representación canónica de todos los M-mers en el k-mer, y seleccionamos el lexicográficamente más pequeño de esos M-mers como minimizador del k-mer. En la práctica, los k-mers adyacentes a menudo tendrán el mismo minimizador.

En la base de datos de Kraken, todos los k-mers con el mismo minimizador se almacenan de forma consecutiva, y se ordenan en orden lexicográfico de sus representaciones canónicas. Una consulta para un k-mer R se puede procesar buscando en un índice las posiciones en la base de datos donde se almacenarían los k-mers con el minimizador de R, y luego realizando una búsqueda binaria dentro de esa región (Figura 5). Debido a que los k-mers adyacentes a menudo tienen el mismo minimizador, el rango de búsqueda es a menudo el mismo entre dos consultas consecutivas, y la búsqueda en la primera consulta a menudo puede traer datos a la caché de la CPU que se utilizará en la segunda consulta. Al permitir los accesos a la memoria en consultas posteriores para acceder a los datos de la caché de la CPU en lugar de la RAM, esta estrategia hace que las consultas posteriores sean mucho más rápidas de lo que serían de otra manera. El índice que contiene las compensaciones de cada grupo de k-mers en la base de datos requiere 8 × 4M bytes. De forma predeterminada, Kraken utiliza minimizadores de 15 bp, pero el usuario puede modificar este valor; por ejemplo, al crear MiniKraken, utilizamos minimizadores de 13 bp para garantizar que el tamaño total de la base de datos se mantuviera por debajo de 4 GB.

Estructura de la base de datos de Kraken. Cada k-mer a consultar en la base de datos tiene una subcadena específica que es su minimizador. Para buscar un k-mer en la base de datos, se examinan las posiciones en la base de datos que contienen k-mers con el mismo minimizador. Estas posiciones se encuentran rápidamente examinando la matriz de desplazamiento del minimizador para las posiciones de inicio de los registros con el minimizador de k-mer (naranja) y el siguiente minimizador posible (azul). Dentro de un rango de registros asociados con un minimizador dado, los registros se ordenan por orden lexicográfico de sus k-mers, lo que permite completar una consulta mediante una búsqueda binaria en este rango.

Al implementar Kraken, realizamos optimizaciones adicionales a la estructura y al algoritmo de búsqueda descritos anteriormente. En primer lugar, como señalaron Roberts et al. , un simple ordenamiento lexicográfico de M-mers puede resultar en una distribución sesgada de minimizadores que sobre-representa M-mers de baja complejidad. En Kraken, tal sesgo crearía muchos rangos de búsqueda grandes, que requerirían más tiempo para buscar. Para crear una distribución más uniforme de minimizadores (y así acelerar las búsquedas), utilizamos la operación exclusiva-or (XOR) para alternar la mitad de los bits de la representación canónica de cada M-mer antes de comparar los M-mers entre sí utilizando el orden lexicográfico. Esta operación XOR codifica eficazmente el pedido estándar y evita el gran sesgo hacia minimizadores de baja complejidad.

También aprovechamos el hecho de que el rango de búsqueda es a menudo el mismo entre consultas para hacer que las consultas de Kraken sean más rápidas. En lugar de calcular el minimizador cada vez que realizamos una consulta, primero buscamos el rango anterior. Si nuestro k-mer consultado se encuentra en este rango, la consulta puede regresar inmediatamente. Si no se encuentra el k-mer, entonces se calcula el minimizador; si el minimizador de k-mer es el mismo que el último k-mer consultado, entonces la consulta falla, ya que se ha demostrado que el espacio de búsqueda del minimizador no tiene el k-mer. Solo si el minimizador ha cambiado, Kraken tiene que ajustar el rango de búsqueda y buscar de nuevo el k-mer.

Construcción de metagenomas simulados

Los metagenomas HiSeq y MiSeq se construyeron utilizando 20 juegos de lecturas de escopeta de genoma completo bacteriano. Estas lecturas se encontraron como parte del proyecto GAGE-B o en el Archivo de Lectura de Secuencias NCBI. Cada metagenoma contiene secuencias de diez genomas (Archivo adicional 1: Tabla S1). Tanto para los 10.000 como para los 10 millones de muestras leídas de cada uno de estos metagenomas, el 10% de sus secuencias se seleccionaron de cada uno de los diez conjuntos de datos del genoma (es decir, cada genoma tenía la misma abundancia de secuencias). Todas las secuencias se recortaron para eliminar bases de baja calidad y secuencias adaptadoras.

La composición de estos dos metagenomas plantea ciertos desafíos a nuestros clasificadores. Por ejemplo, Pelosinus fermentans, que se encuentra en nuestro metagenoma HiSeq, no puede ser identificado correctamente a nivel de género por Kraken (o cualquiera de los otros clasificadores descritos anteriormente), porque no hay genomas de Pelosinus en la base de datos de genomas completos de RefSeq; sin embargo, hay siete genomas de este tipo en la base de datos de Kraken-GB, incluidas seis cepas de P. fermentans. De manera similar, en nuestro metagenoma MiSeq, el Proteus vulgaris a menudo se clasifica incorrectamente a nivel de género porque el único genoma de Proteus en la base de datos de Kraken es un único genoma de Proteus mirabilis. Cinco genomas Proteus más están presentes en la base de datos de Kraken-GB, lo que permite a Kraken-GB clasificar mejor las lecturas de ese género. Además, el metagenoma MiSeq contiene cinco genomas de la familia Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus y Salmonella). La alta similitud de secuencias entre los géneros de esta familia puede dificultar la distinción entre géneros para cualquier clasificador.

El metagenoma simBA-5 fue creado simulando lecturas del conjunto completo de genomas bacterianos y arqueales en RefSeq. Los replicones de esos genomas se usaban si estaban asociados con un taxón que tenía una entrada asociada con el rango de género, dando como resultado un conjunto de replicones de 607 géneros. Luego usamos el simulador de lectura Mason con su modelo Illumina para producir 10 millones de lecturas de 100 bp de estos genomas. Primero creamos genomas simulados para cada especie, utilizando una tasa de SNP del 0,1% y una tasa de indel del 0,1% (ambos parámetros predeterminados), a partir de los cuales generamos las lecturas. Para las lecturas simuladas, multiplicamos las tasas de desajuste e indel predeterminados por cinco, lo que resulta en una tasa de desajuste promedio del 2% (que varía del 1% al principio de las lecturas al 6% al final) y una tasa de indel del 1% (probabilidad de inserción del 0,5% y probabilidad de eliminación del 0,5%). Para el metagenoma simBA-5, el conjunto de 10.000 lecturas se generó a partir de una muestra aleatoria del conjunto de 10 millones de lecturas.

Evaluación de precisión y velocidad

Elegimos medir la precisión principalmente a nivel de género, que era el nivel más bajo para el que podíamos determinar fácilmente la información taxonómica para las predicciones de PhymmBL y NBC de manera automatizada. (Esto se debe a la manera en que PhymmBL y NBC informan sus resultados). Debido a que algunos genomas no tienen entradas taxonómicas en los siete rangos (especie, género, familia, orden, clase, filo y reino), definimos la sensibilidad a nivel de género como A/B, donde A es el número de lecturas con un género asignado que se clasificaron correctamente en ese rango, y B es el número total de lecturas con cualquier género asignado. Definimos la sensibilidad de manera similar para otros rangos taxonómicos.

Debido a que Kraken puede clasificar una lectura a niveles superiores a la especie, medir su precisión requiere que definamos el efecto sobre la precisión de asignar el género correcto (por ejemplo) sin asignar una especie en absoluto. Por esta razón, definimos la precisión a nivel de rango como C/(D + E), donde C es el número de lecturas etiquetadas en o por debajo del taxón correcto en el rango medido, D es el número de lecturas etiquetadas en o por debajo del rango medido, y E es el número de lecturas etiquetadas incorrectamente por encima del rango medido. Por ejemplo, dada una lectura de R que debería etiquetarse como Escherichia coli, una etiqueta de R como E. coli, E. fergusonii o Escherichia mejoraría la precisión a nivel de género. Una etiqueta de Enterobacteriaceae (familia correcta) o Proteobacterias (filo correcto) no tendría ningún efecto en la precisión a nivel de género. Una etiqueta para R de Bacillus (género incorrecto) o Firmicutes (filo incorrecto) disminuiría la precisión a nivel de género.

Al evaluar la precisión de PhymmBL , siguiendo los consejos de sus desarrolladores, seleccionamos un umbral de confianza de género para nuestras comparaciones. Seleccionamos 3.333 lecturas del conjunto de datos simulado de complejidad media (simMC), que abarca 31 géneros diferentes. Para simular lecturas cortas de los datos de la secuencia de Sanger en el conjunto simMC, seleccionamos los últimos 100 pb de cada una de las lecturas. Luego comparamos PhymmBL con esas lecturas de 100 pb y evaluamos la sensibilidad a nivel de género y la precisión de las predicciones de PhymmBL con umbrales de confianza de género de 0 a 1, en incrementos de 0,05. Encontramos que un umbral de 0,65 produjo la puntuación F más alta (la media armónica de sensibilidad y precisión), con 0,60 y 0,70 también teniendo puntuaciones F dentro de 0.5 puntos porcentuales del máximo (Expediente adicional 1: Tabla S2). Por lo tanto, utilizamos el umbral de confianza de género de 0,65 en nuestras comparaciones. Aunque la selección de un umbral depende de las necesidades individuales de un usuario, y por lo tanto es hasta cierto punto arbitraria, un umbral seleccionado de esta manera proporciona una comparación más adecuada a un clasificador selectivo como Kraken que ningún umbral en absoluto.

Los resultados de tiempo y precisión al usar Megablast como clasificador se obtuvieron a partir de los datos de registro producidos por PhymmBL, ya que PhymmBL utiliza Megablast para su paso de alineación. Al asignar una etiqueta taxonómica a una lectura con Megablast, utilizamos el taxón asociado con la primera alineación reportada. Megablast se ejecutó con opciones predeterminadas.

La velocidad se evaluó utilizando la operación de un solo hilo de cada programa (excepto para NBC). PhymmBL fue alterado para que su llamada al programa blastn usara un hilo en lugar de dos. El NBC se ejecutó con 36 procesos concurrentes que operaban en conjuntos disjuntos de genomas en su biblioteca genómica, y el tiempo total para el clasificador se determinó sumando la descompresión y los tiempos de puntuación para cada genoma. Se registraron los tiempos de los relojes de pared de todos los clasificadores. Al comparar Kraken con los otros clasificadores, utilizamos BLAST + 2.2.27, PhymmBL 4.0, NBC 1.1 y MetaPhlAn 1.7.6. Todos los clasificadores se ejecutaban en el mismo ordenador, con 48 CPU AMD Opteron 6172 de 2,1 GHz y 252 GB de RAM, con Red Hat Enterprise Linux 5. Los conjuntos de datos utilizados para la evaluación de la velocidad tenían 10.000 lecturas cada uno para todos los programas que no fueran Kraken (y sus variantes) y MetaPhlAn, que usaban 10.000.000 de conjuntos de datos leídos. Se utilizaron números de lectura más altos con estos programas más rápidos para minimizar el efecto de las operaciones iniciales y finales que tienen lugar durante la ejecución de los programas.

Aunque Kraken es el único de los programas que examinamos que realiza operaciones explícitas para garantizar que sus datos estén en memoria física antes de la clasificación, queríamos asegurarnos de que todos los programas se evaluaran de manera similar. Al evaluar la velocidad, para cada programa, leemos todos los archivos de la base de datos (p.ej. Archivos IMM y bases de datos de voladura para PhymmBL, listas de frecuencias k-mer para NBC y el índice Bowtie para MetaPhlAn) en memoria tres veces antes de ejecutar el programa, con el fin de colocar el contenido de la base de datos en la caché del sistema operativo (que se almacena en la memoria física).

Tamaños de base de datos reducidos

Para generar la base de datos de 4 GB para nuestros resultados de MiniKraken, eliminamos los primeros 18 de cada bloque de 19 registros en la base de datos estándar de Kraken. Se seleccionó un factor de reducción de 19, ya que era el factor entero más pequeño que reduciría el tamaño a menos de 4 GB, un tamaño que puede caber fácilmente en la memoria de muchos ordenadores personales comunes. Para los usuarios que tienen más RAM disponible, Kraken permite utilizar un factor de contracción más pequeño, lo que aumentará la sensibilidad.

Uso de genomas borrador

Al construir la base de datos Kraken-GB, notamos que había varios contigs con secuencias de adaptadores conocidas en los extremos. En pruebas posteriores, también encontramos que algunas secuencias en muestras con grandes cantidades de secuencias humanas fueron sistemáticamente clasificadas erróneamente por esta base de datos, lo que nos llevó a concluir que la contaminación estaba probablemente presente en los genomas preliminares. En un intento de contrarrestar esta contaminación, eliminamos de la base de datos los k-mers de secuencias de adaptadores conocidas, así como los primeros y los últimos 20 k-mers de cada uno de los contiguos de borrador. Si bien esto mejoró la clasificación, no eliminó el problema de clasificación errónea. Por esta razón, creemos que si se usan genomas preliminares en una base de datos de Kraken, se deben usar medidas muy estrictas para eliminar las secuencias de contaminantes de la biblioteca genómica.

Experimentos de exclusión de clados

Al volver a analizar el conjunto de datos simBA-5 para nuestros experimentos de exclusión de clados, algunas lecturas no se utilizaron para ciertos pares de rangos medidos y excluidos. Si el origen de una lectura carecía de una entrada taxonómica en cualquiera de los rangos medidos o excluidos, no se usó para ese experimento en particular.

Además, una lectura no se utilizó en un experimento a menos que al menos otros dos taxones representados en nuestra base de datos (aparte del clado excluido) en el rango excluido compartieran el clado del taxón de origen en el rango medido. Por ejemplo, una lectura del género G no se usaría en un experimento que midiera la precisión en el rango de clase y excluyera el rango de género a menos que la clase casera de G tuviera al menos otros dos géneros con genomas en la biblioteca genómica de Kraken. Sin este paso de filtrado, si se excluyera un género cuando fuera el único género de su clase, Kraken no podría nombrar la clase correcta, ya que todas las entradas en la base de datos de esa clase también serían excluidas. Este es el mismo enfoque adoptado en experimentos similares que se utilizaron para evaluar PhymmBL .

Clasificación de datos de microbioma humano

Clasificamos los datos del Proyecto de Microbioma Humano utilizando una base de datos de Kraken hecha de genomas bacterianos, arqueales y virales RefSeq completos, junto con el genoma humano GRCh37. Recuperamos las secuencias de tres adhesiones (SRS019120, SRS014468 y SRS015055) del Archivo de Lectura de Secuencias de NCBI, y cada adhesión tuvo dos ejecuciones enviadas. Todas las lecturas se recortaron para eliminar bases de baja calidad y secuencias de adaptadores. La corona se utilizó para generar todas las parcelas de distribución taxonómica.

Debido a que todas las secuencias eran lecturas emparejadas, unimos las lecturas concatenando las parejas con una secuencia de ‘NNNNN’ entre ellas. Kraken ignora los k-mers con nucleótidos ambiguos, por lo que los k-mers que abarcan estos caracteres ‘N’ no afectan la clasificación. Esta operación permitió a Kraken clasificar un par de lecturas como una sola unidad en lugar de tener que clasificar a las parejas por separado.

Disponibilidad de software y datos