Bio Core Blast Searching

BLAST: Sequence Similarity Searching

When to Use

• Identifying an unknown protein or DNA sequence by homology to known sequences
• Annotating genes/ESTs in a new genome or transcriptome assembly
• Finding orthologs/paralogs of a gene across species
• Deciding which BLAST program (blastn/blastp/blastx/tblastn/tblastx) fits a query-database combination
• Running searches locally (no internet, custom/private database) instead of against NCBI's servers

Version Compatibility

• NCBI BLAST+ ≥2.14 (blastp, blastn, blastx, tblastn, makeblastdb CLI tools)
• Biopython ≥1.81 for Bio.Blast.NCBIWWW/NCBIXML. The NcbiblastpCommandline wrapper is deprecated since Biopython 1.82 — build the command and run it via subprocess instead (shown below)
• Python ≥3.9, pandas ≥1.5, numpy, matplotlib

Prerequisites

• pip install biopython pandas numpy matplotlib
• Local searches need the BLAST+ binaries: conda install -c bioconda blast (remote-only workflows just need biopython + internet)
• Familiarity with FASTA format and substitution matrices (BLOSUM62)

Choosing a Program

Within blastn: megablast (word size 28, >95% identity, same species) → blastn (word size 11, cross-species) → discontiguous megablast (~80% identity) as sequences diverge further.

Goal: pick the right program/database combo programmatically instead of guessing. Approach: a small lookup table keyed on (query_type, db_type, sensitivity).

def recommend_blast(query_type, db_type, sensitivity="standard"):
    """Recommend the appropriate BLAST program for a query/database combo.

    query_type, db_type: "nucleotide" or "protein"
    sensitivity: "standard" | "high_identity" | "sensitive" | "very_sensitive"
    """
    recommendations = {
        ("nucleotide", "nucleotide"): {
            "high_identity": ("megablast", "Best for >95% identity, same species"),
            "standard": ("blastn", "Standard nucleotide search"),
            "sensitive": ("discontiguous megablast", "Cross-species, ~80% identity"),
            "very_sensitive": ("tblastx", "Both translated -- slowest but most sensitive"),
        },
        ("nucleotide", "protein"): {"standard": ("blastx", "Translates query in 6 frames")},
        ("protein", "protein"): {
            "standard": ("blastp", "Standard protein search"),
            "sensitive": ("PSI-BLAST", "Iterative, for remote homologs"),
        },
        ("protein", "nucleotide"): {"standard": ("tblastn", "Translates DB in 6 frames")},
    }
    options = recommendations.get((query_type, db_type))
    if options is None:
        return "Invalid combination"
    prog, desc = options.get(sensitivity, options["standard"])
    return f"{prog}: {desc}"


print(recommend_blast("protein", "protein", "sensitive"))
# -> "PSI-BLAST: Iterative, for remote homologs"

Running a Remote Search and Parsing Hits

Goal: search a query sequence against an NCBI database and get a tidy hit table. Approach: NCBIWWW.qblast submits the search and blocks until results are ready (can take 30-120s); NCBIXML.read parses the returned XML into a Bio.Blast.Record with one alignment per hit and one or more hsps (high-scoring pairs) per alignment.

from Bio.Blast import NCBIWWW, NCBIXML
import pandas as pd


def run_blast(sequence, program="blastp", database="swissprot",
              evalue=0.001, hitlist_size=20):
    """Run a BLAST search against NCBI and return a parsed Bio.Blast.Record.

    program: blastn, blastp, blastx, tblastn, tblastx
    database: nr, swissprot, nt, pdb, refseq_protein, ...
    """
    result_handle = NCBIWWW.qblast(
        program=program, database=database, sequence=sequence,
        expect=evalue, hitlist_size=hitlist_size,
    )
    blast_record = NCBIXML.read(result_handle)
    result_handle.close()
    return blast_record


def blast_results_to_dataframe(blast_record):
    """Flatten a BLAST record's alignments/HSPs into a pandas DataFrame."""
    rows = []
    for alignment in blast_record.alignments:
        for hsp in alignment.hsps:
            identity_pct = 100 * hsp.identities / hsp.align_length
            coverage = 100 * (hsp.query_end - hsp.query_start + 1) / blast_record.query_length
            rows.append({
                "hit_id": alignment.hit_id,
                "description": alignment.hit_def[:70],
                "e_value": hsp.expect,
                "bit_score": round(hsp.bits, 1),
                "identity_pct": round(identity_pct, 1),
                "coverage_pct": round(coverage, 1),
                "gaps": hsp.gaps,
                "align_length": hsp.align_length,
            })
    return pd.DataFrame(rows)


# record = run_blast(query_sequence, "blastp", "swissprot", hitlist_size=15)
# df = blast_results_to_dataframe(record)
# df.sort_values("e_value").head(10)

Interpretation guide once you have df: E < 1e-10 = confident homolog; identity > 30% = above the "twilight zone" for proteins; coverage > 70% = full-length match rather than a shared domain only.

Running BLAST Locally (custom/private database)

Goal: search against a private FASTA collection without querying NCBI, and get tab-separated output ready for pandas. Approach: makeblastdb builds the index once; blastp -outfmt 6 emits one row per HSP in a fixed column order — call both via subprocess (the old NcbiblastpCommandline wrapper is deprecated).

import subprocess
import pandas as pd
from io import StringIO

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


def build_blast_db(fasta_path, db_type="prot"):
    """Index a FASTA file for local BLAST searches (dbtype: 'prot' or 'nucl')."""
    subprocess.run(
        ["makeblastdb", "-in", fasta_path, "-dbtype", db_type, "-out", fasta_path],
        check=True, capture_output=True, text=True,
    )


def run_local_blastp(query_fasta, db_path, evalue=1e-3, threads=4):
    """Run blastp against a local database and return hits as a DataFrame."""
    result = subprocess.run(
        ["blastp", "-query", query_fasta, "-db", db_path,
         "-evalue", str(evalue), "-outfmt", "6", "-num_threads", str(threads)],
        check=True, capture_output=True, text=True,
    )
    return pd.read_csv(StringIO(result.stdout), sep="\t", names=BLAST_TABULAR_COLUMNS)


# build_blast_db("proteins.fasta", db_type="prot")
# hits = run_local_blastp("query.fasta", "proteins.fasta")
# hits.sort_values("evalue").head(10)

How BLAST Achieves Its Speed

BLAST skips full Smith-Waterman on every query-database pair by seeding then extending: (1) break the query into words (W=3 for blastp, W=11 for blastn) and, for protein BLAST, find all database words scoring ≥ threshold T against BLOSUM62 (the "word neighborhood"); (2) extend each seed ungapped in both directions, stopping once the score drops more than X (the X-drop) below the best score seen; (3) run full gapped dynamic-programming alignment only on the surviving high-scoring seeds. This is why BLAST can miss true homologs: no seed forms in a highly diverged region, or the X-drop cuts extension short.

Pitfalls

• E-value depends on database size: E=1e-5 in SwissProt (2e8 residues) is far more significant than E=1e-5 in nr (9e10 residues) — always report which database you searched.
• E-value is not a p-value: it is the expected number of chance hits of that score, not a false-positive probability. Use E < 1e-5 as a rough default, but weigh alignment length and percent identity too.
• Low-complexity filtering (SEG/DUST): poly-A stretches, coiled-coils, and other low-complexity regions generate spurious hits — leave the default filter on unless you have a specific reason not to.
• PSI-BLAST convergence/contamination: a false positive entering the PSSM in an early iteration biases every later iteration — verify hits biologically before trusting later rounds.
• NCBI qblast rate limits: remote searches are throttled and can take 30-120s each; batch many queries by building a local database instead of looping NCBIWWW.qblast.
• Coordinate off-by-one: BLAST query/subject coordinates in XML and outfmt 6 are 1-based inclusive, unlike BED's 0-based half-open — mixing these causes off-by-one errors downstream.

See Also

• bio-database-access-local-blast
• bio-database-access-blast-searches
• bio-database-access-sequence-similarity
• bio-alignment-pairwise-alignment