What BLAST's max-target-sequences doesn't do

This is a short post to highlight a scary BLAST+ -max_target_seqs bug found and reported by Sujai Kumar, which he discovered in the course of working on some puzzling Blobtools output while analysing the tardigrade genome.

In essence we, and probably most BLAST+ users, had assumed the -max_target_seqs command line option was applied after the search was finished prior to output, and likewise for the plain text output's -num_descriptions and -num_alignments (see also #5 on my Christmas 2014 BLAST wish list).

$ blastp -help
...
-max_target_seqs <Integer, >=1>
Maximum number of aligned sequences to keep
Not applicable for outfmt <= 4
Default = `500'
* Incompatible with: num_descriptions, num_alignments
...

Then Sujai found a reproducible counter example contradicting this mental model, where the top 20 hits switched from all being Eukaryotes to all Bacteria depending on -max_target_seqs.

The NCBI reply by email was:

Hello,
Thank you for the report. 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. In some cases a final HSP will improve enough in the later gapped phase to rise to the top hits. In your case, relaxing the limit to 200 appears to have allowed hits that would have been excluded in the ungapped phase at 100 max target sequences to rise.

So basically while most of the time -max_target_seqs does seem to be a final filter before the output is written, this setting actually happens much earlier in the search and can (unexpectedly) cull what turns out to be the best hit.

I agree with the NCBI BLAST team that this needs better documentation, but believe this is should be treated as bug not a feature.

Many thanks to Sujai for openly reporting this via Twitter and posting his BLAST+ bug report as a gist on GitHub where we could see and comment on the issue. Sadly I'm still waiting for an official NCBI public BLAST+ bug tracker.

Update 26 September 2018

This post was cited in Shah et al. (2018) Misunderstood parameter of NCBI BLAST impacts the correctness of bioinformatics workflows. I also fixed a typo.

Update 2 November 2018

I've posted an initial followup, BLAST max alignment limits repartee - part one, which introduces a small self-contained test case, and emphasises that this issue is not specific to -max_target_seqs but also affects the limits -num_descriptions and -num_alignments used with the human readable plain text or HTML output.

In BLAST max alignment limits repartee - part two, I focus on the question of if database order is important (as claimed in Shar et al. 2018), and how exactly the internal alignment number limit works.

Update 13 November 2018

I've published BLAST max alignment limits - part three looking at the test case from Nidhi Shah, and BLAST max alignment limits - part four looking at the internal alignment number limit in the context of nucleotide databases (where composition based statistics are not used).

Update (7 December 2018)

I've published a fifth follow-up post, BLAST tie breaking by database order, looking at how the BLAST database order is defined in comparison to the FASTA file used to build a database.

Update (8 January 2019)

I've published a sixth follow-up post, "An overly aggressive optimization in BLASTN and MegaBLAST", which comes after the BLAST team's formal reply to the Shah et al. (2018) letter and BLAST+ 2.8.1 were published in late December. This update fixed oddities reported in Shah et al. (2018), which turn out to have been a complex interaction of multiple issues.