alterlab-blast

Files (5)

• SKILL.md
• evals.json · evals/
• blast_cli.md · references/
• diamond.md · references/
• parse_blast_tab.py · scripts/

SKILL.md (first lines)

# 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

```bash
# 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

1. **`-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.
2. **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.
3. **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.
4. **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.
5. **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`](references/blast_cli.md).

## Taxonomic Scoping

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

```bash
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`](references/blast_cli.md#taxonomy).

## 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:

```bash
diamond makedb --in nr.faa -d nr_diamond

From AlterLab-IEU/AlterLab-Academic-Skills, MIT.

More from AlterLab-IEU/AlterLab-Academic-Skills