bio-remote-homology

What this skill does

## Version Compatibility

Reference examples tested with: NCBI BLAST+ 2.15+, HMMER 3.4+, MMseqs2 15+, DIAMOND 2.1+, HH-suite3 3.3+, Foldseek 9+

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then introspect signatures

If a flag is unrecognized or behavior changes, introspect with `--help` and adapt the example to match the installed version rather than retrying.

# Remote Homology

**"Find homologs my BLAST missed"** -> Standard BLAST detects similarity reliably down to ~35% pairwise identity (the "twilight zone", Rost 1999 *Protein Eng* 12:85). Below that, profile methods (PSSMs, HMMs) and structure-aware methods (Foldseek) recover homologs that pairwise alignment misses.

This skill covers the decision: which method, when, against what database. The competition has shifted substantially since 2015: PSI-BLAST is no longer the de-facto standard; MMseqs2 and DIAMOND have replaced BLAST in most large-scale workflows; Foldseek (van Kempen et al. 2024 *Nat Biotechnol* 42:243) detects homologs no sequence method can reach by searching with a 3Di structural alphabet derived from AlphaFold/ESMFold predictions.

- CLI: `psiblast`, `jackhmmer`, `hmmsearch`, `hhblits`, `mmseqs`, `diamond`, `foldseek`
- Python: `Bio.SearchIO` for output parsing; tool-specific clients exist but subprocess is preferred
- Web: HHpred (HHblits webserver), Foldseek webserver, ColabFold for paired structure search

## Required Setup

```bash
# Install via conda
conda install -c bioconda hmmer mmseqs2 diamond hhsuite foldseek
# BLAST+ (separate)
conda install -c bioconda blast

# Verify
hmmsearch -h | head -3       # HMMER 3.4+
mmseqs version               # MMseqs2 15+
diamond --version            # DIAMOND 2.1+
hhblits -h | head -3         # HH-suite3 3.3+
foldseek --version           # Foldseek 9+
```

## Decision matrix: which method when

| Question | Best tool | Why | Sensitivity / Speed |
|---|---|---|---|
| Quick all-vs-all proteome | MMseqs2 or DIAMOND | 100-10,000x faster than BLAST at comparable sensitivity | Highest throughput, near-BLAST sensitivity |
| Identify distant protein homolog (single query) | jackhmmer | Iterative HMM; usually beats PSI-BLAST | Higher sensitivity than PSI-BLAST |
| Distant homology where structure available | Foldseek | 3Di alphabet finds homologs sequence misses | Finds hits PSI-BLAST/HMMER cannot |
| Profile-profile comparison (PDB70 / Pfam) | HHblits + HHsearch | Profile vs profile is most sensitive when target also has profile | Best sensitivity for very-deep homology |
| Domain assignment | hmmscan against Pfam-A | Curated, calibrated thresholds | Standard practice |
| Metagenomic protein clustering | MMseqs2 `easy-cluster` | Scales to >1B sequences | Production-grade |
| ORF search vs metagenome | DIAMOND `blastx --frameshift` | Frameshift-aware; long reads | Best for noisy long reads |
| Structure-aware homology (no AF2 prediction available) | Foldseek + ProstT5 | Predicts 3Di alphabet from sequence via PLM | Skip the AF2 step |

## Foldseek: the 2024 revolution

Foldseek (van Kempen, Kim, Tumescheit et al. 2024 *Nat Biotechnol* 42:243) searches protein structures by representing each residue's local geometry as a 21-letter "3Di" alphabet, then running BLAST-style alignment in this alphabet. Two consequences:
- **4-5 orders of magnitude faster than DALI** (the previous gold-standard structure aligner).
- **Finds homologs that sequence methods cannot**: when sequence has diverged past detection but structure is preserved, Foldseek recovers the homology. Reported sensitivity vs traditional structure search is comparable; sensitivity vs sequence methods is dramatically higher at low identity.

Two access modes:
1. **Have a structure** (PDB or AlphaFold): `foldseek easy-search query.pdb db_dir result.m8 tmp_dir`
2. **Sequence only, no structure**: use **ProstT5** (Heinzinger et al. 2024) to embed sequence to 3Di alphabet directly, skipping AF2 entirely: `foldseek databases ProstT5 prostt5_db tmp` then `foldseek easy-search seq.fa db result.m8 tmp --prostt5-model prostt5_db`

The major prebuilt Foldseek databases (AlphaFoldDB, PDB100, ESMAtlas) are downloadable via `foldseek databases`.

## PSI-BLAST: still useful, but watch the drift

PSI-BLAST (Altschul et al. 1997 *Nucleic Acids Res* 25:3389) builds a PSSM iteratively: each iteration includes hits below `-inclusion_ethresh` (default 0.005) in the next PSSM. Convergence is when no new hits cross the threshold. **Stopping at convergence is often the wrong call** -- iterations 2-3 are usually optimal; iterations 4+ frequently drift into paralog inclusion, contaminating the PSSM.

| Parameter | Default | Postdoc tuning |
|---|---|---|
| `-num_iterations` | 1 | 3 for most workflows; >3 risks drift |
| `-inclusion_ethresh` | 0.005 | 0.002 if specificity matters (Altschul 1997 recommendation) |
| `-evalue` | 10 | 0.01 for reporting cutoff |
| `-num_threads` | 1 | 8 for large DBs |

PSI-BLAST is also **non-deterministic** in detail: different input order or DB version can produce different PSSMs. For reproducibility, save the PSSM (`-out_pssm pssm.asn`) and re-use with `-in_pssm`.

## HMMER 3 (hmmsearch, jackhmmer)

HMMER 3 (Eddy 2011 *PLoS Comput Biol* 7:e1002195) is profile HMM search. Two main workflows:
- **`hmmsearch profile.hmm seqdb`**: search a database with a known HMM (Pfam, custom).
- **`jackhmmer query.fa seqdb`**: iterative search like PSI-BLAST but with full HMM math. Typically higher sensitivity than PSI-BLAST at the same number of iterations.

For domain assignment, `hmmscan query.fa Pfam-A.hmm` is the canonical pipeline. Pfam HMMs come with calibrated gathering thresholds (`-gathering`) -- use them instead of arbitrary E-value cutoffs.

```bash
# Build domain database once
wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
gunzip Pfam-A.hmm.gz
hmmpress Pfam-A.hmm

# Annotate query against Pfam-A with gathering threshold (calibrated cutoff)
hmmscan --cut_ga --domtblout query.domtbl Pfam-A.hmm query.fa
```

## HHblits / HHsearch (HH-suite3)

HHblits (Remmert et al. 2012 *Nat Methods* 9:173; HH-suite3: Steinegger et al. 2019 *BMC Bioinformatics* 20:473) is profile-profile alignment. The most sensitive method when both query and target have an HMM representation. Standard pipeline:
1. Build query MSA with `hhblits` against UniRef30 (or HHblits' default DB).
2. Convert MSA to query HMM.
3. Search against a profile DB (PDB70 for structure, Pfam-A for domains) with `hhsearch`.

Output is in HHM format. For very deep homology (the structural twilight zone), HHsearch vs PDB70 is still the gold standard.

## MMseqs2 (the modern protein search workhorse)

MMseqs2 (Steinegger & Soding 2017 *Nat Biotechnol* 35:1026) replaces BLAST in nearly all large-scale workflows.

Key advantages:
- **Speed**: 400-10,000x faster than blastp at similar sensitivity.
- **Sensitivity**: At `-s 7.5`, matches HMMER sensitivity (but on raw sequence, not profiles).
- **Iterative profile search**: `mmseqs search --num-iterations 3` matches PSI-BLAST behavior at much higher speed.

```bash
# Build target DB
mmseqs createdb target.fasta targetDB
mmseqs createindex targetDB tmp

# Sensitive search
mmseqs easy-search query.fa targetDB results.m8 tmp -s 7.5 --num-iterations 3

# All-vs-all clustering at 50% sequence identity
mmseqs easy-cluster all_proteins.fa cluster tmp --min-seq-id 0.5 -c 0.8
```

The `-s` parameter trades sensitivity for speed: 1.0 (fast), 4.0 (default), 7.5 (HMMER-like sensitivity).

## DIAMOND (the modern blastp replacement)

DIAMOND (Buchfink et al. 2015 *Nat Methods* 12:59; v2: Buchfink et al. 2021 *Nat Methods* 18:366) is the de-facto replacement for blastp on large-scale workflows.

| Feature | DIAMOND v2 | blastp |
|---|---|---|
| Speed | 100-10,000x faster | baseline |
| Sensitivity (default) | ~95% of blastp | baseline |
| `--ultra-sens