Blast hits disappearing after changing -evalue

Question

I was teaching an introduction to bioinformatics when the students and I noticed strange Blast behavior that we couldn't explain.

With the default evalue parameter, our best hit shows an evalue of 3e-06, while that hit disappears with a evalue threshold of 1e-5 (which is bigger!)

(after a few tries it seems that the evalue of the top hit starts changing with an evalue threshold of 2)

Input proteins:

wget ftp://ftp.ncbi.nih.gov/refseq/B_taurus/mRNA_Prot/cow.1.protein.faa.gz
wget ftp://ftp.ncbi.nih.gov/refseq/H_sapiens/mRNA_Prot/human.1.protein.faa.gz 
gzip -d *.gz
head -199 cow.1.protein.faa > medium_cow.faa

Database building:

makeblastdb -in human.1.protein.faa -dbtype prot

Blast command 1 (default):

blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast

output:

# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 13 hits found
NP_001180091.1  XP_011539315.1  29.41   102     61      3       5       100     144     240     3e-06   42.0

Blast command 2 (-evalue 1e-05):

blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 1e-05

output:

# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# 0 hits found

Blast command 3 (-evalue 2):

blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 2

output:

# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 6 hits found
NP_001180091.1  XP_005244948.1  40.62   32      18      1       4       35      394     424     0.076   29.3

I reproduced this behavior with Blast 2.2.31+, 2.7.1+ and the web blast at https://blast.ncbi.nlm.nih.gov (screenshots attached)

Anything obvious that my students and I are missing? To my knowledge changing the evalue threshold should no change the values themselves, it should just change the reporting.

Answer 1

I am the author of the question / answer here Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?

To my knowledge changing the evalue threshold should no change the values themselves, it should just change the reporting.

That is sadly not the case, though thats what I would say is typically expected.

Yours was the exact situation that lead me to question the implementation of NCBI blast. I got results with an evalue above a certain value with a relaxed threshold, but with a stricter threshold (still lower than the value of the results I got) those results disappeared.

From what I found all parameters imposing some threshold show this behavior and are applied in an early step of the program leading to prematurely removed results that later on might score higher than those thresholds. (Say your hit first has an evalue of 5: threshold 10 = ok and kept, threshold 1e-5 = gets removed, but that same hit gets improved further along in the program and then gets an evalue of 1e-6. Now this is only shown when you use an evalue threshold of 5 or higher, though 1e-6 is clearly also a result you would want to see when using a threshold of 3, 1, 0.01, ... etc.)

I was relatively shocked by this and the fact that its basically nowhere officially documented and generally a not known fact what influence these options have.

Sadly for now the stance of the developers is basically all is working as intended and the only feedback I got amounts to

The safest approach is to relax (increase) the threshold and deal with a larger result set.

aka Suck it up and deal with it.