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.
XP_005244948.1
hit is in the first result) $\endgroup$