Bio.Blast

Core Principle: Bio.Blast interfaces with NCBI QBLAST servers for remote searches and provides structured XML parsers for local alignment outputs, where local binary execution is orchestrated via subprocess.

What is gone

Removed in Biopython 1.86 along with Bio.Application:

from Bio.Blast.Applications import NcbiblastnCommandline     # ImportError
NcbiblastnCommandline(query="q.fa", db="nt", out="r.xml")     # gone

Every pre-2026 BLAST tutorial uses this. It is dead code. See Biopython Deprecations.

Running local BLAST now

import subprocess
from pathlib import Path

def blastp(query, db, out, evalue=1e-5, threads=4, outfmt=5):
    subprocess.run([
        "blastp",
        "-query", str(query),
        "-db", str(db),
        "-out", str(out),
        "-evalue", str(evalue),
        "-num_threads", str(threads),
        "-outfmt", str(outfmt),
    ], check=True)
    return Path(out)

check=True raises on failure. Without it, a BLAST that errored leaves you parsing an empty file and finding no hits: which looks like a biological result.

Useful -outfmt values:

outfmt	Format
5	XML: richest, parseable by Biopython
6	tabular: parse with pandas
7	tabular with comment lines
0	pairwise, human-readable

For most work, use -outfmt 6 and pandas. It is faster to parse and easier to filter than XML.

BLAST6 = ["qseqid", "sseqid", "pident", "length", "mismatch", "gapopen",
          "qstart", "qend", "sstart", "send", "evalue", "bitscore"]

hits = pd.read_csv("results.tsv", sep="\t", names=BLAST6)
best = hits.sort_values("bitscore", ascending=False).drop_duplicates("qseqid")

You can request extra fields:

-outfmt "6 qseqid sseqid pident length evalue bitscore qcovs staxids sscinames"

qcovs (query coverage) is essential: a 100%-identity hit over 20 amino acids is not homology, and percent identity alone will not tell you that.

Parsing XML

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 132, in run_cell_capture
    exec(compile(code, "<cell>", "exec"), namespace)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<cell>", line 3, in <module>
FileNotFoundError: [Errno 2] No such file or directory: 'results.xml'

The three-level hierarchy: record → alignment → HSP: is the thing to internalise. One query produces many subject alignments; each alignment can contain several high-scoring segment pairs (separate local alignments, e.g. one per exon or per repeated domain).

Newer Biopython also has Bio.Blast.parse / Bio.Blast.read with an updated object model, and Bio.SearchIO provides a unified interface across BLAST, HMMER, BLAT and others:

Traceback (most recent call last):
  File "/Users/leviosman/Developer/python-for-bio/webapp/scripts/curriculum/generate_cell_outputs.py", line 132, in run_cell_capture
    exec(compile(code, "<cell>", "exec"), namespace)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<cell>", line 2, in <module>
  File "/Library/Frameworks/Python.framework/Versions/3.14/lib/python3.14/site-packages/Bio/SearchIO/__init__.py", line 309, in parse
    with as_handle(handle, **handle_kwargs) as source_file:
         ~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Library/Frameworks/Python.framework/Versions/3.14/lib/python3.14/contextlib.py", line 141, in __enter__
    return next(self.gen)
  File "/Library/Frameworks/Python.framework/Versions/3.14/lib/python3.14/site-packages/Bio/File.py", line 72, in as_handle
    with open(handleish, mode, **kwargs) as fp:
         ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^
FileNotFoundError: [Errno 2] No such file or directory: 'results.xml'

SearchIO is the more consistent API if you work with several search tools.

Remote BLAST

from Bio.Blast import NCBIWWW

with NCBIWWW.qblast("blastn", "nt", sequence_string, expect=1e-5,
                    hitlist_size=50) as h:
    xml = h.read()
Path("results.xml").write_text(xml)

Slow: minutes per query: and rate-limited. Acceptable for one or two exploratory queries. Never put it in a loop. For anything systematic, download the database and run locally, or use a faster tool.

Building a local database

makeblastdb -in proteins.faa -dbtype prot -out mydb -parse_seqids
blastp -query query.faa -db mydb -outfmt 6 -evalue 1e-5 -num_threads 8 -out hits.tsv

-parse_seqids lets you later retrieve sequences from the database with blastdbcmd.

Filtering hits sensibly

hits = pd.read_csv("hits.tsv", sep="\t", names=BLAST6 + ["qcovs"])

good = hits.query("evalue < 1e-10 and pident > 30 and qcovs > 50")

# best hit per query
best = hits.sort_values(["qseqid", "bitscore"], ascending=[True, False]) \
           .drop_duplicates("qseqid")

# reciprocal best hits: the standard orthology heuristic
def rbh(ab, ba):
    a = ab.sort_values("bitscore", ascending=False).drop_duplicates("qseqid")
    b = ba.sort_values("bitscore", ascending=False).drop_duplicates("qseqid")
    m = a.merge(b, left_on=["qseqid", "sseqid"], right_on=["sseqid", "qseqid"],
                suffixes=("_ab", "_ba"))
    return m[["qseqid_ab", "sseqid_ab", "bitscore_ab", "bitscore_ba"]]

Filter on E-value AND coverage, not identity alone. The 30%/50% thresholds above are conventional starting points, not laws.

E-value depends on database size. The same alignment gets a worse E-value in a larger database. Bit scores are database-independent and therefore comparable across searches: prefer them when ranking across different runs.

Faster alternatives

BLAST is slow by modern standards. For large searches:

Tool	Use
DIAMOND	protein search, 100–10,000× faster than BLASTP, BLAST-compatible output
MMseqs2	very fast search and clustering
minimap2	long reads and whole-genome nucleotide alignment
HMMER	profile searches: far more sensitive for remote homology
foldseek	structure-based search; finds homologues sequence methods miss

DIAMOND is a drop-in replacement for most BLASTP use:

diamond makedb --in proteins.faa -d mydb
diamond blastp -q query.faa -d mydb -o hits.tsv --outfmt 6 --evalue 1e-5 --threads 8

Same output format, so the pandas code above works unchanged.

Bioinformatics examples

# annotate assembled contigs against a reference proteome
subprocess.run(["diamond", "blastx", "-q", "contigs.fasta", "-d", "swissprot",
                "-o", "annot.tsv", "--outfmt", "6", "qseqid", "sseqid",
                "pident", "length", "evalue", "bitscore", "qcovhsp",
                "--evalue", "1e-10", "--threads", "8"], check=True)

annot = pd.read_csv("annot.tsv", sep="\t",
                    names=["qseqid","sseqid","pident","length","evalue","bitscore","qcov"])
best = annot.sort_values("bitscore", ascending=False).drop_duplicates("qseqid")
print(f"{best['qseqid'].nunique():,} of {n_contigs:,} contigs annotated")

# check for contamination: what organism do unannotated reads hit?
hits = pd.read_csv("nt_hits.tsv", sep="\t",
                   names=BLAST6 + ["staxids", "sscinames"])
hits["sscinames"].value_counts().head(20)

# verify a primer is specific
subprocess.run(["blastn", "-query", "primers.fa", "-db", "hg38",
                "-task", "blastn-short", "-outfmt", "6",
                "-evalue", "1000", "-out", "primer_hits.tsv"], check=True)

-task blastn-short is required for sequences under ~30 bp: the default blastn word size will find nothing for a 20-mer primer.

Common mistakes

• Bio.Blast.Applications. Removed in 1.86.
• NCBIWWW.qblast in a loop. Minutes per query and rate-limited.
• Filtering on percent identity alone, without coverage.
• Comparing E-values across differently-sized databases. Use bit scores.
• Taking the top hit as the answer: check whether hits 2 and 3 score nearly as well.
• subprocess.run without check=True, so a failed BLAST looks like zero hits.
• Parsing XML when tabular would do.
• Forgetting -task blastn-short for primers and short oligos.
• Not using DIAMOND for large protein searches.

See also

PairwiseAligner · Substitution Matrices · Bio.Entrez · Biopython Deprecations · SeqIO · DataFrame

Exercise

The lesson recommends -outfmt 6 and a table reader. Biopython is not installed here, and neither is pandas needed: parse the tabular output with str.split so you can see exactly which column is which.

Write best_hits(cutoff, lines). Each data line holds the twelve standard BLAST tabular fields separated by tabs: qseqid, sseqid, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore. Skip blank lines, lines beginning with #, and any line with fewer than twelve fields. Discard every hit whose E-value is greater than cutoff, then keep for each query the surviving hit with the highest bit score, breaking ties in favour of the subject id that sorts first. Print one line per query in ascending query-id order as qseqid, sseqid and the bit score rounded to one decimal place, tab separated. If no hit survives, print no hits instead. The cutoff arrives as the first line of input and the hit table follows.