6
$\begingroup$

Very recently a BLAST parameter -max_target_seqs n got a lot of attention. 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 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?

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ relevant bioinformatics.stackexchange.com/questions/2846/… $\endgroup$
    – Chris_Rands
    Commented Sep 28, 2018 at 10:36
  • $\begingroup$ 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) $\endgroup$
    – llrs
    Commented Sep 28, 2018 at 12:26

2 Answers 2

Reset to default
3
$\begingroup$

The -max_target_seqs is also related to -evalue which does not work as one would think. Looking at the blast news 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. 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.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ For a manual filtering approach look at my answer. $\endgroup$
    – Kamil S Jaron
    Commented Feb 27, 2019 at 13:33
1
$\begingroup$

As a follow up to @Llopis' answer (no native blast solution), I wrote a script blast_filter.py and I added it into my collection of useful scripts 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).

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.