I'm trying to figure out why. This will be a longer read, Tl;dr at the end:
Doing a match
There is a good match of Q against LC074724 in 20170824:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724
Q LC074724.1 2581
Note that the bitscore of the LC074724 match is 2581.
But there is no match of LC074724 at all in 20170821:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724
The option -outfmt '6 qaccver saccver bitscore' tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.
What are the best hits for 20170821
So what are the bitscores of the Q matches in 20170821?
$ blastn -outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head
3301
2616
2590
2590
2563
2554
2547
2544
2540
2540
The bitscore of Q against LC074724 in 20170824 (2581) would have placed it in 5th position in the top 10 matches of Q, but LC074724 is not matched at all in 20170821! The lowest bitscores matched in 20170821 is 2354.
What's the difference between the databases?
Get the ids of the subjects in each database
$ 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 has 643 sequences that are not in 20170824:
$ comm -23 20170821.ids 20170824.ids | wc -l
643
20170824 has no sequences that are not in 20170821:
$ comm -13 20170821.ids 20170824.ids | wc -l
0
and the two have 9111 sequences in common:
$ comm -12 20170821.ids 20170824.ids | wc -l
9111
These numbers match with the numbers from wc -l above.
So 20170821 is just 20170824 with an extra 643 sequences.
Maybe LC074724 matches, but with a low bit score?
BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000-matches.txt
And....
$ head 20170821-10000-matches.txt
Q Q 3301
Q X234A 2616
Q CS388973.1 2590
Q DM059402.1 2590
Q LC074724.1 2581
Q KJ843188.1 2576
Q AB697496.1 2576
Q AB697499.1 2576
Q AB697503.1 2576
Q AB697508.1 2576
There it is! And it's in 5th position, just as it should be (based on the bitscores discussed above).
Just to confirm, taking away the -max_target_seqs option gets no match:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724
$
Three confusingly named BLAST options
I've always found the names and descriptions of the following three BLAST command-line options very confusing.
From the Formatting options section of blastn -help:
-num_descriptions =0> Number of database sequences to show one-line descriptions for Not applicable for outfmt > 4
Default = `500' * Incompatible with: max_target_seqs
-num_alignments =0> Number of database sequences to show alignments for Default = `250' * Incompatible with:
max_target_seqs
And from the Restrict search or results section:
-max_target_seqs =1> Maximum number of aligned sequences to keep Not applicable for outfmt <= 4 Default = `500' *
Incompatible with: num_descriptions, num_alignments
From the section names, it seems like only the latter option (-max_target_seqs) would have any effect on the search, and that the former two are just about what is displayed.
But... using either just -num_descriptions or -num_alignments with a big value also gets a match:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_descriptions 10000 | grep LC074724
LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621 2581 0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_alignments 10000 | grep LC074724
LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621 2581 0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621
So BLAST is apparently finding the LC074724 match even without the high -max_target_seqs option, but it's not displaying it unless there is a high -num_alignments or -num_descriptions value. (Keep in mind, this is the 5th best hit! Not somewhere >500)
Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$
Tl;dr:
In summary, the options
-evalue
-max_target_seqs
-num_alignments
-num_descriptions
are NOT just filters applied to the results, but instead directly influence the results found.
This means you can End up with a bacterial sequence as BEST hit when using a value x, and with no bacterial sequences at all in your top20 when using value y.
To quote the NCBI user service:
We don't consider this a bug, but I agree that we should document this
possibility better. This can happen because limits, including max
target sequences, are applied in an early ungapped phase of the
algorithm, as well as later.
And
The cutoffs are applied during an ungapped extension phase, and some
alignments may improve in score during gapped extension. The safest
approach is to relax (increase) the threshold and deal with a larger
result set.
This was totally unexpected behaviour for me and is as far as I know is nowhere documented in the BLAST help or similar.
