Bio Remote Homology

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

# 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

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.

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.

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

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

# Build DIAMOND DB
diamond makedb --in nr.fa -d nr

# Sensitive search
diamond blastp -d nr -q query.fa -o results.tsv \
        --more-sensitive -e 1e-10 -p 16 \
        --outfmt 6 qseqid sseqid pident length qcovhsp evalue bitscore stitle

# Long-read frameshift-aware (for nanopore/PacBio metagenomics)
diamond blastx -d nr -q longreads.fa --frameshift 15 -o reads.tsv --outfmt 6

For protein remote homology in 2026, DIAMOND --ultra-sensitive or MMseqs2 -s 7.5 should be the default before reaching for BLAST.

Iterative HMMER (jackhmmer)

# 3 iterations, save checkpoint HMM at each iteration
jackhmmer -N 3 --chkhmm iter.hmm --tblout hits.tbl query.fa uniref90.fa

# After convergence (no new hits below threshold) use the final HMM for downstream searches
hmmsearch iter-3.hmm target.fa > hits.txt

Code patterns

Foldseek search against AlphaFoldDB

Goal: Find structural homologs of a protein structure (or sequence via ProstT5) in AlphaFold's predicted structure database.

Approach: Download AlphaFoldDB structure DB (or ProstT5 for sequence-only); search with foldseek easy-search; parse m8 tabular output.

Reference (Foldseek 9+):

#!/bin/bash
# Reference: foldseek 9+ | Verify API if version differs

mkdir -p foldseek_dbs tmp
# Download the AlphaFoldDB Swiss-Prot subset (~few GB; full AFDB is much larger)
foldseek databases Alphafold/Swiss-Prot afdb_sp foldseek_dbs/tmp

# Structure-vs-structure search
foldseek easy-search query.pdb foldseek_dbs/afdb_sp results.m8 tmp \
         --format-output query,target,fident,alnlen,evalue,bits,prob,qtmscore,ttmscore

head -5 results.m8
# qtmscore/ttmscore are TM-score equivalents from local Foldseek alignment.
# Hits with prob > 0.9 are confidently structurally homologous.

Sequence-only Foldseek via ProstT5

foldseek databases ProstT5 prostt5_model tmp
foldseek easy-search query.fa foldseek_dbs/afdb_sp seq_results.m8 tmp \
         --prostt5-model prostt5_model --threads 8

This is the path to take when only a protein sequence is available -- ProstT5 (a protein language model) predicts the 3Di alphabet directly from sequence.

PSI-BLAST with saved PSSM

Goal: Build a position-specific scoring matrix iteratively, then re-use it for downstream searches.

Approach: 3 iterations against UniRef90 (or nr); save ASN.1 + ASCII PSSM; subsequent searches use -in_pssm.

Reference (NCBI BLAST+ 2.15+):

psiblast -query distant_protein.fa -db uniref90 \
         -num_iterations 3 \
         -inclusion_ethresh 0.002 \
         -evalue 0.01 \
         -num_threads 8 \
         -out_pssm distant.pssm.asn \
         -out_ascii_pssm distant.pssm.txt \
         -out psiblast_results.txt

# Reuse saved PSSM in subsequent searches against a different DB
psiblast -in_pssm distant.pssm.asn -db swissprot \
         -out swissprot_via_pssm.txt