<p>I'm trying to figure out why. This will be a longer read, <strong>Tl;dr at the end</strong>:</p>

<h2>Doing a match</h2>

<p>There is a good match of <code>Q</code> against <code>LC074724</code> in <code>20170824</code>:</p>

<pre><code>$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
Q   LC074724.1  2581
</code></pre>

<p>Note that the bitscore of the <code>LC074724</code> match is <code>2581</code>.</p>

<p>But there is <strong>no match</strong> of <code>LC074724</code> <strong>at all in</strong> <code>20170821</code>:</p>

<pre><code>$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724
</code></pre>

<p>The option <code>-outfmt '6 qaccver saccver bitscore'</code> tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.</p>

<h2>What are the best hits for 20170821</h2>

<p>So what are the bitscores of the Q matches in 20170821?</p>

<pre><code>$ 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
</code></pre>

<p>The bitscore of <code>Q</code> against <code>LC074724</code> in <code>20170824</code> (2581) would have placed it <strong>in 5th position in the top 10 matches of Q</strong>, but <code>LC074724</code> is <strong>not matched at all</strong> in <code>20170821</code>! The lowest bitscores matched in <code>20170821</code> is 2354.</p>

<h2>What's the difference between the databases?</h2>

<p>Get the ids of the subjects in each database</p>

<pre><code>$ blastdbcmd -entry all -db 20170824 | egrep '^&gt;' | cut -c2- | sort &gt; 20170824.ids
$ blastdbcmd -entry all -db 20170821 | egrep '^&gt;' | cut -c2- | sort &gt; 20170821.ids

$ wc -l 2017082?.ids
    9754 20170821.ids
    9111 20170824.ids
   18865 total
</code></pre>

<p><code>20170821</code> has 643 sequences that are not in <code>20170824</code>:</p>

<pre><code>$ comm -23 20170821.ids 20170824.ids | wc -l
     643
</code></pre>

<p><code>20170824</code> has no sequences that are not in <code>20170821</code>:</p>

<pre><code>$ comm -13 20170821.ids 20170824.ids | wc -l
       0
</code></pre>

<p>and the two have 9111 sequences in common:</p>

<pre><code>$ comm -12 20170821.ids 20170824.ids | wc -l
    9111
</code></pre>

<p>These numbers match with the numbers from <code>wc -l</code> above.</p>

<p>So <code>20170821</code> is just <code>20170824</code> with an extra 643 sequences.</p>

<h2>Maybe LC074724 matches, but with a low bit score?</h2>

<p>BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:</p>

<pre><code>$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 -max_target_seqs 10000 &gt; 20170821-10000-matches.txt
</code></pre>

<p>And....</p>

<pre><code>$ 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
</code></pre>

<p>There it is! And <strong>it's in 5th position, just as it should be</strong> (based on the bitscores discussed above).</p>

<p>Just to confirm, taking away the <code>-max_target_seqs</code> option gets no match:</p>

<pre><code>$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
$
</code></pre>

<h2>Three confusingly named BLAST options</h2>

<p>I've always found the names and descriptions of the following three BLAST command-line options very confusing.</p>

<p>From the Formatting options section of blastn -help:</p>

<blockquote>
  <p>-num_descriptions =0>   Number of database sequences to show one-line descriptions for   Not applicable for outfmt > 4<br>
  Default = `500'    * Incompatible with:  max_target_seqs</p>
  
  <p>-num_alignments =0>   Number of database sequences to show alignments for   Default = `250'    * Incompatible with: 
  max_target_seqs</p>
</blockquote>

<p>And from the Restrict search or results section:</p>

<blockquote>
  <p>-max_target_seqs =1>   Maximum number of aligned sequences to keep   Not applicable for outfmt &lt;= 4   Default = `500'    *
  Incompatible with:  num_descriptions, num_alignments</p>
</blockquote>

<p>From the section names, it seems like only the latter option (<code>-max_target_seqs</code>) would have any effect on the search, and that the former two are <strong>just about what is displayed</strong>.</p>

<p>But... using either just <code>-num_descriptions</code> or <code>-num_alignments</code> with a big value also gets a match:</p>

<pre><code>$ 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
&gt;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
&gt;LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621
</code></pre>

<p>So BLAST is apparently finding the <code>LC074724</code> match even without the high <code>-max_target_seqs</code> option, but it's not displaying it unless there is a high <code>-num_alignments</code> or <code>-num_descriptions</code> value. (<strong>Keep in mind, this is the 5th best hit! Not somewhere >500</strong>)</p>

<p>Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:</p>

<pre><code>$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$
</code></pre>

<hr>

<p><strong>Tl;dr:</strong></p>

<p>In summary, the options </p>

<pre><code> -evalue
 -max_target_seqs
 -num_alignments
 -num_descriptions
</code></pre>

<p>are <strong>NOT</strong> just filters applied to the results, but instead directly influence the results <strong>found</strong>. </p>

<p>This means you can End up with a <strong>bacterial</strong> sequence as <strong>BEST hit</strong> when using a value <code>x</code>, and with no bacterial sequences <strong>at all</strong> in your top20 when using value <code>y</code>.</p>

<p>To quote the NCBI user service:</p>

<blockquote>
  <p>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.</p>
</blockquote>

<p>And</p>

<blockquote>
  <p>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.</p>
</blockquote>

<p>This was totally unexpected behaviour for me and is as far as I know is nowhere documented in the BLAST help or similar.</p>

CC BY-SA 3.0
2017-11-14T15:36:38.227Z
2017-11-14T15:45:06.537Z
2017-11-14T15:45:06.537Z
voiDnyx
<p>Disclaimer: This is a self answered question for documentation purpose and I adapted this from the following <a href="https://gist.github.com/sujaikumar/504b3b7024eaf3a04ef5" rel="noreferrer">github gist</a>. Especially from users <strong>terrycojones</strong> and <strong>peterjc</strong> as well as <strong>sujaikumar</strong> who raised the issue.</p>

<hr>

<p>I have a strange situation. I have two BLAST databases (<code>20170821</code> and <code>20170824</code>) that both contain a certain subject sequence <code>LC074724</code>. When I query them both with a given query, <code>Q</code>, in one case (<code>20170824</code>) there is a <strong>very strong match</strong> and in the other (<code>20170821</code>) <strong>no match at all</strong>.</p>

<p>Why is there not even a lower quality match found although the databases are nearly the same?</p>

<pre><code>BLAST version

$ blastn -version
blastn: 2.6.0+
 Package: blast 2.6.0, build May 11 2017 22:22:40
</code></pre>

<p>Both databases contain an identical subject sequence <code>LC074724</code>:</p>

<pre><code>$ blastdbcmd -entry LC074724 -db 20170824 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

$ blastdbcmd -entry LC074724 -db 20170821 | md5
fcfbc167c7f81dfd75aad1bf3db6220b
</code></pre>

<p>Only difference for the databases is <code>20170821</code> is just <code>20170824</code> with an extra 643 sequences. Rest is completely identical.</p>

<blast><phylogenetics><filtering>
Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?
# Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?

(Asked by: username_0 on 2017-11-14)

Disclaimer: This is a self answered question for documentation purpose and I adapted this from the following github gist. Especially from users terrycojones and peterjc as well as sujaikumar who raised the issue.

I have a strange situation. I have two BLAST databases (20170821 and 20170824) that both contain a certain subject sequence LC074724. When I query them both with a given query, Q, in one case (20170824) there is a very strong match and in the other (20170821) no match at all.
Why is there not even a lower quality match found although the databases are nearly the same?
BLAST version

$ blastn -version
blastn: 2.6.0+
 Package: blast 2.6.0, build May 11 2017 22:22:40

Both databases contain an identical subject sequence LC074724:
$ blastdbcmd -entry LC074724 -db 20170824 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

$ blastdbcmd -entry LC074724 -db 20170821 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

Only difference for the databases is 20170821 is just 20170824 with an extra 643 sequences. Rest is completely identical.

**Tags:** blast, phylogenetics, filtering

---

## 1 answers

### Answer by username_0 on 2017-11-14. Score: 11 | Accepted answer

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.


---

