Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-14T15:36:38.227Z
Answers.LastActivityDate: 2017-11-14T15:45:06.537Z
Answers.LastEditDate: 2017-11-14T15:45:06.537Z
Answers.UserDisplayName: voiDnyx
Body: <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>
  
Tags: <blast><phylogenetics><filtering>
Title: Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?
ThreadText: # 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](https://gist.github.com/sujaikumar/504b3b7024eaf3a04ef5). 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.
  
  
  ---
  
  