{
  "skill": "alterlab-blast",
  "evals": [
    {
      "id": "makeblastdb-then-blastp",
      "prompt": "I have a FASTA of ~5,000 bacterial proteins and a handful of query proteins. I want to build a local protein database on my machine and BLAST my queries against it, and get a tab-separated hit table I can load into pandas. How do I do this with the command-line BLAST tools?",
      "expected_output": "Invokes alterlab-blast for the local CLI workflow: builds the database with `makeblastdb -dbtype prot -parse_seqids`, runs `blastp -query ... -db ... -outfmt 6` with an explicit -evalue, and hands the tabular file to scripts/parse_blast_tab.py. Notes that -outfmt 6 yields the std 12 columns and that -parse_seqids is needed for later blastdbcmd retrieval. Does not reach for the Biopython web API.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "makeblastdb" },
        { "type": "behavior", "value": "Builds the DB with makeblastdb -parse_seqids and searches with blastp -outfmt 6, rather than using the web NCBIWWW API." }
      ]
    },
    {
      "id": "max-target-seqs-best-hit-pitfall",
      "prompt": "I ran blastn with -max_target_seqs 1 to get the single best hit per query, but a colleague said my best hits look wrong. Is -max_target_seqs 1 the right way to get the top hit?",
      "expected_output": "Invokes alterlab-blast and corrects the misconception: explains that -max_target_seqs is a heuristic keep-count applied during the search (ties broken by database order), NOT a guaranteed top-N best-hits filter, so -max_target_seqs 1 does not reliably return the best hit. Recommends keeping a generous value and selecting the top row by post-hoc bitscore sort via scripts/parse_blast_tab.py --best-hit.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "max_target_seqs" },
        { "type": "behavior", "value": "States that -max_target_seqs is a heuristic cutoff, not a top-N filter, and recommends sorting by bitscore after the search to get the true best hit." }
      ]
    },
    {
      "id": "diamond-large-protein-search",
      "prompt": "blastp against the full nr database is taking forever for my 2 million predicted metagenome proteins. Is there a faster command-line option that still gives me BLAST-like protein hits in tabular form?",
      "expected_output": "Invokes alterlab-blast and routes to the DIAMOND fast path: `diamond makedb` then `diamond blastp --ultra-sensitive --outfmt 6 ...` for BLAST-comparable recall at large scale. Flags the DIAMOND quoting gotcha (custom --outfmt 6 columns are UNQUOTED, the inverse of BLAST+) and notes DIAMOND is protein-space only. Output stays parseable by scripts/parse_blast_tab.py.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "diamond" },
        { "type": "behavior", "value": "Recommends DIAMOND blastp --ultra-sensitive for the large protein search and warns the custom --outfmt column list must be unquoted for DIAMOND." }
      ]
    },
    {
      "id": "taxid-scoped-blastn",
      "prompt": "I'm BLASTing some contigs against nt but I only care about hits to primates, and I want to exclude bacterial contamination. How do I restrict a command-line blastn search by taxonomy?",
      "expected_output": "Invokes alterlab-blast and uses taxonomic scoping: `-taxids` to restrict to a clade (primates) and `-negative_taxids` to exclude bacteria, run against a taxonomy-aware DB such as NCBI's pre-formatted nt. Explains that subtree expansion is automatic and that a taxonomy-aware database (or build-time taxids) is required.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "taxids" },
        { "type": "behavior", "value": "Uses -taxids/-negative_taxids against a taxonomy-aware DB and notes the database must carry taxonomy mapping." }
      ]
    },
    {
      "id": "near-miss-biopython-web-blast",
      "prompt": "I don't want to install or download any databases. Can you script a quick BLAST of one protein sequence against NCBI's servers from Python using Bio.Blast.NCBIWWW and parse the XML result?",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-biopython. The user explicitly wants the Biopython web BLAST API (Bio.Blast.NCBIWWW.qblast against NCBI servers, no local DB), not local BLAST+ binaries or a local database. alterlab-blast is the command-line/local-database skill and does not run the web NCBIWWW API.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-biopython" }
      ]
    },
    {
      "id": "near-miss-gget-oneliner",
      "prompt": "I just want a fast one-liner to BLAST a single short sequence and also grab the gene's Ensembl info, without setting up any databases. What's the quickest tool for that interactive lookup?",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-gget. The user wants a quick interactive one-liner lookup (gget blast plus Ensembl gene info) with no local database setup, which is gget's niche. alterlab-blast targets local BLAST+ databases and CLI searches, not throwaway one-liner web queries.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-gget" }
      ]
    },
    {
      "id": "near-miss-phylogenetics-tree",
      "prompt": "I have 60 homologous protein sequences and I want to align them and build a maximum-likelihood phylogenetic tree with bootstrap support, then visualize it. How should I do this?",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-phylogenetics. The user wants multiple sequence alignment plus tree inference and visualization, not a similarity/homology search against a database. alterlab-blast only performs BLAST+/DIAMOND sequence search and does not build phylogenetic trees.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-phylogenetics" }
      ]
    }
  ]
}

# BLAST+ 2.17.0 — Command-Line Reference

Deep reference for the AlterLab BLAST+ skill. Pins to NCBI BLAST+ **2.17.0**.
Primary source: NCBI BLAST Command Line Applications User Manual
(https://www.ncbi.nlm.nih.gov/books/NBK569856/) and its options appendix
(https://www.ncbi.nlm.nih.gov/books/NBK279684/).

## Programs

| Program | Query | Subject DB (`-dbtype`) | Typical use |
|---------|-------|------------------------|-------------|
| `blastn` | nucleotide | `nucl` | nucleotide vs nucleotide |
| `blastp` | protein | `prot` | protein vs protein |
| `blastx` | nucleotide (6-frame translated) | `prot` | annotate an unknown ORF/transcript |
| `tblastn` | protein | `nucl` (translated) | find a protein in a genome/transcriptome |

All take `-query`, `-db`, `-out`, `-evalue`, `-outfmt`, `-num_threads`.

## Building a database: `makeblastdb`

```bash
makeblastdb -in seqs.fasta -dbtype {nucl|prot} -parse_seqids \
  -out mydb -title "human description"
```

- **`-parse_seqids`** parses the FASTA deflines into the DB's index so you can
  later retrieve specific records with `blastdbcmd -entry <id>` and so the IDs
  stay clean for DIAMOND reuse. It cannot be added retroactively — rebuild if you
  forgot it.
- `-dbtype nucl` for nucleotide subjects, `prot` for protein.
- BLAST+ **2.17.0** (released 21 July 2025) can read compressed FASTA input
  directly — **gzip (`.gz`), bzip2 (`.bz2`) and zstd (`.zst`)** — choosing the
  decompressor from the file extension of `-in`. Building from source needs the
  matching zlib/bzip2/zstd libraries.
- Multi-FASTA on stdin: `... -in - ...`.

### Retrieving sequences back out

```bash
blastdbcmd -db mydb -entry sp|P01308|INS_HUMAN          # one record
blastdbcmd -db mydb -entry all -outfmt "%f"             # dump all as FASTA
blastdbcmd -db mydb -info                               # DB stats
```

`blastdbcmd -entry` only works if the DB was built with `-parse_seqids`.

## Output formats (`-outfmt`)

| `-outfmt` | Meaning |
|-----------|---------|
| `0` | pairwise (human-readable; default) |
| `6` | **tabular** — the canonical machine-parseable format |
| `7` | tabular with comment lines (`# ` headers + `# N hits found`) |
| `5` | BLAST XML |
| `15` | BLAST JSON (single file) |

Format `6`/`7` default columns (the `std` keyword) are exactly:

```
qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore
```

Customize by appending a **quoted** space-separated keyword list (BLAST+ syntax):

```bash
blastp -query q.faa -db nr -outfmt '6 qseqid sseqid pident length evalue bitscore staxids sscinames' -out hits.tsv
```

Useful extra keywords: `staxids`, `sscinames`, `qcovs` (query coverage per
subject), `qlen`, `slen`, `stitle`. (DIAMOND uses the *unquoted* form — see
`diamond.md`.)

## `-task` (blastn)

`blastn` ships several task presets; the default is `megablast`:

| `-task` | When |
|---------|------|
| `megablast` (default) | highly similar sequences (same species, sequencing-error scale) |
| `dc-megablast` | discontiguous megablast — inter-species comparison |
| `blastn` | traditional, more sensitive to divergent hits |
| `blastn-short` | short queries (< ~30 nt): primers, sgRNA, oligos |

`blastp` likewise has `blastp` (default), `blastp-fast`, `blastp-short`.
Picking the wrong task is a common cause of "BLAST missed an obvious hit".

## `-max_target_seqs` — the heuristic trap

From the options appendix, `-max_target_seqs` is the **"Number of aligned
sequences to keep"** (default **500**), applied with report formats above
`-outfmt 4`. Critically:

- It is a **heuristic limit applied during the search**, not a post-hoc "return
  the N best" filter. The set of sequences kept is **not guaranteed** to be the N
  highest-scoring ones.
- **Ties are broken by order of sequences in the database**, not by score.
- Therefore `-max_target_seqs 1` does **not** reliably return the single best
  hit. To get a true best-hit, leave `-max_target_seqs` generous and select the
  top row *after* sorting by `bitscore` (see `scripts/parse_blast_tab.py`).

## Taxonomy scoping {#taxonomy}

```bash
blastn -query q.fna -db nt -taxids 9606            # restrict to Homo sapiens
blastp -query q.faa -db nr -taxids 2,4751          # bacteria + fungi
blastp -query q.faa -db nr -negative_taxids 2      # exclude bacteria
```

- `-taxids` / `-negative_taxids` take comma-separated NCBI taxids; subtree
  expansion is automatic (a clade taxid covers its descendants).
- `-taxidlist` / `-negative_taxidlist` read taxids from a file.
- Requires a **taxonomy-aware DB**: NCBI's pre-formatted `nt` / `nr` ship with
  the needed taxid mapping; for a custom DB supply taxids at build time
  (`-taxid` / `-taxid_map`) so scoping works.

## Threading: `-num_threads` and `-mt_mode`

- `-num_threads N` — number of CPU threads.
- `-mt_mode` (integer) controls *how* work is split across threads:
  - `0` (**default**) — BLAST selects the method for you from query size, database
    size, program and task. Since the 2.15 release this auto-selection is what NCBI
    recommends, and it is where the "2–10x faster with many queries against a small
    database" speedup comes from.
  - `1` — **ThreadByQuery**: each thread takes a batch of queries and searches the
    whole database. Good for many queries against a relatively small database.
  - `2` — **ThreadByDatabase**: split by database volume. Suits larger databases and
    any number of queries.
- Override the default only when you have measured a reason to; NCBI's manual
  explicitly calls overriding "not recommended".

## Common e-value / filter flags

- `-evalue 1e-5` — expectation-value threshold (lower = stricter).
- `-qcov_hsp_perc 80` — minimum % query coverage per HSP.
- `-perc_identity 90` — (blastn) minimum percent identity.
- `-word_size`, `-gapopen`, `-gapextend`, `-reward`, `-penalty` — alignment
  scoring knobs; defaults differ per `-task`.

## Worked example: annotate unknown transcripts against nr

```bash
# Translated nucleotide query vs protein nr, taxonomy + coverage columns
blastx -query transcripts.fna -db nr \
  -task blastx \
  -evalue 1e-10 -max_target_seqs 50 -num_threads 8 -mt_mode 1 \
  -outfmt '6 qseqid sseqid pident length evalue bitscore staxids sscinames qcovs' \
  -out transcripts_vs_nr.tsv

uv run python scripts/parse_blast_tab.py transcripts_vs_nr.tsv \
  --columns 'qseqid sseqid pident length evalue bitscore staxids sscinames qcovs' \
  --best-hit --min-identity 30 --min-qcov 50
```

# DIAMOND — Fast Protein-Space Search

Reference for the DIAMOND fast path of the AlterLab BLAST+ skill.
Primary source: https://github.com/bbuchfink/diamond and its command-line wiki
(https://github.com/bbuchfink/diamond/wiki/3.-Command-line-options).

## What it is

DIAMOND does pairwise alignment of proteins and translated DNA "at 100x-10,000x
speed of BLAST" (project's own claim). It is a drop-in replacement for `blastp`
and `blastx` (protein subject databases only) when BLAST+ is too slow on large
inputs — think millions of query proteins or nr-scale subject DBs. It does
**not** do nucleotide-vs-nucleotide search; for that stay with `blastn`.

## Build a database

```bash
diamond makedb --in reference.faa -d reference_diamond
# --in / -d ; output gets a .dmnd extension
```

Add `--taxonmap`, `--taxonnodes`, `--taxonnames` (NCBI taxonomy dumps) at build
time to enable taxonomic columns/filters in the output.

## Search

```bash
diamond blastp -d reference_diamond -q query.faa -o hits.tsv \
  --ultra-sensitive \
  --outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore
```

- `--db/-d`, `--query/-q`, `--out/-o` mirror BLAST+.
- `blastx` is also available (`diamond blastx`) for translated-nucleotide queries
  against a protein DB.

### Sensitivity ladder

Fast → most sensitive (mutually exclusive):

`--fast` · `--mid-sensitive` · `--sensitive` · `--more-sensitive` ·
`--very-sensitive` · `--ultra-sensitive`

- Default (no flag) is the fastest mode — fine for high-identity hits, but it
  misses divergent matches.
- `--ultra-sensitive` is the closest to BLAST recall; use it when you would
  otherwise reach for `blastp` and just need it to finish.
- `--sensitive` / `--more-sensitive` are common middle grounds for homology
  searches.

## Output format — the quoting gotcha

DIAMOND's tabular format is `--outfmt 6` (alias `-f 6`) with the same default
columns as BLAST+. **Custom column keywords must be UNQUOTED and
space-separated**, directly after the `6`:

```bash
# CORRECT (DIAMOND):
diamond blastp ... --outfmt 6 qseqid sseqid pident evalue bitscore

# WRONG for DIAMOND (this is the BLAST+ syntax):
diamond blastp ... --outfmt '6 qseqid sseqid pident evalue bitscore'
```

This is the inverse of BLAST+, which wants the spec quoted. Mixing them up is a
frequent silent failure. The resulting `.tsv` is compatible with
`scripts/parse_blast_tab.py`.

## When to choose DIAMOND vs BLAST+

| Situation | Tool |
|-----------|------|
| Nucleotide-vs-nucleotide | BLAST+ `blastn` (DIAMOND can't) |
| Small protein query set, need exact BLAST behavior | BLAST+ `blastp` |
| Millions of proteins / metagenomic ORFs vs nr | DIAMOND `blastp --ultra-sensitive` |
| Translated DNA vs protein at scale | DIAMOND `blastx` |
| Need taxonomy columns | either, but build the DB with the taxonomy maps |

## Version

Pin to a current DIAMOND release — **2.2.8** is the latest on bioconda as of
2026-09. Run `diamond version` and record it alongside results; a DIAMOND
database built by one minor version is not guaranteed readable by another, so
rebuild the `.dmnd` after upgrading.

#!/usr/bin/env python3
"""Parse and QC BLAST+/DIAMOND tabular (`-outfmt 6` / `--outfmt 6`) output.

Self-contained: Python standard library only (no requests/pandas needed), so it
runs in a bare `uv run python` environment. Reads a BLAST+ or DIAMOND tabular
file (or `-` for stdin), applies identity / coverage / e-value filters, and can
extract the true best hit *per query* by sorting on bitscore.

Why this exists: `-max_target_seqs` is a heuristic keep-count, NOT a top-N
best-hits filter (ties are broken by database order, not score). So the only
reliable way to get the best hit is to sort the tabular output by bitscore here,
after the search. This script does that and warns when the column layout looks
inconsistent with the declared schema.

Default schema is the BLAST+/DIAMOND `std` 12 columns:
    qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore
Override with --columns if you searched with a custom -outfmt spec.

Usage:
    parse_blast_tab.py hits.tsv --best-hit
    parse_blast_tab.py hits.tsv --min-identity 90 --min-qcov 80 --max-evalue 1e-10
    parse_blast_tab.py hits.tsv --columns 'qseqid sseqid pident length evalue bitscore qcovs' --best-hit
    blastp ... -outfmt 6 | parse_blast_tab.py - --best-hit --json
"""

from __future__ import annotations

import argparse
import json
import sys
from dataclasses import dataclass, field

STD_COLUMNS = (
    "qseqid sseqid pident length mismatch gapopen "
    "qstart qend sstart send evalue bitscore"
).split()

# Columns we know how to coerce to numbers for filtering / sorting.
NUMERIC = {
    "pident", "length", "mismatch", "gapopen", "qstart", "qend",
    "sstart", "send", "evalue", "bitscore", "qcovs", "qcovhsp",
    "qlen", "slen", "score", "nident", "positive", "gaps",
}


@dataclass
class Result:
    columns: list[str]
    rows: list[dict] = field(default_factory=list)
    warnings: list[str] = field(default_factory=list)
    n_input: int = 0
    n_kept: int = 0


def _coerce(col: str, value: str):
    if col not in NUMERIC:
        return value
    try:
        return float(value)
    except ValueError:
        return None


def parse(handle, columns: list[str]) -> Result:
    res = Result(columns=columns)
    ncol = len(columns)
    for raw in handle:
        line = raw.rstrip("\n")
        if not line or line.startswith("#"):  # outfmt 7 comment lines
            continue
        fields = line.split("\t")
        res.n_input += 1
        if len(fields) != ncol:
            res.warnings.append(
                f"row {res.n_input}: expected {ncol} columns, got {len(fields)} "
                f"-- check your --columns matches the search -outfmt spec"
            )
            # Still record what we can, padding/truncating defensively.
            fields = (fields + [""] * ncol)[:ncol]
        row = {col: _coerce(col, val) for col, val in zip(columns, fields)}
        res.rows.append(row)
    return res


def _pass_filters(row, args) -> bool:
    if args.min_identity is not None:
        v = row.get("pident")
        if not isinstance(v, float) or v < args.min_identity:
            return False
    if args.min_qcov is not None:
        # accept either qcovs or qcovhsp, whichever was requested
        v = row.get("qcovs", row.get("qcovhsp"))
        if not isinstance(v, float) or v < args.min_qcov:
            return False
    if args.max_evalue is not None:
        v = row.get("evalue")
        if not isinstance(v, float) or v > args.max_evalue:
            return False
    return True


def best_hit_per_query(rows: list[dict]) -> list[dict]:
    """One row per qseqid, the highest bitscore (e-value tiebreak). This is the
    correct way to get 'the best hit' -- do NOT rely on -max_target_seqs 1."""
    best: dict[str, dict] = {}
    for row in rows:
        q = row.get("qseqid")
        if q is None:
            continue
        bs = row.get("bitscore")
        ev = row.get("evalue")
        cur = best.get(q)
        if cur is None:
            best[q] = row
            continue
        cur_bs = cur.get("bitscore")
        better = isinstance(bs, float) and (
            not isinstance(cur_bs, float)
            or bs > cur_bs
            or (bs == cur_bs and isinstance(ev, float)
                and isinstance(cur.get("evalue"), float) and ev < cur["evalue"])
        )
        if better:
            best[q] = row
    return list(best.values())


def main(argv=None) -> int:
    p = argparse.ArgumentParser(description=__doc__.split("\n")[0])
    p.add_argument("path", help="tabular BLAST+/DIAMOND file, or '-' for stdin")
    p.add_argument("--columns", default=" ".join(STD_COLUMNS),
                   help="space-separated -outfmt column names (default: std 12)")
    p.add_argument("--best-hit", action="store_true",
                   help="keep only the top-bitscore hit per query (correct best-hit selection)")
    p.add_argument("--min-identity", type=float, metavar="PCT",
                   help="drop hits below this percent identity (pident)")
    p.add_argument("--min-qcov", type=float, metavar="PCT",
                   help="drop hits below this query coverage (needs qcovs/qcovhsp column)")
    p.add_argument("--max-evalue", type=float, metavar="E",
                   help="drop hits with e-value above this threshold")
    p.add_argument("--json", action="store_true", help="emit JSON instead of TSV")
    args = p.parse_args(argv)

    columns = args.columns.split()

    handle = sys.stdin if args.path == "-" else open(args.path, "r", encoding="utf-8")
    try:
        res = parse(handle, columns)
    finally:
        if handle is not sys.stdin:
            handle.close()

    rows = [r for r in res.rows if _pass_filters(r, args)]
    if args.best_hit:
        rows = best_hit_per_query(rows)
        rows.sort(key=lambda r: (-(r.get("bitscore") or 0.0), str(r.get("qseqid"))))
    res.n_kept = len(rows)

    for w in res.warnings:
        print(f"WARNING: {w}", file=sys.stderr)
    print(
        f"# parsed {res.n_input} hit(s); kept {res.n_kept} after filters"
        + (" (best-hit-per-query)" if args.best_hit else ""),
        file=sys.stderr,
    )

    if args.json:
        json.dump(
            {"columns": columns, "n_input": res.n_input, "n_kept": res.n_kept,
             "warnings": res.warnings, "rows": rows},
            sys.stdout, indent=2, default=str,
        )
        sys.stdout.write("\n")
    else:
        sys.stdout.write("\t".join(columns) + "\n")
        for r in rows:
            sys.stdout.write("\t".join(str(r.get(c, "")) for c in columns) + "\n")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())