Bio Blast Searches

Version Compatibility

Reference examples tested with: BioPython 1.83+, NCBI BLAST+ 2.15+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show biopython then help(Bio.Blast.NCBIWWW.qblast) to check signatures
• CLI: blastn -version then blastn -help

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

BLAST Searches (Remote)

"Find similar sequences in NCBI's database" -> Submit a query to NCBI's remote BLAST servers; receive a Request ID (RID); poll for completion; parse the XML hit table. Best for one-off identification of a few sequences. For >50 sequences, switch to local-blast or DIAMOND/MMseqs2 in remote-homology.

The two most consequential decisions: which program (defines query+target molecule types and word-size defaults) and which database (defines the search space and therefore E-value baselines). The third most important: do NOT misuse max_target_seqs -- it is an early-termination heuristic, not a "give me the top N hits" filter (Shah et al. 2019).

• Python: NCBIWWW.qblast(program, db, sequence) + NCBIXML.read(handle) (BioPython)
• CLI: blastn -remote -db nt -query seq.fa -out hits.xml -outfmt 5 (BLAST+)
• Web: https://blast.ncbi.nlm.nih.gov/Blast.cgi (RID lookup)

Required Setup

from Bio.Blast import NCBIWWW, NCBIXML
from Bio import SeqIO

No API key needed for remote BLAST itself, but NCBI's general rate-limit ethic still applies -- one search at a time, polite waiting, no parallelism.

Program decision (query vs database molecule)

The misuse to avoid: using default blastn (word=11) for cross-species DNA where dc-megablast is the right tool. Or using megablast (word=28) for cross-species homology where it will miss every divergent hit. The most-misused BLAST parameter according to literature.

Database decision (search space)

For publication reproducibility, never search nt or nr without recording the snapshot date and ideally archiving a frozen copy. Default to refseq_select for any cross-species homology question; switch to nt/nr only when curated coverage is insufficient.

E-value interpretation (Karlin-Altschul)

E-value = K * m * n * exp(-lambda * S), where m = effective query length, n = effective database size, lambda and K are scoring-matrix-dependent constants (Karlin & Altschul 1990 PNAS 87:2264).

Key implication of E = K * m * n * exp(-lambda * S): the same alignment against a 100x larger database has a 100x larger E-value. Cross-database E-value comparison is meaningless. Bit-score is database-size normalized and is the right cross-database metric.

For protein remote homology where E is marginal (10^-3 to 10^-1), reach for profile methods: PSI-BLAST, jackhmmer, HHblits, or Foldseek -- see remote-homology skill.

Composition-Based Statistics (CBS)

Compositional bias inflates significance for low-complexity proteins. The CBS modes (Yu et al. 2006 Nucleic Acids Res 34:5966):

For protein queries under 30 aa, switch to CBS=3. For protein with known compositional bias (e.g. coiled-coil regions, signal peptides), CBS=2 is appropriate but consider hard-masking with SEG.

The max_target_seqs trap

The misuse: max_target_seqs=10 is interpreted as "return the 10 most significant hits". It is not. The flag is an early termination parameter that affects which hits the search ever considers, not which it ultimately reports (Shah N, Nute MG, Warnow T, Pop M. (2019) Misunderstood parameter of NCBI BLAST impacts the correctness of bioinformatics workflows. Bioinformatics 35:1613-1614).

Consequences:

• Setting max_target_seqs=10 can return entirely different hits than max_target_seqs=500 then filtering to top 10 by E-value.
• The "top 10" by E-value as reported may not be the actual top 10.

Correct pattern: set hitlist_size (Bio.Blast parameter name) large (1000+), then post-filter to the top N by E-value or bit-score in Python.

Word size, gap costs, and matrix

For very short query proteins (e.g. proteomics-identified peptides), BLOSUM45 + word=2 + PAM30 substitution matrix is more sensitive than the default. Use matrix='PAM30' for searches against swissprot.

RID lifecycle

NCBIWWW.qblast() handles polling internally with a fixed retry interval. For long-running searches (>5 min) or batches, submit and capture the RID, then poll independently to avoid blocking. The RID is visible at https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=... for 24-36 hours.

Code patterns

Standard remote BLASTN with reproducible parameters

Goal: Run BLASTN with explicit, paper-quality parameters.

Approach: Specify program, database (refseq_select for stability), word size, expect, and a large hitlist_size to dodge the max_target_seqs trap.

Reference (BioPython 1.83+):

from Bio.Blast import NCBIWWW, NCBIXML

handle = NCBIWWW.qblast(
    program='blastn',
    database='refseq_select_rna',
    sequence=query_seq,
    expect=1e-10,
    word_size=11,
    hitlist_size=500,  # large; filter top-N downstream
    format_type='XML',
)
record = NCBIXML.read(handle); handle.close()
top10 = sorted(record.alignments, key=lambda a: a.hsps[0].expect)[:10]