doi: 10.1371/journal.pcbi.1002195.
Epub 2011 Oct 20.
Accelerated Profile HMM Searches
Affiliations
- PMID: 22039361
- PMCID: PMC3197634
- DOI: 10.1371/journal.pcbi.1002195
Item in Clipboard
Accelerated Profile HMM Searches
PLoS Comput Biol.
2011 Oct.
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.
Conflict of interest statement
The author has declared that no competing interests exist.
Figures
A: Profile HMM architecture used by HMMER3 , , . 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 . 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.
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 . 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.
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.
. 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.
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 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).
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 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).
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.
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.
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.
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” 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.
Similar articles
-
Hidden Markov model speed heuristic and iterative HMM search procedure.BMC Bioinformatics. 2010 Aug 18;11:431. doi: 10.1186/1471-2105-11-431. BMC Bioinformatics. 2010. PMID: 20718988 Free PMC article.
-
ParAlign: a parallel sequence alignment algorithm for rapid and sensitive database searches.Nucleic Acids Res. 2001 Apr 1;29(7):1647-52. doi: 10.1093/nar/29.7.1647. Nucleic Acids Res. 2001. PMID: 11266569 Free PMC article.
-
A new generation of homology search tools based on probabilistic inference.Genome Inform. 2009 Oct;23(1):205-11. Genome Inform. 2009. PMID: 20180275
-
Profile hidden Markov models.Bioinformatics. 1998;14(9):755-63. doi: 10.1093/bioinformatics/14.9.755. Bioinformatics. 1998. PMID: 9918945 Review.
-
A Review of Parallel Implementations for the Smith-Waterman Algorithm.Interdiscip Sci. 2022 Mar;14(1):1-14. doi: 10.1007/s12539-021-00473-0. Epub 2021 Sep 6. Interdiscip Sci. 2022. PMID: 34487327 Free PMC article. Review.
Cited by
-
Sticking together: independent evolution of biofilm formation in different species of staphylococci has occurred multiple times via different pathways.BMC Genomics. 2024 Aug 28;25(1):812. doi: 10.1186/s12864-024-10719-y. BMC Genomics. 2024. PMID: 39198733
-
Aspergillus brasiliensis E_15.1: A Novel Thermophilic Endophyte from a Volcanic Crater Unveiled through Comprehensive Genome-Wide, Phenotypic Analysis, and Plant Growth-Promoting Trails.J Fungi (Basel). 2024 Jul 25;10(8):517. doi: 10.3390/jof10080517. J Fungi (Basel). 2024. PMID: 39194843 Free PMC article.
-
Seqrutinator: scrutiny of large protein superfamily sequence datasets for the identification and elimination of non-functional homologues.Genome Biol. 2024 Aug 26;25(1):230. doi: 10.1186/s13059-024-03371-y. Genome Biol. 2024. PMID: 39187866 Free PMC article.
-
Structure and rational engineering of the PglX methyltransferase and specificity factor for BREX phage defence.Nat Commun. 2024 Aug 22;15(1):7236. doi: 10.1038/s41467-024-51629-7. Nat Commun. 2024. PMID: 39174540 Free PMC article.
-
Navigating the human-monkeypox virus interactome: HuPoxNET atlas reveals functional insights.Front Microbiol. 2024 Aug 2;15:1399555. doi: 10.3389/fmicb.2024.1399555. eCollection 2024. Front Microbiol. 2024. PMID: 39155985 Free PMC article.
References
-
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–410. - PubMed
-
- Krogh A, Brown M, Mian IS, Sjölander K, Haussler D. Hidden Markov models in computational biology: Applications to protein modeling. J Mol Biol. 1994;235:1501–1531. - PubMed
-
- Durbin R, Eddy SR, Krogh A, Mitchison GJ. Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge UK: Cambridge University Press; 1998.
Publication types
MeSH terms
Grants and funding
LinkOut - more resources
Full Text Sources
Other Literature Sources
Research Materials
Miscellaneous
