6
$\begingroup$

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)

hits

no hits

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.

$\endgroup$
6
  • 1
    $\begingroup$ Possible duplicate of Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?. See the end of the answer to this question. e-value doesn't filter only, it influence the results. $\endgroup$
    – llrs
    Jan 26 '18 at 15:07
  • $\begingroup$ @Llopis that answer discusses parameters affecting the number of sequences to show. I don't see why the e-value would change here. Note that the actual e-value is different depending on the threshold set. The size of the DB, however, remains the same so it is very surprising that both the e-value and the bit-scores are different when using different e-value thresholds. These are two completely different hits, not the same hit shown or not shown depending on the options. Are you sure it's the same issue? What am I missing? $\endgroup$
    – terdon
    Jan 26 '18 at 16:34
  • 1
    $\begingroup$ That answers says that evalue parameter is not a threshold is an internal filter!! so using them alters the results and thus the sequences and the e-values: " -evalue -max_target_seqs -num_alignments -num_descriptions are NOT just filters applied to the results, but instead directly influence the results found. " $\endgroup$
    – llrs
    Jan 26 '18 at 16:39
  • $\begingroup$ Because the space searched is different it changes everything of the hits as far as I understood, not only if a sequence is a hit but also the e-value of the hits. But probably the best way to be sure it is a duplicate or not would be to contact them. (Don't cast your binding vote of duplicate if you aren't sure, perhaps it is not an exact duplicate but it is a highly related question) $\endgroup$
    – llrs
    Jan 26 '18 at 16:44
  • 1
    $\begingroup$ If this is just a case of sequences being filtered out, I think it is the same issue. The filter is applied to evalues calculated during the ungapped extension phase, so it's not surprising that they are different. If the actual reported evalue is different, that's another issue, but I don't see any evidence from the question that it is though (it would be nice to know what the evalue of the XP_005244948.1 hit is in the first result) $\endgroup$ Jan 26 '18 at 16:58
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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