Bio.Blast

In one line: parse BLAST output and (slowly) run remote searches — but Bio.Blast.Applications was removed in 1.86, so you run local BLAST with 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

from Bio.Blast import NCBIXML

with open("results.xml") as h:
    for record in NCBIXML.parse(h):          # one per query
        for alignment in record.alignments:   # one per subject
            for hsp in alignment.hsps:        # one per local alignment
                if hsp.expect < 1e-10:
                    print(alignment.title, hsp.expect, hsp.identities, hsp.align_length)
                    print(hsp.query[:60])
                    print(hsp.match[:60])
                    print(hsp.sbjct[:60])

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:

from Bio import SearchIO
for qresult in SearchIO.parse("results.xml", "blast-xml"):
    for hit in qresult:
        for hsp in hit:
            print(qresult.id, hit.id, hsp.evalue, hsp.ident_num)

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