Kraken: ultra-rápida metagenomic sequência de classificação usando exata alinhamentos

Sequência de classificação algoritmo

Para classificar uma sequência de DNA S, recolhemos todos os k-mers dentro dessa sequência em um conjunto, denotado como K(S). Em seguida, mapeamos cada K-mer em K (S), usando o algoritmo descrito abaixo, para o táxon LCA de todos os genomas que contêm esse K-mer. Estes LCA táxons e seus ancestrais na árvore de taxonomia, formam o que chamamos de árvore de classificação, uma podada subárvore que é usado para classificar, S. Cada nó na árvore de classificação é ponderada com o número de k-mers em K(S) que mapeada para o táxon associada a esse nó. Em seguida, cada caminho raiz a folha (RTL) na árvore de classificação é pontuado calculando a soma de todos os pesos dos nós ao longo do caminho. O caminho RTL de pontuação máxima na árvore de classificação é o caminho de classificação e S recebe o rótulo correspondente à sua folha (se houver vários caminhos de pontuação máxima, o LCA de todas as folhas desses caminhos é selecionado). Este algoritmo, ilustrado na Figura 1, permite que Kraken considere cada K-mer dentro de uma sequência como uma peça separada de evidência e, em seguida, tente resolver qualquer evidência conflitante, se necessário. Observe que, para uma escolha apropriada de k, a maioria dos k-mers mapeará exclusivamente para uma única espécie, simplificando muito o processo de classificação. Sequências para as quais nenhum dos k-mers em K(S) são encontrados em qualquer genoma não são classificadas por este algoritmo.

o uso da pontuação do caminho RTL na árvore de classificação é necessário à luz das inevitáveis diferenças entre as sequências a serem classificadas e as sequências presentes em qualquer biblioteca de genomas. Tais diferenças podem, mesmo para grandes valores de k, resultar em um K-mer que está presente na biblioteca, mas associado a uma espécie distante da verdadeira espécie-fonte. Ao marcar os vários caminhos RTL na árvore de classificação, podemos compensar essas diferenças e classificar corretamente as sequências, mesmo quando uma pequena minoria de K-mers em uma sequência indica que a sequência deve receber um rótulo taxonômico incorreto.

criação de banco de dados

a implementação eficiente do algoritmo de classificação de Kraken requer que o mapeamento de K-mers para taxa seja realizado consultando um banco de dados pré-calculado. Kraken cria esse banco de dados por meio de um processo de várias etapas, começando com a seleção de uma biblioteca de sequências genômicas. Kraken inclui uma biblioteca padrão, baseada em genomas microbianos concluídos no banco de dados RefSeq do National Center for Biotechnology Information (NCBI), mas a biblioteca pode ser personalizada conforme necessário por usuários individuais .

uma vez que a Biblioteca é escolhida, usamos o contador K-mer multithread Medusa para criar um banco de dados contendo todos os 31-mer distintos na biblioteca. Assim que o banco de dados for concluído, os espaços de 4 bytes usados para armazenar as contagens k-mer no arquivo de banco de dados são usados pelo Kraken para armazenar os números de identificação taxonômica dos valores LCA do k-mers. Depois que o banco de dados foi criado pela Jellyfish, as sequências genômicas na biblioteca são processadas uma de cada vez. Para cada sequência, o táxon associado a ele é usado para definir os valores LCA armazenados de todos os K-mers na sequência. À medida que as sequências são processadas, se um k-mer de uma sequência tiver seu valor LCA definido anteriormente, o LCA do valor armazenado e o táxon da sequência atual são calculados e esse LCA é armazenado para o K-mer. As informações do táxon são obtidas no banco de dados de taxonomia do NCBI.

estrutura de Banco de dados e algoritmo de busca

Porque Kraken muito frequentemente usa um k-mer, uma consulta de banco de dados imediatamente após consultar um adjacentes a k-mer, e porque adjacentes a k-mers compartilhar uma quantidade substancial de sequência, nós utilizamos o minimizer conceito de grupo semelhante k-mers juntos. Para explicar nossa aplicação desse conceito, definimos aqui a representação canônica de uma sequência de DNA S como lexicograficamente menor de S e o complemento reverso de S. Para determinar o k-mer do minimizer de comprimento M, consideramos a representação canônica de todos os M-mers no k-mer, e selecione o lexicographically menor de todos os M-mers, como o k-mer do minimizer. Na prática, os K-mers adjacentes geralmente terão o mesmo minimizador.

no banco de dados de Kraken, todos os K-mers com o mesmo minimizador são armazenados consecutivamente e são classificados em ordem lexicográfica de suas representações canônicas. Uma consulta para um K-mer R pode então ser processada procurando em um índice as posições no banco de dados onde os k-mers com o minimizador de R seriam armazenados e, em seguida, realizando uma pesquisa binária nessa região (Figura 5). Como os K-mers adjacentes geralmente têm o mesmo minimizador, o intervalo de pesquisa costuma ser o mesmo entre duas consultas consecutivas, e a pesquisa na primeira consulta geralmente pode trazer dados para o cache da CPU que será usado na segunda consulta. Ao permitir que acessos de memória em consultas subsequentes acessem dados no cache da CPU em vez de RAM, essa estratégia torna as consultas subsequentes muito mais rápidas do que seriam de outra forma. O índice que contém as compensações de cada grupo de K-mers no banco de dados requer 8 × 4m bytes. Por padrão, o Kraken usa minimizadores de 15 bp, mas o usuário pode modificar esse valor; por exemplo, na criação do MiniKraken, usamos minimizadores de 13 bp para garantir que o tamanho total do banco de dados permanecesse abaixo de 4 GB.

Kraken estrutura de banco de dados. Cada K-mer a ser consultado contra o banco de dados tem uma substring específica que é seu minimizador. Para procurar um K-mer no banco de dados, as posições no banco de dados que contêm k-mers com o mesmo minimizador são examinadas. Essas posições são encontradas rapidamente examinando a matriz de deslocamento do minimizador para as posições iniciais dos registros com o minimizador k-mer (laranja) e o próximo minimizador possível (azul). Dentro de um intervalo de registros associados a um determinado minimizador, os registros são classificados por ordenação lexicográfica de seus K-mers, permitindo que uma consulta seja concluída usando uma pesquisa binária nesse intervalo.

na implementação do Kraken, fizemos outras otimizações para a estrutura e algoritmo de pesquisa descrito acima. Primeiro, conforme observado por Roberts et al. , uma simples ordenação lexicográfica de M-mers pode resultar em uma distribuição distorcida de minimizadores que representam excessivamente M-mers de baixa complexidade. Em Kraken, esse viés criaria muitos grandes intervalos de pesquisa, o que exigiria mais tempo para pesquisar. Para criar uma distribuição mais uniforme de minimizadores (e, assim, acelerar as pesquisas), usamos a operação exclusive-or (XOR) para alternar metade dos bits da representação canônica de cada M-mer antes de comparar os m-mers entre si usando ordenação lexicográfica. Esta operação XOR efetivamente embaralha a ordenação padrão e evita o grande viés em direção a minimizadores de baixa complexidade.

também aproveitamos o fato de que o intervalo de pesquisa costuma ser o mesmo entre as consultas para tornar as consultas de Kraken mais rápidas. Em vez de calcular o minimizador cada vez que realizamos uma consulta, primeiro pesquisamos o intervalo anterior. Se o nosso K-mer consultado for encontrado neste intervalo, a consulta pode retornar imediatamente. Se o k-mer não for encontrado, o minimizador será calculado; se o minimizador do K-mer for o mesmo que o último K-mer consultado, a consulta falhará, pois o espaço de pesquisa do minimizador foi mostrado para não ter o K-mer. Somente se o minimizador tiver mudado, o Kraken terá que ajustar o intervalo de pesquisa e pesquisar novamente pelo K-mer.

construindo metagenomas simulados

os metagenomas HiSeq e MiSeq foram construídos usando 20 conjuntos de leituras bacterianas de espingarda de genoma inteiro. Essas leituras foram encontradas como parte do projeto GAGE-B ou no arquivo de leitura de sequência do NCBI. Cada metagenoma contém sequências de dez genomas (arquivo adicional 1: Tabela S1). Para as 10.000 e 10 milhões de amostras lidas de cada um desses metagenomas, 10% de suas sequências foram selecionadas de cada um dos dez conjuntos de dados do genoma do componente (ou seja, cada genoma tinha abundância de sequência igual). Todas as sequências foram cortadas para remover bases de baixa qualidade e sequências de adaptadores.

a composição desses dois metagenomas apresenta certos desafios aos nossos classificadores. Por exemplo, Pelosinus fermentans, encontradas em nossa HiSeq metagenome, não pode ser corretamente identificados em nível de gênero nível de Kraken (ou qualquer uma das outras descritas anteriormente classificadores), porque não há Pelosinus genomas do RefSeq completa de genomas de banco de dados; no entanto, existem sete genomas em Kraken-GB do banco de dados, incluindo seis cepas de P. fermentans. Da mesma forma, em nosso Metagenoma MiSeq, Proteus vulgaris é frequentemente classificado incorretamente no nível do gênero porque o único genoma Proteus no banco de dados de Kraken é um único genoma Proteus mirabilis. Mais cinco genomas Proteus estão presentes no banco de dados do Kraken-GB, permitindo que o Kraken-GB classifique as leituras melhor desse gênero. Além disso, o Metagenoma MiSeq contém cinco genomas da família Enterobacteriaceae (Citrobacter, Enterobacter, Klebsiella, Proteus e Salmonella). A Alta semelhança de sequência entre os gêneros desta família pode dificultar a distinção entre gêneros para qualquer classificador.

o metagenoma simBA-5 foi criado simulando leituras do conjunto de genomas bacterianos e archaeal completos em RefSeq. Replicons desses genomas foram usados se estivessem associados a um táxon que tinha uma entrada associada ao gênero classificação, resultando em um conjunto de replicons de 607 gêneros. Em seguida, usamos o Mason read simulator com seu modelo Illumina para produzir 10 milhões de leituras de 100 bp desses genomas. Primeiro, criamos genomas simulados para cada espécie, usando uma taxa SNP de 0,1% e uma taxa indel de 0,1% (ambos os parâmetros padrão), a partir da qual geramos as leituras. Para as leituras simuladas, multiplicamos a incompatibilidade padrão e as taxas de indel por cinco, resultando em uma taxa média de incompatibilidade de 2% (variando de 1% no início das leituras a 6% nas extremidades) e uma taxa de indel de 1% (0,5% probabilidade de inserção e 0,5% probabilidade de exclusão). Para o metagenoma simBA-5, o conjunto de leitura 10,000 foi gerado a partir de uma amostra aleatória do conjunto de leitura 10 milhões.

avaliação de precisão e velocidade

optamos por medir a precisão principalmente no nível do gênero, que era o nível mais baixo para o qual poderíamos facilmente determinar as informações de taxonomia para phymmbl e as previsões da NBC de forma automatizada. (Isso se deve à maneira como PhymmBL e NBC relatam seus resultados). Porque alguns genomas não tem taxonômica entradas em todos os sete fileiras (espécie, gênero, família, ordem, classe, filo e reino), definimos gênero em nível de sensibilidade como A/B, onde a é o número de leituras com um determinado gênero que foram corretamente classificados em ordem, e B é o número total de leituras com qualquer atribuída gênero. Definimos sensibilidade de forma semelhante para outras classificações taxonômicas.Como Kraken pode classificar uma leitura em níveis acima da espécie, medir sua precisão requer que definamos o efeito na precisão de atribuir o gênero correto (por exemplo) sem atribuir uma espécie. Por essa razão, definimos a classificação do nível de precisão, como C/(D + E), onde C é o número de leituras rotulado na ou abaixo do correto táxon na medida classificação, D é o número de leituras rotulado na ou abaixo da medida de classificação, e e é o número de leituras incorretamente rotulado acima de medida de valor. Por exemplo, dado um R lido que deve ser rotulado como Escherichia coli, uma rotulagem de R como E. coli, E. fergusonii ou Escherichia melhoraria a precisão no nível do gênero. Um rótulo de Enterobacteriaceae (família correta) ou proteobactérias (filo correto) não teria efeito na precisão do nível do gênero. Um rótulo para R de Bacillus (gênero incorreto) ou Firmicutes (filo incorreto) diminuiria a precisão no nível do gênero.

ao avaliar a precisão do PhymmBL, seguindo o conselho de seus desenvolvedores , selecionamos um limite de confiança do gênero para nossas comparações. Selecionamos 3.333 leituras do conjunto de dados simulado de média complexidade (simMC), abrangendo 31 gêneros diferentes. Para simular leituras curtas dos dados da sequência Sanger no conjunto simMC, selecionamos os últimos 100 bp de cada uma das leituras. Em seguida, executamos o PhymmBL em relação a essas leituras de 100 bp e avaliamos a sensibilidade e a precisão do nível de gênero das previsões do PhymmBL com limiares de confiança do gênero de 0 a 1, em incrementos de 0,05. Descobrimos que um limiar de 0,65 produziu a maior pontuação F (A média harmônica de sensibilidade e precisão), com 0,60 e 0,70 também tendo pontuações F dentro de 0.5 pontos percentuais do máximo (arquivo adicional 1: Tabela S2). Portanto, usamos o limiar de confiança de 0,65 gênero em nossas comparações. Embora a seleção de um limite dependa das necessidades individuais de um usuário e, portanto, seja até certo ponto arbitrária, um limite selecionado dessa maneira fornece uma comparação mais adequada a um classificador seletivo, como Kraken, do que nenhum limite.

os resultados de tempo e precisão ao usar o Megablast como classificador foram obtidos a partir dos dados de log produzidos pelo PhymmBL, pois o PhymmBL usa o Megablast para sua etapa de alinhamento. Ao atribuir um rótulo taxonômico a uma leitura com Megablast, usamos o táxon associado ao primeiro alinhamento relatado. Megablast foi executado com opções padrão.

velocidade foi avaliada usando a operação single-threaded de cada programa (exceto para NBC). PhymmBL foi alterado para que sua chamada para o programa blastn usasse um thread em vez de dois. A NBC foi executada com 36 processos simultâneos operando em conjuntos disjuntos de genomas em sua biblioteca genômica, e o tempo total para o classificador foi determinado somando os tempos de descompressão e pontuação para cada genoma. Os horários dos relógios de parede foram registrados para todos os classificadores. Ao comparar Kraken com os outros classificadores, usamos BLAST+ 2.2.27, PhymmBL 4.0, NBC 1.1 e MetaPhlAn 1.7.6. Os classificadores eram todos executados no mesmo computador, com 48 CPUs AMD Opteron 6172 2.1 GHz e 252 GB de RAM, executando o Red Hat Enterprise Linux 5. Os conjuntos de dados usados para avaliação de velocidade tinham 10.000 leituras cada para todos os programas além de Kraken (e suas variantes) e MetaPhlAn, que usava 10.000.000 conjuntos de dados de leitura. Números de leitura mais altos foram usados com esses programas mais rápidos para minimizar o efeito das operações iniciais e finais que ocorrem durante a execução dos programas.

embora Kraken seja o único dos programas que examinamos que realiza explicitamente operações para garantir que seus dados estejam na memória física antes da classificação, queríamos ter certeza de que todos os programas foram avaliados de maneira semelhante. Ao avaliar a velocidade, para cada programa, lemos todos os arquivos de banco de dados (por exemplo. Arquivos IMM e bancos de dados BLAST para PhymmBL, listas de frequência k-mer para NBC e o índice Bowtie para MetaPhlAn) na memória três vezes antes de executar o programa, a fim de colocar o conteúdo do banco de dados no cache do sistema operacional (que é armazenado na memória física).

tamanhos de banco de dados reduzidos

para gerar o banco de dados de 4 GB para nossos resultados MiniKraken, removemos os primeiros 18 de cada bloco de 19 registros no banco de dados Kraken padrão. Um fator de encolhimento de 19 foi selecionado, pois era o menor fator inteiro que reduziria o tamanho para menos de 4 GB, um tamanho que pode caber facilmente na memória de muitos computadores pessoais comuns. Para usuários que têm mais RAM disponível, Kraken permite que um fator de encolhimento menor seja usado, o que dará maior sensibilidade.

uso de genomas de rascunho

ao construir o banco de dados Kraken-GB, notamos que havia vários contigs com sequências de adaptadores conhecidas nas extremidades. Em testes subsequentes, também descobrimos que algumas sequências em amostras com grandes quantidades de sequência humana foram consistentemente mal classificadas por este banco de dados, levando-nos a concluir que a contaminação provavelmente estava presente nos genomas de rascunho. Na tentativa de neutralizar essa contaminação, removemos do banco de dados os K-mers de sequências de adaptadores conhecidas, bem como os primeiros e últimos 20 K-mers de cada um dos contigs de rascunho. Embora isso tenha melhorado a classificação, não eliminou o problema de classificação incorreta. Por esse motivo, acreditamos que, se genomas de rascunho forem usados em um banco de dados Kraken, medidas muito rigorosas devem ser usadas para remover sequências de contaminantes da biblioteca genômica.

experimentos de exclusão de clados

ao re-analisar o conjunto de dados simBA-5 para nossos experimentos de exclusão de clados, algumas leituras não foram usadas para certos pares de classificações medidas e excluídas. Se a origem de uma leitura não tivesse uma entrada taxonômica em nenhuma das classificações medidas ou excluídas, ela não era usada para esse experimento específico.Além disso, uma leitura não foi usada em um experimento, a menos que pelo menos dois outros táxons representados em nosso banco de dados (além do clado excluído) na classificação excluída compartilhassem o táxon do clado de origem na classificação medida. Por exemplo, uma leitura do gênero G não seria usada em um experimento medindo a precisão na classificação da classe e excluindo a classificação do gênero, a menos que a classe doméstica de G tivesse pelo menos dois outros gêneros com genomas na biblioteca genômica de Kraken. Sem essa etapa de Filtragem, um gênero foi excluído quando era o único gênero em sua classe, Kraken não poderia nomear a classe correta, pois todas as entradas no banco de dados dessa classe também seriam excluídas. Esta é a mesma abordagem adotada em experimentos semelhantes que foram usados para avaliar o PhymmBL .

Classificação dos dados do microbioma humano

classificamos os dados do projeto microbioma humano usando um banco de dados Kraken feito de genomas bacterianos, archaeal e virais RefSeq completos, juntamente com o genoma humano GRCh37. Recuperamos as sequências de três acessos (SRS019120, SRS014468 e SRS015055) do arquivo de leitura de sequência do NCBI, e cada adesão teve duas execuções enviadas. Todas as leituras foram cortadas para remover bases de baixa qualidade e sequências de adaptadores. A coroa foi usada para gerar todas as parcelas de distribuição taxonômica.

como as sequências foram todas leituras emparelhadas, juntamos as leituras concatenando os companheiros com uma sequência de ‘NNNNN’ entre eles. Kraken ignora k-mers com nucleotídeos ambíguos, então os k-mers que abrangem esses caracteres ‘N’ não afetam a classificação. Esta operação permitiu que Kraken classificasse um par de leituras como uma única unidade, em vez de ter que classificar os companheiros separadamente.

disponibilidade de Software e dados