CoChat
Sign Up
Cited by 6,975Genomics and Phylogenetic StudiesAlgorithms and Data CompressionFractal and DNA sequence analysis

The DOI matches its Crossref record; 46 of the 46 references it cites were found in a registry or an index.

Accelerated Profile HMM Searches

Open access
Eddy, Sean R.·2011·PLoS Computational Biology
DOI 10.1371/journal.pcbi.1002195PMC PMC3197634
Open this paper in CoChat — free accountOpen PDF

Abstract

Profile hidden Markov models (profile HMMs) and probabilistic inference methods have made important contributions to the theory of sequence database homology search. However, practical use of profile HMM methods has been hindered by the computational expense of existing software implementations. Here I describe an acceleration heuristic for profile HMMs, the "multiple segment Viterbi" (MSV) algorithm. The MSV algorithm computes an optimal sum of multiple ungapped local alignment segments using a striped vector-parallel approach previously described for fast Smith/Waterman alignment. MSV scores follow the same statistical distribution as gapped optimal local alignment scores, allowing rapid evaluation of significance of an MSV score and thus facilitating its use as a heuristic filter. I also describe a 20-fold acceleration of the standard profile HMM Forward/Backward algorithms using a method I call "sparse rescaling". These methods are assembled in a pipeline in which high-scoring MSV hits are passed on for reanalysis with the full HMM Forward/Backward algorithm. This accelerated pipeline is implemented in the freely available HMMER3 software package. Performance benchmarks show that the use of the heuristic MSV filter sacrifices negligible sensitivity compared to unaccelerated profile HMM searches. HMMER3 is substantially more sensitive and 100- to 1000-fold faster than HMMER2. HMMER3 is now about as fast as BLAST for protein searches.

Overview

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Integrity signals

Registry

Crossref · 10.1371/journal.pcbi.1002195

Cites

46 of 46 references verified0 citation errors · 0 not checkable

Statistics

Figures

Figure 1 The MSV profile. A: Profile HMM architecture used by HMMER3 [4], [5], [31]. Regions homologously aligned to the query are represented by a linear core model consisting of consensus positions (in this example, ), each consisting of a match, a delete, and an insert state (shown as boxes marked M, circles marked D, and diamonds marked I), connected by state transition probabilities (arrows). Match states carry position-specific emission probabilities for scoring residues at each consensus position. Insert states emit residues with emission probabilities identical to a background distribution. Additional flanking states (marked N, C, and J) emit zero or more residues from the background distribution, modeling nonhomologous regions preceding, following, or joining homologous regions aligned to the core model. Start (S), begin (B), end (E) and termination (T) states do not emit. B: The MSV profile is formed by implicitly treating all match-match transition probabilities as 1.0. This corresponds to the virtual removal of the delete and insert states. The rest of the profile parameterization stays the same. This model generates sequences containing one or more ungapped local alignment segments. Note that both models appear to be improperly normalized; for example, each match state in the MSV model has probability 1.0 local exit transition (orange arrows) in addition to the probability 1.0 match-match transition. This is because of a trick used to establish a uniform local fragment length distribution, in which these profiles are collapsed representations of a much larger (and properly normalized) “implicit probability model”, as explained in [31]. C: An example of what an alignment of a larger MSV profile (of length ) to a target sequence (of length ) might look like, as a path through a dynamic programming (DP) matrix. Here, the model identifies two high-scoring ungapped alignment segments (black dots, indicating residues aligned to profile match states), and assigns all other residues to N, J, and C states in the model (orange dots; unfilled indicates a “mute” nonemitting state or state transition). Note that the ungapped diagonals are not enforced to be consistent with a single gapped alignment.
Figure 1 The MSV profile. A: Profile HMM architecture used by HMMER3 [4], [5], [31]. Regions homologously aligned to the query are represented by a linear core model consisting of consensus positions (in this example, ), each consisting of a match, a delete, and an insert state (shown as boxes marked M, circles marked D, and diamonds marked I), connected by state transition probabilities (arrows). Match states carry position-specific emission probabilities for scoring residues at each consensus position. Insert states emit residues with emission probabilities identical to a background distribution. Additional flanking states (marked N, C, and J) emit zero or more residues from the background distribution, modeling nonhomologous regions preceding, following, or joining homologous regions aligned to the core model. Start (S), begin (B), end (E) and termination (T) states do not emit. B: The MSV profile is formed by implicitly treating all match-match transition probabilities as 1.0. This corresponds to the virtual removal of the delete and insert states. The rest of the profile parameterization stays the same. This model generates sequences containing one or more ungapped local alignment segments. Note that both models appear to be improperly normalized; for example, each match state in the MSV model has probability 1.0 local exit transition (orange arrows) in addition to the probability 1.0 match-match transition. This is because of a trick used to establish a uniform local fragment length distribution, in which these profiles are collapsed representations of a much larger (and properly normalized) “implicit probability model”, as explained in [31]. C: An example of what an alignment of a larger MSV profile (of length ) to a target sequence (of length ) might look like, as a path through a dynamic programming (DP) matrix. Here, the model identifies two high-scoring ungapped alignment segments (black dots, indicating residues aligned to profile match states), and assigns all other residues to N, J, and C states in the model (orange dots; unfilled indicates a “mute” nonemitting state or state transition). Note that the ungapped diagonals are not enforced to be consistent with a single gapped alignment.CC BY · PubMed Central
Figure 2 Illustration of striped indexing for SIMD vector calculations. The top row (magenta outline) shows one row of the dynamic programming lattice for a model of length . Assuming an example of vectors containing cells each, the 14 cells are contained in vectors numbered . (Two unused cells, marked x, are set to a sentinel value.) In the dynamic programming recursion, when we calculate each new cell in a new row , we access the value in cell in the previous row . With striped indexing, vector contains exactly the four cells needed to calculate the four cells in a new vector on a new row of the dynamic programming matrix (turquoise outline). For example, when we calculate cells in vector , we access the previous row's vector which contains the cells we need in the order we need them, (dashed lines and box). If instead we indexed cells into vectors in the obvious way, in linear order ( in vector and so on), there is no such correspondence of with four 's, and each calculation of a new vector would require expensive meddling with the order of cells in the previous row's vectors. With striped indexing, only one shift operation is needed per row, outside the innermost loop: the last vector on each finished row is rightshifted (mpv, in grey with red cell indices) and used to initialize the next row calculation.
Figure 2 Illustration of striped indexing for SIMD vector calculations. The top row (magenta outline) shows one row of the dynamic programming lattice for a model of length . Assuming an example of vectors containing cells each, the 14 cells are contained in vectors numbered . (Two unused cells, marked x, are set to a sentinel value.) In the dynamic programming recursion, when we calculate each new cell in a new row , we access the value in cell in the previous row . With striped indexing, vector contains exactly the four cells needed to calculate the four cells in a new vector on a new row of the dynamic programming matrix (turquoise outline). For example, when we calculate cells in vector , we access the previous row's vector which contains the cells we need in the order we need them, (dashed lines and box). If instead we indexed cells into vectors in the obvious way, in linear order ( in vector and so on), there is no such correspondence of with four 's, and each calculation of a new vector would require expensive meddling with the order of cells in the previous row's vectors. With striped indexing, only one shift operation is needed per row, outside the innermost loop: the last vector on each finished row is rightshifted (mpv, in grey with red cell indices) and used to initialize the next row calculation.CC BY · PubMed Central
Figure 3 MSV scores follow a predictable distribution. A: example MSV score distributions for a typical Pfam model, CNP1, on random i.i.d. sequences of varying lengths from 25 to 25,600, with the shortest, typical, and longest lengths highlighted as red, black, and blue lines, respectively. The predicted distribution, following the procedure of [31] including an edge correction on the slope , is shown in orange (though largely obscured by the data lines right on top of it). B: Histogram of maximum likelihood values obtained from score distributions of 11,912 Pfam models, showing that most are tolerably close to the conjectured , albeit with more dispersion for default entropy-weighted models (black line) than high relative entropy models without entropy-weighting (gray line). C: The observed fraction of nonhomologous sequences that pass the filter at a P-value of 0.02 should be 0.02. Histograms of the actual filter fraction for 11,912 different Pfam 24 models are shown, for a range of random sequence lengths from 25 to 25,600, for both default models (black lines) and high relative entropy models with no entropy weighting (gray lines).
Figure 3 MSV scores follow a predictable distribution. A: example MSV score distributions for a typical Pfam model, CNP1, on random i.i.d. sequences of varying lengths from 25 to 25,600, with the shortest, typical, and longest lengths highlighted as red, black, and blue lines, respectively. The predicted distribution, following the procedure of [31] including an edge correction on the slope , is shown in orange (though largely obscured by the data lines right on top of it). B: Histogram of maximum likelihood values obtained from score distributions of 11,912 Pfam models, showing that most are tolerably close to the conjectured , albeit with more dispersion for default entropy-weighted models (black line) than high relative entropy models without entropy-weighting (gray line). C: The observed fraction of nonhomologous sequences that pass the filter at a P-value of 0.02 should be 0.02. Histograms of the actual filter fraction for 11,912 different Pfam 24 models are shown, for a range of random sequence lengths from 25 to 25,600, for both default models (black lines) and high relative entropy models with no entropy weighting (gray lines).CC BY · PubMed Central
Figure 4 The HMMER3 acceleration pipeline. Representative calculation speeds are shown in red, in units of millions of dynamic programming cells per second (Mc/s). Default P-value thresholds for MSV, Viterbi, and Forward filtering steps are shown in orange. The bias filter and the domain definition steps are not described in detail in this manuscript, and are shown in gray.
Figure 4 The HMMER3 acceleration pipeline. Representative calculation speeds are shown in red, in units of millions of dynamic programming cells per second (Mc/s). Default P-value thresholds for MSV, Viterbi, and Forward filtering steps are shown in orange. The bias filter and the domain definition steps are not described in detail in this manuscript, and are shown in gray.CC BY · PubMed Central
Figure 5 Speed benchmarks. Each point represents a speed measurement for one search with one query against target sequences ( for the slow HMMER2 and SAM programs, for FASTA and SSEARCH), on a single CPU core (see Methods for more details). Both axes are logarithmic, for speed in millions of dynamic programming cells per second (Mc/s) on the y-axis and query length in residues on the x-axis. Panel A shows “typical best performance” speed measurements for several different programs including HMMER3, for 76 queries of varying consensus lengths, chosen from Pfam 24, for searches of randomized (shuffled) target sequences. Panel B shows a wider range of more realistic speed measurements for all 11,912 profiles in Pfam 24, on searches of real target protein sequences from UniProt TrEMBL.
Figure 5 Speed benchmarks. Each point represents a speed measurement for one search with one query against target sequences ( for the slow HMMER2 and SAM programs, for FASTA and SSEARCH), on a single CPU core (see Methods for more details). Both axes are logarithmic, for speed in millions of dynamic programming cells per second (Mc/s) on the y-axis and query length in residues on the x-axis. Panel A shows “typical best performance” speed measurements for several different programs including HMMER3, for 76 queries of varying consensus lengths, chosen from Pfam 24, for searches of randomized (shuffled) target sequences. Panel B shows a wider range of more realistic speed measurements for all 11,912 profiles in Pfam 24, on searches of real target protein sequences from UniProt TrEMBL.CC BY · PubMed Central
Figure 6 Benchmark of search sensitivity and specificity. For different programs, searches are performed either by constructing a single profile from the query alignment (HMMER3, HMMER2, SAM, PSI-BLAST), or by using “family pairwise search” [41] in which each individual sequence is used as a query and the best E-value per target sequence is recorded (BLASTP, SSEARCH, FASTA). In each benchmark, true positive subsequences have been selected to be no more than 25% identical to any sequence in the query alignment (see Methods). Panel A shows results where nonhomologous sequence has been synthesized by a simple random model, and each true positive sequence contains a single embedded homologous subsequence (a total of 2,141 query multiple alignments, 11,547 true positive sequences, and 200,000 decoys). Panel B shows results where nonhomologous sequence is synthesized by shuffling randomly chosen subsequences from UniProt, and each true positive contains two embedded homologous subsequences (a total of 2,141 query alignments, 24,040 true positive sequences, and 200,000 decoys). The Y-axis is the fraction of true positives detected with an E-value better than the number of false positives per query specified on the X-axis.
Figure 6 Benchmark of search sensitivity and specificity. For different programs, searches are performed either by constructing a single profile from the query alignment (HMMER3, HMMER2, SAM, PSI-BLAST), or by using “family pairwise search” [41] in which each individual sequence is used as a query and the best E-value per target sequence is recorded (BLASTP, SSEARCH, FASTA). In each benchmark, true positive subsequences have been selected to be no more than 25% identical to any sequence in the query alignment (see Methods). Panel A shows results where nonhomologous sequence has been synthesized by a simple random model, and each true positive sequence contains a single embedded homologous subsequence (a total of 2,141 query multiple alignments, 11,547 true positive sequences, and 200,000 decoys). Panel B shows results where nonhomologous sequence is synthesized by shuffling randomly chosen subsequences from UniProt, and each true positive contains two embedded homologous subsequences (a total of 2,141 query alignments, 24,040 true positive sequences, and 200,000 decoys). The Y-axis is the fraction of true positives detected with an E-value better than the number of false positives per query specified on the X-axis.CC BY · PubMed Central

Verification

Claims the paper makes

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Claims other papers make about this one

Each is checked against this paper’s own text: does it say what the citing paper says it does?

Self-check

Percentages, sums and reported differences re-computed from the paper’s own numbers.

The self-check needs the paper’s full text.

Citations

Cites

Show all 46 references

Research works printed but not linked: 1

  1. S Johnson, 2006, Remote Protein Homology Detection Using Hidden Markov Models [Ph.D. thesis]

Cited by

Sorted by citations

Related

See all

Metadata

Authors
Eddy, S. R.
Affiliations
Howard Hughes Medical Institute (US)
Corresponding author
Eddy, Sean R.
ORCID
Eddy, Sean R.
Venue
PLoS Computational Biology
Publisher
Public Library of Science (PLoS)
Published
2011-10-20
Volume
7
Issue
10
Pages
e1002195
ISSN
  • 1553-734X (ISSN-L)
  • 1553-7358
DOI
10.1371/journal.pcbi.1002195
PMC
PMC3197634
OpenAlex
W2138122982
Cited by
6,975
References
46
Citations by year
2026: 791 · 2025: 804 · 2024: 768 · 2023: 693 · 2022: 671 · 2021: 706 · 2020: 496 · 2019: 411 · 2018: 386 · 2017: 312 · 2016: 271 · 2015: 258 · 2014: 206 · 2013: 145 · 2012: 52
FWCI
52.51
Type
journal article
Language
en
Open access
Licence
Locations
  • PLoS Computational Biology (publishedVersion, cc-by)
  • Europe PMC (PubMed Central) (publishedVersion, cc-by)
  • PubMed Central (publishedVersion, cc-by)
  • PubMed (publishedVersion)
  • RePEc: Research Papers in Economics (submittedVersion)
  • DOAJ (DOAJ: Directory of Open Access Journals) (submittedVersion, cc-by-sa)
Funders
  • Howard Hughes Medical Institute
  • Hebrew University of Jerusalem
Topics
  • Life Sciences › Biochemistry, Genetics and Molecular Biology › Molecular Biology › Genomics and Phylogenetic Studies
  • Physical Sciences › Mathematics › Theoretical Computer Science › Algorithms and Data Compression
  • Life Sciences › Biochemistry, Genetics and Molecular Biology › Molecular Biology › Fractal and DNA sequence analysis
Keywords
profile hidden Markov models, sequence homology search, protein sequence alignment, computational acceleration, heuristic filtering, forward-backward algorithm
MeSH terms
Algorithms, Markov Chains, Software, Information Storage and Retrieval, Sequence Alignment, Computational Biology, Databases, Protein

About these checks

Registry, retraction, Cites and statistics results are computed, not judged, with one exception: a printed reference that no index lists is read once by a model before it counts as a citation error. Claims are extracted by a model and kept only when their quote appears verbatim in the text; each cited source is then checked, with a model’s help, for whether it says what the claim relies on. Those verdicts can be wrong — read the claim’s quote and the cited source before relying on one.

Statistics and arithmetic are recomputed, not judged.

Confidence bands are shown only for fields calibrated on our annotated test set. None is calibrated yet, so no band is shown.

What counts as a claim ›How we measure confidence ›How we check this ›

Are you an author, or do you see a problem with this page? Request a review