PyHMMER

PyHMMER is a Python module implemented using Cython that provides direct bindings to HMMER3. Unlike CLI wrappers, it interacts directly with HMMER internals, offering several advantages:

• Batteries-included: Can be installed as a pip or conda dependency without requiring external HMMER binaries.
• Flexible: Allows creating pyhmmer.easel.Sequence and pyhmmer.plan7.HMM objects via the API or loading them from files.
• Practical: Nested results can be retrieved as pyhmmer.plan7.TopHits objects for file writing or Python-based analysis.
• Fast: Supports parallel execution of hmmsearch using an efficient threading model that can outperform the original HMMER in certain use cases.
• Shareable: Supports distributing and loading pyhmmer.plan7.HMM objects directly from within Python packages.
• Transparent: Provides access to pyhmmer.plan7.HMM internals, allowing inspection and manual editing of transition or emission scores.

The pyhmmer.easel module provides a high-performance framework for handling biological data, including sequences, alignments, and linear algebra operations. It is organized into several functional areas:

• Biological Alphabets: Defines Alphabet types such as DNA, RNA, and AA (Amino Acids).
• Sequences: Core sequence representations including Sequence, TextSequence, and DigitalSequence.
• Sequence Blocks: Specialized containers for sequence segments like SequenceBlock, TextSequenceBlock, and DigitalSequenceBlock.
• Alignments: Support for Multiple Sequence Alignments (MSA) via MSA, TextMSA, and DigitalMSA.
• Parsers: Tools for reading biological files, specifically SequenceFile and MSAFile.
• Linear Algebra: High-performance math primitives including Vector and Matrix types with various precision levels (D for double, F for float, U8 for unsigned 8-bit).
• Data Structures: Low-level utilities like Bitfield and KeyHash.

When performing HMM searches or scans, hits are reported with three primary statistics. Understanding these is essential for filtering results:

• Hit Bit-score (Hit.score): The score of the alignment of a sequence to an HMM, expressed in bits. This is calculated after the sequence passes all pipeline filters. If a domain definition workflow is run, corrected bit-scores are calculated for both the sequence level and individual domains.
• Hit P-value: The probability $P(X \ge \text{score})$ that a score at least as high as the observed score would occur by chance. It is computed using the HMM's calibration parameters ($\tau$ and $\lambda$).
• Hit E-value (Hit.evalue): The expected number of hits with a score at least as high as the observed score. It is computed by multiplying the Hit.pvalue by the number of testing hypotheses $Z$ (usually the number of targets) using Bonferroni's method to correct for multiple testing.

When the 'domain definition workflow' identifies multiple HMM matches to the same sequence, PyHMMER reports statistics at the domain level:

• Domain Bit-score: The bit score computed based on the local alignment of the specific domain only.
• Domain P-value: Computed identically to the hit P-value, but using the domain bit-score and the HMM's E-value parameters.
• Domain E-values: Because only hits passing pipeline filters (defined by Pipeline parameters $F_1$, $F_2$, and $F_3$) reach the domain workflow, HMMER provides two types of E-value corrections:
  • Conditional E-value (Domain.c_evalue): Takes into account only the number of domains passing the first pipeline stages (denoted as domZ). This is generally more statistically accurate but harder to replicate without knowing domZ a-priori.
  • Independent E-value (Domain.i_evalue): Takes into account all sequence/HMM comparisons independently (denoted as $Z$).

PyHMMER uses a different parallelization strategy than the standard HMMER C binaries to achieve better CPU utilization, especially when searching multiple HMMs against multiple targets.

HMMER (Standard Binary) Strategy

• Serial mode: Reads HMMs and sequences one by one. Low memory usage but becomes I/O-bound because sequence files must be re-read for every HMM.
• Threaded mode: Clones HMM profiles for worker threads and reads sequences in fixed blocks (e.g., 1000 sequences). This can suffer from load imbalance where some threads idle while others process large blocks, and it only achieves true parallelism if the number of targets is very large ($1000 imes n_{cpus}$).

PyHMMER Strategy

PyHMMER pre-fetches target sequences into memory and passes them by reference to worker threads. Each worker receives an HMM from the main thread and processes the entire set of targets for that HMM. This allows the pipeline to run multiple HMMs in parallel simultaneously, preventing threads from idling and ensuring better scaling across available CPUs.

Note: This optimized parallelization strategy is most effective when performing multiple queries (e.g., annotating a large protein dataset with a subset of the Pfam HMM library).