Answers.Body: <p>As a follow up to @Llopis' answer (no native blast solution), I wrote a script <a href="https://github.com/KamilSJaron/generic_genomics/blob/master/blast_filter.py" rel="nofollow noreferrer">blast_filter.py</a> and I added it into my collection of useful scripts <a href="https://github.com/KamilSJaron/generic_genomics" rel="nofollow noreferrer">generic_genomics</a>. The script filters out the blast output and it works with streams. To get the <code>n</code> (5) best blast hits by e-value you can run:</p>
  
  <pre><code>blastp -query <span class="math-container">$QUERY -db $</span>PROTEINS -evalue 1e-10 -outfmt 6 | blast_filter.py 5 &gt; $BLASTOUT
  </code></pre>
  
  <p>I assume that the blast output is sorted primarily by query name and secondarily by e-value. Also, note that default <code>max_target_seqs</code> in blast is 500, if there will be a query with 500 alignments, script will drop a warning (the query has no guarantee of the best hit found).</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-03T16:24:07.547Z
Answers.LastActivityDate: 2018-10-03T16:24:07.547Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <p>The <code>-max_target_seqs</code> is also related to <code>-evalue</code> which <a href="https://bioinformatics.stackexchange.com/questions/3409/blast-hits-disappearing-after-changing-evalue?rq=1">does not work</a> as one would think. Looking at the <a href="https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMDWeb&amp;PAGE_TYPE=BlastNews" rel="nofollow noreferrer">blast news</a> one finds that since version BLAST 2.2.27+ released, one should use the <code>-max_target_seqs</code>: </p>
  
  <blockquote>
    <p>4.) The output for reports without separate descriptions and alignments sections (all –outfmt greater than 4) should now use –max_target_seqs to control the output rather than –num_descriptions and –num_alignments.</p>
  </blockquote>
  
  <p>So it seems that <code>num_alignments</code> is also <a href="https://bioinformatics.stackexchange.com/questions/2846/why-does-a-very-strong-blast-hit-get-lost-when-i-change-num-alignments-num-desc">misunderstood</a>. I came to the conclusion that one cannot expect to get the best <code>n</code> alignments, unless you get all and filter yourself the best ones. </p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-27T13:33:04.760Z
Answers.Text: For a manual filtering approach look at my answer.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-28T10:32:45.987Z
Answers.LastActivityDate: 2018-09-28T10:45:31.567Z
Answers.LastEditDate: 2018-09-28T10:45:31.567Z
Answers.UserDisplayName: llrs
Body: <p>Very recently a <a href="https://www.ncbi.nlm.nih.gov/books/NBK279690/" rel="nofollow noreferrer">BLAST</a> parameter <code>-max_target_seqs n</code> got a <a href="https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/bty833/5106166" rel="nofollow noreferrer">lot of attention</a>. Instead of the intuitive interpretation (return the best <code>n</code> sequences) the parameters asks blast to return the first <code>n</code> sequences that pass the e-value threshold. This is affecting thousands of workflows and analyses that assumed the intuitive interpretation.</p>
  
  <p>I also found that there is another parameter called <code>-num_alignments</code>. What is the difference? Related question on <a href="https://www.biostars.org/p/120287/#340449" rel="nofollow noreferrer">BioStars</a> seems that is misinterpreting the <code>max_target_seqs</code> parameter. Is <code>num_alignments</code> also misunderstood?</p>
  
  <p>Now the main question. How to actually run blast to give me the best <code>n</code> hits in the database?</p>
  
  <p>Till now I thought that this example will do the job:</p>
  
  <pre><code>blastp -query <span class="math-container">$QUERY -db $</span>PROTEINS -out $BLASTOUT -evalue 1e-10 -outfmt 6 -num_alignments 5
  </code></pre>
  
  <p>-- edit --</p>
  
  <p>Now I understood that i can not achieve the desired output only by parameters of blast, but I still don't see an answer to how can I get the best n hits. So I wanted to check out how to get all the results and that's not intuitive either. The parameter <code>-max_target_seqs</code> is set by default to 500. Does it mean that if there are more than 500 significant hits, I don't have a guarantee for the best one? Do I have to specify <code>max_target_seqs</code> to some crazy high number to be sure that I got them all?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-28T10:36:41.797Z
Comments.Text: relevant https://bioinformatics.stackexchange.com/questions/2846/why-does-a-very-strong-blast-hit-get-lost-when-i-change-num-alignments-num-desc
Comments.UserDisplayName: Chris_Rands
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-09-28T12:26:29.957Z
Comments.Text: To your last edit: I would think so, without knowing the source code it might be needed. However, I think that the default value is the same as what they use for the web, so if you want to have the same results on your inline command and on the web it might better to leave it as is. (In case you wanted to modify the c++ code I tried to see the license of the software, but I couldn't find anything)
Comments.UserDisplayName: llrs
Tags: <blast>
Title: How to achieve blast results according to the intuitive interpretation of `-max_target_seqs`?
ThreadText: # How to achieve blast results according to the intuitive interpretation of `-max_target_seqs`?
  
  (Asked by: username_0 on 2018-09-28)
  
  Very recently a [BLAST](https://www.ncbi.nlm.nih.gov/books/NBK279690/) parameter `-max_target_seqs n` got a [lot of attention](https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/bty833/5106166). Instead of the intuitive interpretation (return the best `n` sequences) the parameters asks blast to return the first `n` sequences that pass the e-value threshold. This is affecting thousands of workflows and analyses that assumed the intuitive interpretation.
  
  I also found that there is another parameter called `-num_alignments`. What is the difference? Related question on [BioStars](https://www.biostars.org/p/120287/#340449) seems that is misinterpreting the `max_target_seqs` parameter. Is `num_alignments` also misunderstood?
  
  Now the main question. How to actually run blast to give me the best `n` hits in the database?
  
  Till now I thought that this example will do the job:
  
  ```
  blastp -query $QUERY -db $PROTEINS -out $BLASTOUT -evalue 1e-10 -outfmt 6 -num_alignments 5
  
  ```
  
  -- edit --
  
  Now I understood that i can not achieve the desired output only by parameters of blast, but I still don't see an answer to how can I get the best n hits. So I wanted to check out how to get all the results and that's not intuitive either. The parameter `-max_target_seqs` is set by default to 500. Does it mean that if there are more than 500 significant hits, I don't have a guarantee for the best one? Do I have to specify `max_target_seqs` to some crazy high number to be sure that I got them all?
  
  **Tags:** blast
  
  > **Comments:**
  > **username_1**: relevant https://bioinformatics.stackexchange.com/questions/2846/why-does-a-very-strong-blast-hit-get-lost-when-i-change-num-alignments-num-desc  
  > **username_2**: To your last edit: I would think so, without knowing the source code it might be needed. However, I think that the default value is the same as what they use for the web, so if you want to have the same results on your inline command and on the web it might better to leave it as is. (In case you wanted to modify the c++ code I tried to see the license of the software, but I couldn't find anything)  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2018-09-28. Score: 3 | Accepted answer
  
  The `-max_target_seqs` is also related to `-evalue` which [does not work](https://bioinformatics.stackexchange.com/questions/3409/blast-hits-disappearing-after-changing-evalue?rq=1) as one would think. Looking at the [blast news](https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMDWeb&PAGE_TYPE=BlastNews) one finds that since version BLAST 2.2.27+ released, one should use the `-max_target_seqs`:
  
  > 4.) The output for reports without separate descriptions and alignments sections (all –outfmt greater than 4) should now use –max\_target\_seqs to control the output rather than –num\_descriptions and –num\_alignments.
  
  So it seems that `num_alignments` is also [misunderstood](https://bioinformatics.stackexchange.com/questions/2846/why-does-a-very-strong-blast-hit-get-lost-when-i-change-num-alignments-num-desc). I came to the conclusion that one cannot expect to get the best `n` alignments, unless you get all and filter yourself the best ones.
  
  > **Comments:**
  > **username_0**: For a manual filtering approach look at my answer.  
  
  ---
  
  ### Answer by username_0 on 2018-10-03. Score: 1
  
  As a follow up to @Llopis' answer (no native blast solution), I wrote a script [blast\_filter.py](https://github.com/username_0/generic_genomics/blob/master/blast_filter.py) and I added it into my collection of useful scripts [generic\_genomics](https://github.com/username_0/generic_genomics). The script filters out the blast output and it works with streams. To get the `n` (5) best blast hits by e-value you can run:
  
  ```
  blastp -query $QUERY -db $PROTEINS -evalue 1e-10 -outfmt 6 | blast_filter.py 5 > $BLASTOUT
  
  ```
  
  I assume that the blast output is sorted primarily by query name and secondarily by e-value. Also, note that default `max_target_seqs` in blast is 500, if there will be a query with 500 alignments, script will drop a warning (the query has no guarantee of the best hit found).
  
  
  ---
  
  