alterlab-blast

BLAST+ — Command-Line Sequence Search

Run local NCBI BLAST+ 2.17.0 searches end-to-end: build a database with makeblastdb, search it with blastn / blastp / blastx / tblastn, emit machine-parseable tabular output, and scope by taxonomy. For very large protein searches, hand off to DIAMOND blastp --ultra-sensitive (100x–10,000x the speed of BLAST, per the DIAMOND project). This is the CLI / local-database skill; it is deliberately distinct from the Biopython web API and the gget one-liner (see routing table below).

Bulk DB builds and large searches are CPU/IO-heavy and fully offline — good candidates to run on local compute rather than burning API calls.

When to Use This Skill

Use this skill when the request involves any of:

"BLAST these sequences", "run blastn/blastp/blastx/tblastn", "command-line BLAST"
"build a local BLAST database", "makeblastdb", "index this FASTA for BLAST"
"search my reads against a local nt/nr database", "get tabular BLAST hits I can parse"
"scope the BLAST search to a taxon" (-taxids / -negative_taxids)
"BLAST is too slow on millions of proteins" → DIAMOND blastp
retrieving sequences out of a BLAST DB (blastdbcmd, requires -parse_seqids)

Does NOT Trigger

Route adjacent requests to the right sibling skill instead of forcing BLAST+:

The request is really about…	Route to
The web BLAST API (`Bio.Blast.NCBIWWW.qblast`), or scripting BLAST inside a Python pipeline with `Bio.Blast` parsing	`alterlab-biopython`
A quick one-liner BLAST/database lookup (`gget blast`, gene/structure/enrichment lookups)	`alterlab-gget`
Unified programmatic access to many bio web services (UniProt, KEGG, Ensembl REST, NCBI eUtils)	`alterlab-bioservices`
Building/searching a phylogenetic tree from sequences, not a similarity search	`alterlab-phylogenetics`
Read alignment to a reference genome (BWA/minimap2 → BAM) and SAM/BAM handling	`alterlab-pysam`
FASTQ→VCF variant calling pipeline	`alterlab-nf-core-sarek`
Transcript-level RNA-seq quantification (salmon/kallisto)	`alterlab-rnaseq-quant`
16S/ITS amplicon classification (QIIME 2)	`alterlab-qiime2-amplicon`
Protein structure prediction / embeddings (ESM, AlphaFold)	`alterlab-esm`

If the user explicitly says "web BLAST", "NCBIWWW", or "without installing anything", they want alterlab-biopython, not this skill.

Quick Start

# 1. Build a protein DB (‑parse_seqids enables blastdbcmd retrieval + DIAMOND reuse)
makeblastdb -in proteins.fasta -dbtype prot -parse_seqids -out mydb -title "my proteins"

# 2. Search, tabular output you can parse, std 12 columns
blastp -query query.faa -db mydb -outfmt 6 -evalue 1e-5 -out hits.tsv

# 3. QC / summarize the tabular output (stdlib only)
uv run python scripts/parse_blast_tab.py hits.tsv --best-hit

-outfmt 6 is the canonical machine-readable format; its default columns are the std set: qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore. Use -outfmt 7 for the same columns plus comment lines.

Choosing the Right Program

Query	Subject DB	Program
nucleotide	nucleotide	`blastn`
protein	protein	`blastp`
nucleotide (translated)	protein	`blastx`
protein	nucleotide (translated)	`tblastn`

-dbtype for makeblastdb is nucl for nucleotide subjects, prot for protein.

The Five Things People Get Wrong

-max_target_seqs is NOT a "top N best hits" filter. It is the number of aligned sequences to keep, applied during the search as a heuristic cutoff; ties are broken "by order of sequences in the database", not by score. Setting -max_target_seqs 1 does not reliably return the single best hit. To get the best hit, keep a generous value and pick the top row after sorting by bitscore (see scripts/parse_blast_tab.py --best-hit). Default is 500.
Wrong -task for blastn. megablast (default) is for highly similar sequences; use blastn for cross-species / more divergent hits and blastn-short for queries < ~30 nt (primers, sgRNAs). dc-megablast is the discontiguous option for inter-species comparison.
Forgetting -parse_seqids at DB-build time. Without it you cannot pull sequences back out with blastdbcmd -entry, and DIAMOND cannot reuse the sequence IDs cleanly. You cannot add it later without rebuilding.
Quoting the -outfmt custom column list for DIAMOND. BLAST+ wants the spec quoted (-outfmt '6 qseqid sseqid pident evalue'); DIAMOND wants it unquoted (--outfmt 6 qseqid sseqid pident evalue). Mixing these up is a common silent error.
Multithreading. Use -num_threads N. Since BLAST+ 2.15 the default -mt_mode 0 means BLAST picks the split for you from query and database size, which NCBI recommends leaving alone. Override only deliberately: -mt_mode 1 = ThreadByQuery (many small queries), -mt_mode 2 = ThreadByDatabase (few large queries, big DB).

Full option reference, taxonomy scoping, and DB-prep details: references/blast_cli.md.

Taxonomic Scoping

Restrict a search to (or away from) clades by NCBI taxid:

blastn -query q.fna -db nt -taxids 9606 -outfmt 6 -out human_only.tsv
blastp -query q.faa -db nr -negative_taxids 2 -outfmt 6 -out no_bacteria.tsv

Scoping by taxid requires a taxonomy-aware database (one built/downloaded with its *.taxid mapping, e.g. NCBI's pre-formatted nt / nr). See references/blast_cli.md.

DIAMOND — Fast Path for Large Protein Searches

When blastp / blastx against millions of proteins is too slow, DIAMOND is a drop-in for protein-space search:

diamond makedb --in nr.faa -d nr_diamond
diamond blastp -d nr_diamond -q query.faa -o hits.tsv \
  --ultra-sensitive --outfmt 6 qseqid sseqid pident length evalue bitscore

Sensitivity ladder (fast → most sensitive): --fast, --mid-sensitive, --sensitive, --more-sensitive, --very-sensitive, --ultra-sensitive. With no sensitivity flag DIAMOND runs its default mode, which sits between --fast and --mid-sensitive. Use --ultra-sensitive when you need BLAST-comparable recall. DIAMOND's --outfmt 6 is compatible with the BLAST+ tabular parser below. Details and tradeoffs: references/diamond.md.

Recommended Workflow

Pick the program from the query/subject table above.
Build the DB with makeblastdb -parse_seqids (or download a pre-formatted NCBI DB). For >~1M proteins, build a DIAMOND DB instead.
Search with -outfmt 6, an explicit -evalue threshold, the right -task (blastn), and -num_threads. Add -taxids if scoping.
Parse & QC with scripts/parse_blast_tab.py — it sorts by bitscore, extracts best-hit-per-query, applies identity/coverage/e-value filters, and flags the -max_target_seqs pitfall if the column count looks truncated.
Retrieve any hit sequence with blastdbcmd -db mydb -entry <id> (needs -parse_seqids).

Verify Before Reporting

Confirm blastn -version / diamond version actually ran — never report hits you did not produce.
State the program, -task, -evalue, and DB used; results are meaningless without them.
If you used -max_target_seqs, confirm best-hit selection was done by post-hoc bitscore sort, not by trusting the keep-count as a top-N.
For DIAMOND results, note the sensitivity level used.

References

references/blast_cli.md — full BLAST+ 2.17.0 option reference: programs, makeblastdb, -outfmt columns, -task, taxonomy scoping, -mt_mode, blastdbcmd retrieval, and the -max_target_seqs caveat.
references/diamond.md — DIAMOND DB build, sensitivity modes, output formats, and when to choose it over BLAST+.
NCBI BLAST+ manual: https://www.ncbi.nlm.nih.gov/books/NBK569856/
DIAMOND: https://github.com/bbuchfink/diamond

Part of the AlterLab Academic Skills suite.