# Bio Core Blast Searching

> Search protein/nucleotide sequences for homologs with NCBI BLAST+ (blastp/blastn/blastx/tblastn) via Biopython qblast or a local blastdb; parse E-value/bit-score/identity. Use for BLAST search, homology/similarity search, or FASTA annotation.

- Skill: `pavel-kravchenko/bio-core-blast-searching` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-blast-searching`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-blast-searching/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Program | Query | Database | Translation | Use Case |
|---|---|---|---|---|
| **blastn** | Nucleotide | Nucleotide | None | Gene finding, species ID, primer design |
| **blastp** | Protein | Protein | None | Protein homology, function prediction |
| **blastx** | Nucleotide | Protein | Query in 6 frames | Annotate a novel nucleotide sequence (EST/mRNA) |
| **tblastn** | Protein | Nucleotide | DB in 6 frames | Find unannotated protein-coding genes in a genome |
| **tblastx** | Nucleotide | Nucleotide | Both in 6 frames | Compare unannotated genomes (slowest, rarely needed) |

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).

```python
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.

```python
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).

```python
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`