Estou tentando descobrir o porquê. Esta será uma leitura mais longa, Tl; dr no final :
Fazendo uma correspondência
Há uma boa correspondência de Q
contra LC074724
em 20170824
:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \ -task blastn -max_hsps 1 | grep LC074724Q LC074724.1 2581
Observe que o bitcore da correspondência LC074724
é 2581
.
Mas não há nenhuma correspondência de LC074724
em 20170821
:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \ -task blastn -max_hsps 1 | grep LC074724
A opção -outfmt '6 qaccver saccver bitscore'
diz ao BLAST para imprimir o número de acesso da consulta e versão, o acesso do assunto e versão, e o bitscore da partida.
Quais são os melhores resultados para 20170821
Então, quais são os bitscores das partidas Q em 20170821?
$ blastn - outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head3301261625902590256325542547254425402540
O bitcore de Q
contra LC074724
em 20170824
(2581) o teria colocado como na 5ª posição entre as 10 melhores correspondências de Q , mas LC074724
não é correspondido de forma alguma em 20170821
! O menor escore de bits correspondente em 20170821
é 2354.
Qual é a diferença entre os bancos de dados?
Obtenha os ids dos assuntos em cada banco de dados
$ blastdbcmd -entry all -db 20170824 | egrep '^ >' | cut -c2- | sort > 20170824.ids $ blastdbcmd -entry all -db 20170821 | egrep '^ >' | cut -c2- | sort > 20170821.ids $ wc -l 2017082? .ids 9754 20170821.ids 9111 20170824.ids 18865 total
20170821
tem 643 sequências que não estão em 20170824
:
$ comm -23 20170821.ids 20170824.ids | wc -l 643
20170824
não tem sequências que não estejam em 20170821
:
$ comm -13 20170821.ids 20170824.ids | wc -l 0
e os dois têm 9111 sequências em comum:
$ comm -12 20170821.ids 20170824.ids | wc -l 9111
Esses números correspondem aos números de wc -l
acima.
Portanto, 20170821
é apenas 20170824
com 643 sequências extras.
Talvez LC074724 corresponda, mas com uma pontuação de bits baixa?
O BLAST mostra apenas os 500 principais assuntos correspondidos por padrão. Vamos pedir para mostrar mais:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \ -task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000- Match.txt
E ....
$ head 20170821-10000-Match.txtQ Q 3301Q X234A 2616Q CS388973.1 2590Q DM059402.1 2590Q LC074724.1 2581Q KJ843188.1 2576Q AB697496.1 2576Q AB697499.1 2576Q AB697503.1 2576Q AB697508.1 2576
Pronto! E está na 5ª posição, como deveria estar (com base nos bitscores discutidos acima).
Apenas para confirmar, retirando o -max_target_seqs
opção não obtém correspondência:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \ -task blastn -max_hsps 1 | grep LC074724 $
Três opções BLAST com nomes confusos
Sempre achei os nomes e descrições das três opções de linha de comando do BLAST a seguir muito confusos.
Na seção Opções de formatação de blastn -help:
-num_descriptions = 0> Número de sequências de banco de dados para mostrar descrições de uma linha para Não aplicável para outfmt> 4
Padrão = `500 '* Incompatível com: max_target_seqs
-num_alignments = 0> Número de sequências de banco de dados para mostrar alinhamentos para Default = `250 '* Incompatível com: max_target_seqs
E da seção Restringir pesquisa ou resultados:
-max_target_seqs = 1> Número máximo de sequências alinhadas a manter Não aplicável para outfmt < = 4 Padrão = `500 '* Incompatível com: num_descriptions, num_alignments
De nomes de seção, parece que apenas a última opção ( -max_target_seqs
) teria algum efeito na pesquisa, e que as duas primeiras são apenas sobre o que é exibido .
Mas ... usar apenas -num_descriptions
ou -num_alignments
com um grande valor também obtém uma correspondência:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_descriptions 10000 | grep LC074724LC074724.1 DNA do vírus da hepatite B, genoma completo, isolado: p621 2581 0.0>LC074724.1 DNA do vírus da hepatite B, genoma completo, isolado: p621 $ blastn -db 20170821 -query Q.fasta -task blastn -max_hs00al | grep LC074724LC074724.1 DNA do vírus da hepatite B, genoma completo, isolado: p621 2581 0.0>LC074724.1 DNA do vírus da hepatite B, genoma completo, isolado: p621
Então, o BLAST está aparentemente encontrando o LC074724 corresponde mesmo sem a opção -max_target_seqs
alta, mas não está exibindo a menos que haja um -num_alignments
ou -num_descriptions code alto > valor. ( Tenha em mente que este é o quinto melhor acerto! Não em algum lugar> 500 )
Apenas para confirmar, quando nenhuma das 3 opções é fornecida e nenhuma opção -outfmt também é , a correspondência não foi encontrada:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724 $
Tl; dr:
Em resumo, as opções
-evalue -max_target_seqs -num_alignments -num_descriptions
são NÃO apenas filtros aplicados aos resultados, mas influenciam diretamente os resultados encontrados .
Isso significa que você pode terminar com uma sequência bacteriana como BEST hit ao usar um valor x
, e sem sequências bacterianas de todo em seu top20 ao usar o valor y
.
Para citar o serviço de usuário NCBI:
Não consideramos isso um bug, mas concordo que devemos documentar melhor essa possibilidade. Isso pode acontecer porque os limites, incluindo as sequências máximas de destino, são aplicados em uma fase inicial sem intervalo do algoritmo, bem como posteriormente.
E
O cortes são aplicados durante uma fase de extensão sem lacuna e alguns alinhamentos podem melhorar na pontuação durante a extensão com lacuna. A abordagem mais segura é relaxar (aumentar) o limite e lidar com um conjunto de resultados maior.
Esse foi um comportamento totalmente inesperado para mim e, até onde sei, não está documentado em nenhum lugar do BLAST ajuda ou semelhante.