Ask BLAST to find homologs of a protein and it does well down to roughly 30 percent identity. Below that, alignment scores decay until true relatives are indistinguishable from noise, even though the structural and functional relationship is real. Profile HMMs solve precisely this problem: instead of comparing one sequence to another, they compare a sequence to a statistical model of a whole protein family.
HMMER is the tool that implements this. It is the engine behind Pfam, InterPro and most domain annotation you have ever read, and once you understand what it is doing, the E-values stop looking arbitrary.
Table of Contents
What a profile HMM actually is
A profile HMM is a position-specific model built from a multiple sequence alignment. At each column of the alignment it stores which residues occur and how often, and how often gaps appear. Match states emit residues, insert states allow extra residues between columns, and delete states skip a column entirely.
| Model feature | What it captures | Consequence for search |
|---|---|---|
| Match emission probabilities | Conserved residue preferences per column | Finds relatives that keep the same chemistry |
| Insert states | Tolerated length variation | Accepts indels inside a domain |
| Delete states | Tolerated missing columns | Handles partial domains and fragmentation |
| Transition probabilities | How often gaps follow gaps | Penalises implausible gap patterns |
| Log-odds scoring | Model versus a null background | Gives a score comparable across queries |
Because conservation is weighted per column, a profile HMM can find a relative that has diverged at 20 percent identity but keeps the few positions that define the family’s fold. That is information a pairwise alignment cannot use.
Install and get the databases
conda create -n hmmer -c bioconda -c conda-forge hmmer
conda activate hmmer
hmmscan -h | head -3 # version check
# Pfam-A is the standard domain database
wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
gunzip Pfam-A.hmm.gz
# Index it — hmmscan needs both, and will fail cryptically without them
hmmpress Pfam-A.hmmThat hmmpress step creates four indexed files alongside the model file. Forgetting it is the most common first error in HMMER; the message you get mentions a missing .h3m file and nothing else.
Annotating a protein with domains
# One protein, all Pfam domains
hmmscan --domtblout protein.domtbl \
--cpu 4 \
Pfam-A.hmm query.faa
# Human-readable summary
hmmscan --tblout protein.tbl Pfam-A.hmm query.faa
grep -v '^#' protein.tbl | headThe domain table is where the useful information lives. Columns you will keep coming back to: the domain name, the accession, the full-sequence E-value, the domain E-value, the alignment start and end on the query, and the hmm coordinates.
# domtblout columns worth extracting
# 1 target name 4 query name 7 full E-value
# 13 domain i-Evalue 18 hmm from 19 hmm to
# 20 ali from 21 ali to 22 env from 23 env to
awk '!/^#/ {printf "%-22s %-10s fullE=%-9s domE=%-9s env=%s-%s\n", $1, $2, $7, $13, $22, $23}' \
protein.domtbl | sort -k3,3g | head -10Note the difference between the full-sequence E-value and the domain E-value. The first asks how surprising it is that this model matches anywhere in the protein; the second asks how surprising it is that it matches this specific region. A protein with ten kinase domains will show a modest full-sequence E-value and ten very strong domain E-values. Filter on the domain value when you care about a region, and on the full-sequence value when you care about the protein.
The reverse direction: which family is this?
The direction of the search matters and beginners mix it up. hmmscan answers which models match my sequence. hmmsearch answers which sequences match my model. They are not interchangeable.
| Command | Query | Database | Question answered |
|---|---|---|---|
hmmscan | One protein | Many HMMs (Pfam) | What domains does this protein contain? |
hmmsearch | One HMM | Many proteins | Which proteins contain this domain? |
phmmer | One sequence | Protein sequences | Sequence similarity search, HMM-based |
jackhmmer | One sequence | Protein sequences | Iterative search, like PSI-BLAST |
# Which proteins in a proteome carry a kinase domain?
hmmsearch --domtblout hits.domtbl \
--cpu 8 \
PF00069.hmm proteome.faaBuilding your own profile HMM
You build a custom HMM when Pfam does not cover your family: a novel domain, a viral protein, or a gene family in a non-model organism. The quality of the resulting model depends almost entirely on the quality of the input alignment.
Step one: a trustworthy alignment
# Start from homologs you trust, then align them
mafft --auto family.faa > family.aln.faa
# Trim badly aligned columns — this matters more than anything else
trimal -in family.aln.faa -out family.trim.faa -automated1Include diversity deliberately. A model trained on ten near-identical sequences will match those ten perfectly and nothing else; that is overfitting, expressed as a search result. Aim for representatives spanning the family’s range, and keep the alignment to the homologous region rather than full-length proteins with variable flanks.
Step two: build and calibrate
hmmbuild family.hmm family.trim.faa
# Calibrate for E-values; without this hmmsearch prints no E-values at all
hmmpress family.hmm # only needed for hmmscan
hmmsearch --tblout out.tbl family.hmm query_proteome.faaIf you plan to use the model with hmmscan — scanning many proteins against your model set — you must run hmmpress. If you only use hmmsearch, pressing is unnecessary but harmless.
Step three: validate before you trust it
Three checks separate a usable model from a decorative one.
- Search a proteome that should not contain the family and confirm you get no strong hits.
- Search a proteome that should contain it and confirm recall matches what you know biologically.
- Re-search the training sequences and confirm they score well — but note that this proves nothing on its own, because you trained on them.
Reading E-values honestly
An E-value is the expected number of hits this good or better that would arise by chance in a database of this size. It is a function of the database, not just the sequence: the same protein against Pfam and against a ten-million-protein metagenome database gets very different E-values for identical evidence.
# The database size is printed in the header — read it
grep '^#' out.tbl | head -5Consequences that follow directly. First, a threshold such as 1e-5 is a convention, not a law; report the value you used. Second, comparing E-values across searches against different databases is meaningless. Third, in enormous databases even true homologs can drift above a strict threshold, so coverage and domain architecture often matter more than the E-value alone.
Coverage and bias correction
Two filters do most of the real work in practice: coverage of the model, and the composition bias filter.
hmmsearch --domtblout hits.domtbl \
--cut_ga \ # use Pfam's curated gathering thresholds
--max \
Pfam-A.hmm proteome.faaFor Pfam models, --cut_ga is the right choice. The gathering threshold is curated per family by Pfam’s curators, and it is more defensible than a blanket E-value cutoff because it accounts for how conserved each family actually is. Low-complexity proteins — collagens, proline-rich regions, coiled coils — produce spurious hits, and HMMER’s bias composition filter handles most of them; if you see a suspicious hit in a repetitive protein, look at the bias correction before believing it.
The failures you will actually hit
| Symptom | Cause | Fix |
|---|---|---|
| No E-values printed | Model not calibrated | Rebuild with standard pipeline, do not hand-edit the HMM |
hmmscan fails on missing index | Forgot hmmpress | Run hmmpress Pfam-A.hmm |
| Everything hits weakly | Database too small or model overfitted | Broaden training set; use a larger database |
| Many spurious hits | Low-complexity query | Rely on bias filter; inspect alignment, not just score |
| Custom model matches nothing | Alignment truncated or misaligned | Re-align with MAFFT and trim with trimAl |
What to report
For any domain annotation you publish, record the HMMER version, the database release, the threshold used (--cut_ga or an explicit E-value), and whether you searched proteins or models. A reader can then reproduce the annotation exactly, which is the difference between a claim and a finding.
Domain annotation is one layer of a larger argument about what a protein does. If you want to build the whole picture — from domain architecture to family membership to functional inference — Ampersand Academy covers the reasoning step by step.
Practice
- Scan a bacterial proteome against Pfam and count how many proteins have at least one domain at the gathering threshold.
- Compare the full-sequence and domain E-values for a multi-domain protein and explain the difference in one sentence.
- Build an HMM from a five-member family and search a proteome that should not contain it.
- Re-run the same scan with a plain E-value cutoff instead of
--cut_gaand quantify how many hits change.
Related tutorials
- Multiple Sequence Alignment: Beginner’s Guide to MSA — the alignment a profile HMM is built from
- Sequence Alignment Algorithms Explained: Global, Local, and BLAST — why profile models detect relationships pairwise scoring cannot
- Protein Structure Prediction Explained: From Amino Acid Sequence to AlphaFold — taking a domain annotation through to a predicted structure
What is the difference between BLAST and HMMER?
BLAST compares one sequence to another and scores pairwise similarity. HMMER compares a sequence to a statistical model built from a family alignment, which lets it detect relatives that have diverged past the point where pairwise scores become unreliable — roughly below 30 percent identity.
When should I use hmmscan and when hmmsearch?
Use hmmscan when you have one or a few proteins and want to know which families they belong to; it scans a sequence against a database of models. Use hmmsearch when you have a model and want to find every protein containing it, which is the usual direction for large-scale annotation.
Why does my hmmsearch output have no E-values?
The model was not calibrated for E-value statistics, which happens when a profile HMM is hand-edited or produced outside the standard pipeline. Rebuild it with hmmbuild from a proper alignment; do not patch the model file.
What does –cut_ga do?
It applies Pfam’s curated gathering thresholds, which are per-family cutoffs chosen by curators to separate true family members from partial or spurious matches. For Pfam models this is generally more defensible than a single global E-value threshold.
How many sequences do I need to build a good profile HMM?
There is no hard minimum, but fewer than about ten diverse, well-aligned members usually produces an overfitted model. Diversity matters more than count: twenty sequences spanning the family outperform a hundred near-duplicates from one genus.
Can I use HMMER on nucleotide sequences?
Yes, with the nucleotide alphabet and appropriate options, but protein-level search is far more sensitive for coding regions because the twenty-letter alphabet retains signal that the four-letter alphabet loses. For non-coding RNA families, specialised covariance models such as Infernal are the better tool.
What E-value should I use as a cutoff?
For Pfam models use –cut_ga. Otherwise 1e-5 is a common convention for domain-level significance, but remember the E-value depends on database size, so state both the threshold and the database. For very large databases, combine an E-value with a coverage requirement.