🧠 HMM Fundamentals

Q: What is a Hidden Markov Model?

A Hidden Markov Model (HMM) is a statistical model where the system being modeled is assumed to be a Markov process with hidden (unobserved) states. Each hidden state emits observable symbols according to a probability distribution. HMMs are used for sequence analysis, gene prediction, speech recognition, and many other tasks.

Q: What is the difference between Viterbi, Forward-Backward, and Baum-Welch?

Viterbi finds the single most likely sequence of hidden states given observations. Forward-Backward computes the probability of being in each state at each position (posterior probabilities). Baum-Welch is an Expectation-Maximization algorithm that learns the model parameters (transition and emission probabilities) from unlabeled data.

Q: What is a profile HMM and what is HMMER used for?

A profile HMM is an HMM built from a multiple sequence alignment of a protein family. It captures position-specific conservation and insertion/deletion patterns. HMMER is the standard tool for building profile HMMs and searching databases for remote homologs (e.g., Pfam, InterPro).

Q: What are the three main problems of HMMs?

The three canonical HMM problems are: (1) Evaluation: compute the probability of an observation sequence given the model (Forward algorithm); (2) Decoding: find the most likely hidden state sequence (Viterbi algorithm); (3) Learning: estimate model parameters from data (Baum-Welch / EM algorithm).

Q: Which Python library should I use for HMMs?

Use hmmlearn for general-purpose discrete and Gaussian HMMs (scikit-learn-style API). Use pomegranate for flexible probabilistic models. For profile HMMs and biological sequence search, use HMMER (pyhmmer for Python bindings). For teaching, implement Viterbi/Forward-Backward from scratch.

🎯 What is a Hidden Markov Model?

A Hidden Markov Model (HMM) is a statistical model where the system being modeled is assumed to be a Markov process with hidden (unobserved) states. Each hidden state emits observable symbols according to a probability distribution. HMMs are used for sequence analysis, gene prediction, speech recognition, and many other tasks. They are a cornerstone of probabilistic sequence modeling in bioinformatics.

Why HMMs Matter

Sequence: Modeling. Profile: Protein families. Gene: Prediction. Speech: Recognition. NLP: Tagging. Finance: Regime switching. Foundation: Probabilistic.

📌 Hidden States📌 Markov Property📌 Emissions📌 Probabilistic

# Hidden Markov Models — the big picture
# Probabilistic models for sequences.

# CONCEPTUAL OVERVIEW
print("""
=== Hidden Markov Models (HMM) ===

What it is:
  - Statistical model
  - Hidden states
  - Observable emissions
  - Markov property

Why it's powerful:
  - Sequence modeling
  - Probabilistic
  - Handles variable length
  - Interpretable
  - Foundation for many tools

Components:
  - States (hidden)
  - Transitions
  - Emissions
  - Initial probabilities
  - Final probabilities

Types:
  - Discrete HMM
  - Gaussian HMM
  - Profile HMM
  - Pair HMM
  - Input-output HMM

Three problems:
  1. Evaluation: P(O | λ)
     → Forward algorithm
  2. Decoding: argmax Q
     → Viterbi algorithm
  3. Learning: argmax λ
     → Baum-Welch (EM)

Applications:
  - Gene prediction
  - Protein family profiles
  - Sequence alignment
  - Speech recognition
  - NLP tagging
  - Financial modeling
  - CpG island detection
  - Splice site prediction
""")

# WHY IT MATTERS
print("=== Why It Matters ===")
print("""
Biological questions:
  - Is this a gene?
  - What family is this protein?
  - Where are the splice sites?
  - What is the secondary structure?
  - Are there CpG islands?

HMMs answer:
  - State sequences
  - Probabilities
  - Likelihoods
  - Annotations
  - Classification
""")

# HISTORICAL CONTEXT
print("=== Historical Context ===")
print("""
1960s: Baum-Welch, Forward-Backward
1966: Viterbi algorithm
1970s: Speech recognition
1980s: HMMs in bioinformatics
1989: Profile HMMs (Gribskov)
1994: HMMER (Eddy)
1998: GENSCAN (Burge & Karlin)
2000s: Pfam, profile HMMs
2010s: Deep learning hybrids
2020s: Still foundational
""")

# CHALLENGES
print("=== Challenges ===")
print("""
1. LOCAL OPTIMA
   - Baum-Welch
   - Initialization
   - Multiple restarts

2. MODEL SELECTION
   - Number of states
   - Topology
   - Complexity

3. COMPUTATIONAL
   - Large sequences
   - Many states
   - Memory

4. INTERPRETATION
   - Hidden states
   - Biological meaning
   - Validation

5. OVERFITTING
   - Too many parameters
   - Small data
   - Regularization
""")

🎯 Why HMMs Matter

🎯 The Value of HMMs

HMMs are foundational to probabilistic sequence analysis, enabling gene prediction, protein family modeling, and much more.

Key Applications

Gene: Prediction. Protein: Families. Alignment: Profile. Speech: Recognition. NLP: Tagging. Finance: Regimes.

🎯 Gene prediction🎯 Protein families🎯 Alignment🎯 Speech

# Why HMMs Matter
print("=== 1. Gene Prediction ===")
print("""
Gene prediction:
  - Find genes in DNA
  - Exons, introns
  - Splice sites

HMM:
  - States: exon, intron, intergenic
  - Transitions: splice sites
  - Emissions: codons, bases

Tools:
  - GENSCAN
  - AUGUSTUS
  - Glimmer
  - GeneMark

Example:
  States: E (exon), I (intron), N (intergenic)
  Transitions: E→I (splice), I→E (splice)
""")

print("\n=== 2. Protein Family Profiles ===")
print("""
Protein families:
  - Profile HMM
  - Position-specific
  - Remote homology

HMM:
  - Match, insert, delete states
  - Position-specific emissions
  - Built from MSA

Tools:
  - HMMER
  - Pfam
  - InterPro
  - SMART

Example:
  - Kinase domain
  - Zinc finger
  - Immunoglobulin
""")

print("\n=== 3. Sequence Alignment ===")
print("""
Alignment:
  - Pair HMM
  - Probabilistic
  - Optimal

HMM:
  - Match, insert, delete
  - Pairwise
  - Posterior

Tools:
  - HMMER (hmmalign)
  - SAM
  - ProbCons

Example:
  - Protein alignment
  - DNA alignment
""")

print("\n=== 4. CpG Island Detection ===")
print("""
CpG islands:
  - High CpG content
  - Promoters
  - Methylation

HMM:
  - States: CpG+, CpG-
  - Emissions: A, C, G, T
  - Transitions

Tools:
  - CpGIE
  - CpGProD
  - HMM (custom)

Example:
  - CpG+ state: high CG
  - CpG- state: low CG
""")

print("\n=== 5. Splice Site Prediction ===")
print("""
Splice sites:
  - Donor: GT
  - Acceptor: AG
  - Branch point

HMM:
  - States: exon, intron
  - Emissions: sequence
  - Transitions: splice

Tools:
  - GENSCAN
  - AUGUSTUS
  - SplicePort

Example:
  - Donor: GTRAGT
  - Acceptor: YYYYYYYYYYNCAG
""")

print("\n=== 6. Secondary Structure ===")
print("""
Secondary structure:
  - α-helix, β-sheet
  - Prediction

HMM:
  - States: H, E, C
  - Emissions: amino acids
  - Transitions

Tools:
  - HMMSTR
  - SAM
  - PSIPRED (not HMM)

Example:
  - H: helix
  - E: strand
  - C: coil
""")

print("\n=== 7. Speech Recognition ===")
print("""
Speech:
  - Acoustic → words
  - HMM states: phonemes

HMM:
  - States: phonemes
  - Emissions: acoustic features
  - Transitions

Tools:
  - HTK
  - Kaldi
  - CMU Sphinx

Example:
  - Word "cat": K, AE, T
""")

print("\n=== 8. NLP Tagging ===")
print("""
NLP:
  - POS tagging
  - NER
  - Chunking

HMM:
  - States: POS tags
  - Emissions: words
  - Transitions

Tools:
  - NLTK
  - hmmlearn
  - pomegranate

Example:
  - "The cat sat" → DET NOUN VERB
""")

🔗 Markov Chains

🎯 The Foundation

A Markov chain is a stochastic process where the future depends only on the present, not the past. It's the foundation of HMMs.

Markov Chains

States: Observable. Transitions: Probabilities. Memoryless: Markov property. Stationary: Time-homogeneous. Matrix: Transition.

🔗 States🔗 Transitions🔗 Memoryless🔗 Matrix

# Markov Chains
print("=== Definition ===")
print(""
Markov chain:
  - States: S = {s1, s2, ..., sN}
  - Transitions: a_ij = P(q_t = s_j | q_{t-1} = s_i)
  - Markov property: P(q_t | q_{t-1}, ..., q_1) = P(q_t | q_{t-1})
  - Initial: π_i = P(q_1 = s_i)

Memoryless:
  - Future depends only on present
  - Not on past
  - Simplifies modeling
""")

print("\n=== Transition Matrix ===")
print(""
Transition matrix A:
  A = [a_ij]

Properties:
  - a_ij ≥ 0
  - Σ_j a_ij = 1 (each row sums to 1)

Example (weather):
  States: Sunny, Rainy
  A = [[0.8, 0.2],
       [0.3, 0.7]]

Meaning:
  P(Sunny | Sunny) = 0.8
  P(Rainy | Sunny) = 0.2
  P(Sunny | Rainy) = 0.3
  P(Rainy | Rainy) = 0.7
""")

print("\n=== Example: Weather ===")
print(""
States: Sunny (S), Rainy (R)
Initial: π = [0.5, 0.5]
Transitions:
  A = [[0.8, 0.2],
       [0.3, 0.7]]

Sequence:
  Day 1: Sunny (π_S = 0.5)
  Day 2: Sunny (0.5 * 0.8 = 0.4)
  Day 3: Rainy (0.4 * 0.2 = 0.08)

Probability:
  P(S, S, R) = 0.5 * 0.8 * 0.2 = 0.08
""")

print("\n=== Example: CpG Island ===")
print(""
States: CpG+ (high CG), CpG- (low CG)
Transitions:
  A = [[0.9, 0.1],
       [0.1, 0.9]]

Emissions (not yet HMM):
  CpG+: A=0.1, C=0.4, G=0.4, T=0.1
  CpG-: A=0.3, C=0.2, G=0.2, T=0.3

Sequence:
  ACGT...
""")

print("\n=== Higher-Order Markov ===")
print(""
Order 1:
  P(q_t | q_{t-1})

Order 2:
  P(q_t | q_{t-1}, q_{t-2})

Order k:
  P(q_t | q_{t-1}, ..., q_{t-k})

Trade-off:
  - Higher order: more context
  - More parameters
  - More data needed

Bioinformatics:
  - Order 1: CpG islands
  - Order 5: splice sites
  - Order k: motifs
""")

print("\n=== Stationary Distribution ===")
print(""
Stationary:
  π = πA
  π_i = Σ_j π_j a_ji

Interpretation:
  - Long-run frequencies
  - Steady state
  - Independent of initial

Example:
  Weather: π ≈ [0.6, 0.4]
""")

print("\n=== Hidden vs Observable ===")
print(""
Markov chain:
  - States observable
  - We see the states

Hidden Markov Model:
  - States hidden
  - We see emissions
  - Need to infer states

Example (weather):
  Markov: We see weather
  HMM: We see umbrella use, infer weather
""")

🎲 HMM Definition

🎯 The Model

An HMM extends a Markov chain with emissions: each hidden state emits observable symbols according to a probability distribution.

HMM

States: Hidden. Transitions: Between states. Emissions: Symbols. Initial: Start. Final: End. Parameters: λ = (A, B, π).

🎲 States🎲 Transitions🎲 Emissions🎲 Parameters

# HMM Definition
print("=== Components ===")
print(""
HMM λ = (A, B, π):

N: number of states
M: number of symbols
Q = {q_1, ..., q_N}: states
V = {v_1, ..., v_M}: symbols

A = {a_ij}: transition matrix
  a_ij = P(q_{t+1} = j | q_t = i)

B = {b_j(k)}: emission matrix
  b_j(k) = P(v_k at t | q_t = j)

π = {π_i}: initial distribution
  π_i = P(q_1 = i)
""")

print("\n=== Graphical Model ===")
print(""
q_1 → q_2 → q_3 → ... → q_T
↓     ↓     ↓           ↓
o_1   o_2   o_3   ...   o_T

Where:
  q_t: hidden state at time t
  o_t: observation at time t

Dependencies:
  q_t depends on q_{t-1}
  o_t depends on q_t
""")

print("\n=== Example: CpG Island ===")
print(""
States: CpG+ (1), CpG- (2)
Symbols: A, C, G, T

Initial:
  π = [0.5, 0.5]

Transitions:
  A = [[0.9, 0.1],
       [0.1, 0.9]]

Emissions:
  B = [[0.1, 0.4, 0.4, 0.1],   # CpG+
       [0.3, 0.2, 0.2, 0.3]]   # CpG-

Sequence: A C G T
States:   ? ? ? ?
""")

print("\n=== Example: Dishonest Casino ===")
print(""
Classic example:
  - Fair die: uniform 1-6
  - Loaded die: 6 with prob 0.5

States:
  - Fair (F)
  - Loaded (L)

Transitions:
  A = [[0.95, 0.05],
       [0.05, 0.95]]

Emissions:
  B_Fair = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6]
  B_Loaded = [0.1, 0.1, 0.1, 0.1, 0.1, 0.5]

Question:
  Given rolls, when was loaded die used?
""")

print("\n=== Parameter Estimation ===")
print(""
Known:
  - States
  - Observations

Unknown:
  - A, B, π

Estimation:
  - Supervised: labeled data
  - Unsupervised: Baum-Welch

Supervised:
  - Count transitions
  - Count emissions
  - Normalize

Unsupervised:
  - Baum-Welch (EM)
  - Iterative
""")

print("\n=== Types of HMMs ===")
print(""
1. Discrete HMM:
   - Discrete emissions
   - Most common

2. Gaussian HMM:
   - Continuous emissions
   - Gaussian distributions

3. Gaussian Mixture HMM:
   - Mixture of Gaussians
   - More flexible

4. Profile HMM:
   - From MSA
   - Match/insert/delete

5. Pair HMM:
   - Two sequences
   - Alignment

6. Input-Output HMM:
   - Input-dependent
   - Conditional
""")

print("\n=== Notation ===")
print(""
T: length of observation sequence
O = o_1 o_2 ... o_T: observations
Q = q_1 q_2 ... q_T: hidden states
λ = (A, B, π): model parameters

Forward: α_t(i) = P(o_1...o_t, q_t = i | λ)
Backward: β_t(i) = P(o_{t+1}...o_T | q_t = i, λ)
Viterbi: δ_t(i) = max P(q_1...q_{t-1}, q_t = i, o_1...o_t | λ)
""")

❓ Three Problems

🎯 The Canonical Problems

There are three canonical problems in HMMs: evaluation, decoding, and learning. Each has a classic algorithm.

Problems

Evaluation: P(O | λ). Decoding: argmax Q. Learning: argmax λ. Forward: Evaluation. Viterbi: Decoding. Baum-Welch: Learning.

❓ Evaluation❓ Decoding❓ Learning❓ Algorithms

# Three Problems
print("=== Problem 1: Evaluation ===")
print(""
Question:
  - Given O and λ
  - Compute P(O | λ)

Algorithm:
  - Forward algorithm
  - Dynamic programming
  - O(N²T)

Naive:
  - Sum over all state sequences
  - O(N^T) — infeasible

Forward:
  - α_t(i) = P(o_1...o_t, q_t = i | λ)
  - Recursion
  - O(N²T)

Use:
  - Likelihood
  - Classification
  - Model comparison
""")

print("\n=== Problem 2: Decoding ===")
print(""
Question:
  - Given O and λ
  - Find most likely Q

Algorithm:
  - Viterbi algorithm
  - Dynamic programming
  - O(N²T)

Goal:
  - argmax_Q P(Q | O, λ)
  - = argmax_Q P(Q, O | λ)

Viterbi:
  - δ_t(i) = max P(q_1...q_{t-1}, q_t = i, o_1...o_t)
  - Backpointers
  - Traceback

Use:
  - Gene prediction
  - POS tagging
  - Speech recognition
""")

print("\n=== Problem 3: Learning ===")
print(""
Question:
  - Given O
  - Find λ = (A, B, π)

Algorithm:
  - Baum-Welch (EM)
  - Iterative
  - O(N²T) per iteration

Steps:
  1. Initialize λ
  2. E-step: compute γ, ξ
  3. M-step: update λ
  4. Repeat until convergence

γ_t(i) = P(q_t = i | O, λ)
ξ_t(i,j) = P(q_t = i, q_{t+1} = j | O, λ)

Use:
  - Unsupervised learning
  - Profile HMMs
  - Clustering
""")

print("\n=== Comparison ===")
print(f"{'Problem':<20} {'Algorithm':<20} {'Complexity':<15} {'Output':<20}")
print("-" * 75)
print(f"{'Evaluation':<20} {'Forward':<20} {'O(N²T)':<15} {'P(O|λ)':<20}")
print(f"{'Decoding':<20} {'Viterbi':<20} {'O(N²T)':<15} {'Best Q':<20}")
print(f"{'Learning':<20} {'Baum-Welch':<20} {'O(N²T)':<15} {'λ*':<20}")

print("\n=== Other Problems ===")
print(""
1. Smoothing:
   - P(q_t = i | O, λ)
   - Forward-Backward

2. Prediction:
   - P(q_{t+1} = i | O, λ)
   - Forward

3. Classification:
   - Which model?
   - Likelihood ratio

4. Segmentation:
   - Change points
   - Viterbi

5. Alignment:
   - Pair HMM
   - Forward-Backward
""")

print("\n=== Complexity ===")
print(""
Forward: O(N²T)
Backward: O(N²T)
Viterbi: O(N²T)
Baum-Welch: O(N²T) per iteration

Space:
  - Forward: O(NT)
  - Viterbi: O(NT)
  - Backward: O(NT)

Optimization:
  - Scaling
  - Log-space
  - Sparse transitions
""")

➡️ Forward & Backward

🎯 Computing Likelihoods

The Forward-Backward algorithm computes the probability of observations and the posterior probability of being in each state at each time.

Forward-Backward

Forward: α. Backward: β. Likelihood: P(O|λ). Posterior: γ. Dynamic: Programming. Complexity: O(N²T).

➡️ Forward⬅️ Backward➡️ Likelihood➡️ Posterior

# Forward & Backward
print("=== Forward Algorithm ===")
print(""
Goal:
  - Compute P(O | λ)
  - Dynamic programming

Definition:
  α_t(i) = P(o_1, o_2, ..., o_t, q_t = i | λ)

Initialization:
  α_1(i) = π_i * b_i(o_1)

Recursion:
  α_{t+1}(j) = [Σ_i α_t(i) * a_ij] * b_j(o_{t+1})

Termination:
  P(O | λ) = Σ_i α_T(i)

Complexity:
  O(N²T)
""")

print("\n=== Forward Example ===")
print(""
Model:
  States: 1, 2
  Symbols: A, B
  π = [0.6, 0.4]
  A = [[0.7, 0.3],
       [0.4, 0.6]]
  B = [[0.5, 0.5],
       [0.3, 0.7]]

Observation: A B

Step 1 (t=1):
  α_1(1) = 0.6 * 0.5 = 0.30
  α_1(2) = 0.4 * 0.3 = 0.12

Step 2 (t=2):
  α_2(1) = (0.30*0.7 + 0.12*0.4) * 0.5 = 0.129
  α_2(2) = (0.30*0.3 + 0.12*0.6) * 0.7 = 0.1134

Likelihood:
  P(A B | λ) = 0.129 + 0.1134 = 0.2424
""")

print("\n=== Backward Algorithm ===")
print(""
Goal:
  - Compute β_t(i)
  - Complement to forward

Definition:
  β_t(i) = P(o_{t+1}, ..., o_T | q_t = i, λ)

Initialization:
  β_T(i) = 1

Recursion:
  β_t(i) = Σ_j a_ij * b_j(o_{t+1}) * β_{t+1}(j)

Termination:
  P(O | λ) = Σ_i π_i * b_i(o_1) * β_1(i)

Complexity:
  O(N²T)
""")

print("\n=== Posterior Probability ===")
print(""
γ_t(i):
  - P(q_t = i | O, λ)
  - Posterior

Formula:
  γ_t(i) = α_t(i) * β_t(i) / P(O | λ)

Interpretation:
  - Probability of state i at time t
  - Given all observations

Use:
  - State labeling
  - Soft classification
  - Confidence
""")

print("\n=== Xi (Transition Posterior) ===")
print(""
ξ_t(i,j):
  - P(q_t = i, q_{t+1} = j | O, λ)

Formula:
  ξ_t(i,j) = α_t(i) * a_ij * b_j(o_{t+1}) * β_{t+1}(j) / P(O | λ)

Use:
  - Baum-Welch
  - Transition estimation
""")

print("\n=== Scaling ===")
print(""
Problem:
  - α_t(i) → 0 exponentially
  - Underflow

Solution:
  - Scale α, β
  - Use log-space

Scaling:
  c_t = 1 / Σ_i α_t(i)
  α_t(i) = c_t * α_t(i)
  β_t(i) = c_t * β_t(i)

Log-space:
  - Use log-probabilities
  - Log-sum-exp trick
  - Avoid underflow
""")

print("\n=== Implementation ===")
print(""
def forward(obs, A, B, pi):
    \"\"\"Forward algorithm.\"\"\"
    T = len(obs)
    N = len(pi)
    alpha = np.zeros((T, N))
    
    # Initialization
    alpha[0] = pi * B[:, obs[0]]
    
    # Recursion
    for t in range(1, T):
        for j in range(N):
            alpha[t, j] = np.sum(
                alpha[t-1] * A[:, j]
            ) * B[j, obs[t]]
    
    # Termination
    return alpha, np.sum(alpha[-1])
""")

print("\n=== Applications ===")
print(""
Uses:
  - Likelihood calculation
  - Model comparison
  - Classification
  - Posterior probabilities
  - Baum-Welch
  - Confidence scores
""")

🏆 Viterbi Algorithm

🎯 Finding the Best Path

The Viterbi algorithm finds the single most likely sequence of hidden states given the observations.

Viterbi

Best path: argmax Q. Dynamic: Programming. Max: Not sum. Backpointers: Traceback. Complexity: O(N²T). Use: Decoding.

🏆 Best path🏆 Dynamic programming🏆 Backpointers🏆 Decoding

# Viterbi Algorithm
print("=== Definition ===")
print(""
Goal:
  - Find best state sequence
  - argmax_Q P(Q | O, λ)

Viterbi variable:
  δ_t(i) = max P(q_1...q_{t-1}, q_t = i, o_1...o_t | λ)

Backpointer:
  ψ_t(i) = argmax_j δ_{t-1}(j) * a_ji

Complexity:
  O(N²T)
""")

print("\n=== Algorithm ===")
print(""
Initialization:
  δ_1(i) = π_i * b_i(o_1)
  ψ_1(i) = 0

Recursion:
  δ_t(j) = max_i [δ_{t-1}(i) * a_ij] * b_j(o_t)
  ψ_t(j) = argmax_i [δ_{t-1}(i) * a_ij]

Termination:
  P* = max_i δ_T(i)
  q*_T = argmax_i δ_T(i)

Traceback:
  q*_t = ψ_{t+1}(q*_{t+1})
""")

print("\n=== Example ===")
print(""
Model:
  π = [0.6, 0.4]
  A = [[0.7, 0.3], [0.4, 0.6]]
  B = [[0.5, 0.5], [0.3, 0.7]]

Observation: A B

Step 1 (t=1):
  δ_1(1) = 0.6 * 0.5 = 0.30
  δ_1(2) = 0.4 * 0.3 = 0.12

Step 2 (t=2):
  δ_2(1) = max(0.30*0.7, 0.12*0.4) * 0.5 = 0.105
  δ_2(2) = max(0.30*0.3, 0.12*0.6) * 0.7 = 0.063

Termination:
  P* = 0.105
  q*_2 = 1

Traceback:
  q*_1 = 1

Best path: [1, 1]
""")

print("\n=== Implementation ===")
print(""
def viterbi(obs, A, B, pi):
    \"\"\"Viterbi algorithm.\"\"\"
    T = len(obs)
    N = len(pi)
    delta = np.zeros((T, N))
    psi = np.zeros((T, N), dtype=int)
    
    # Initialization
    delta[0] = pi * B[:, obs[0]]
    psi[0] = 0
    
    # Recursion
    for t in range(1, T):
        for j in range(N):
            probs = delta[t-1] * A[:, j]
            psi[t, j] = np.argmax(probs)
            delta[t, j] = np.max(probs) * B[j, obs[t]]
    
    # Termination
    path = np.zeros(T, dtype=int)
    path[-1] = np.argmax(delta[-1])
    prob = np.max(delta[-1])
    
    # Traceback
    for t in range(T-2, -1, -1):
        path[t] = psi[t+1, path[t+1]]
    
    return path, prob
""")

print("\n=== Viterbi vs Forward ===")
print(f"{'Aspect':<20} {'Viterbi':<20} {'Forward':<20}")
print("-" * 60)
print(f"{'Operation':<20} {'Max':<20} {'Sum':<20}")
print(f"{'Output':<20} {'Best path':<20} {'Likelihood':<20}")
print(f"{'Probability':<20} {'P(Q*, O)':<20} {'P(O)':<20}")
print(f"{'Use':<20} {'Decoding':<20} {'Evaluation':<20}")
print(f"{'Complexity':<20} {'O(N²T)':<20} {'O(N²T)':<20}")
print(f"{'Backtrack':<20} {'Yes':<20} {'No':<20}")

print("\n=== Applications ===")
print(""
Viterbi:
  - Gene prediction
  - POS tagging
  - Speech recognition
  - CpG detection
  - Splice sites
  - Protein secondary structure
  - Handwriting recognition
  - Error correction
""")

print("\n=== Variants ===")
print(""
1. Standard Viterbi:
   - Max over paths

2. Posterior Viterbi:
   - Max posterior probability

3. List Viterbi:
   - Top K paths

4. Constrained Viterbi:
   - Fixed states

5. Segment Viterbi:
   - Segmental

6. Viterbi training:
   - Approximate EM
""")

print("\n=== Log-Space ===")
print(""
Problem:
  - Underflow
  - Small probabilities

Solution:
  - Log-space
  - Addition instead of multiplication

Log-Viterbi:
  log δ_t(j) = max_i [log δ_{t-1}(i) + log a_ij] + log b_j(o_t)

Numerical:
  - Stable
  - Standard
""")

📚 Baum-Welch (EM)

🎯 Learning the Parameters

The Baum-Welch algorithm is an Expectation-Maximization (EM) algorithm that learns HMM parameters from unlabeled data.

Baum-Welch

EM: Expectation-Maximization. E-step: γ, ξ. M-step: Update A, B, π. Iterate: Until convergence. Local: Optima. Unsupervised: Learning.

📚 EM📚 E-step📚 M-step📚 Convergence

# Baum-Welch (EM)
print("=== Algorithm ===")
print(""
Goal:
  - Learn λ = (A, B, π)
  - From observations O
  - Unsupervised

Steps:
  1. Initialize λ
  2. E-step: compute γ, ξ
  3. M-step: update λ
  4. Repeat until convergence

Convergence:
  - Log-likelihood change
  - Parameter change
  - Max iterations
""")

print("\n=== E-Step ===")
print(""
Compute:
  γ_t(i) = P(q_t = i | O, λ)
  ξ_t(i,j) = P(q_t = i, q_{t+1} = j | O, λ)

Formulas:
  γ_t(i) = α_t(i) β_t(i) / P(O | λ)
  ξ_t(i,j) = α_t(i) a_ij b_j(o_{t+1}) β_{t+1}(j) / P(O | λ)

Interpretation:
  - γ: state occupancy
  - ξ: transition usage
""")

print("\n=== M-Step ===")
print(""
Update:
  π_i = γ_1(i)

  a_ij = Σ_t ξ_t(i,j) / Σ_t γ_t(i)

  b_j(k) = Σ_{t: o_t = k} γ_t(j) / Σ_t γ_t(j)

Interpretation:
  - π: initial state frequency
  - a: transition frequency
  - b: emission frequency
""")

print("\n=== Example ===")
print(""
Model:
  States: 1, 2
  Symbols: A, B

Data:
  O = A B A B

Iteration:
  1. Initialize randomly
  2. E-step: compute γ, ξ
  3. M-step: update
  4. Compute log-likelihood
  5. Repeat

Convergence:
  - Log-likelihood plateaus
  - Parameters stable
""")

print("\n=== Implementation ===")
print(""
def baum_welch(obs, N, M, max_iter=100, tol=1e-6):
    \"\"\"Baum-Welch algorithm.\"\"\"
    T = len(obs)
    
    # Initialize
    A = np.random.rand(N, N)
    A /= A.sum(axis=1, keepdims=True)
    B = np.random.rand(N, M)
    B /= B.sum(axis=1, keepdims=True)
    pi = np.random.rand(N)
    pi /= pi.sum()
    
    prev_ll = -np.inf
    
    for iteration in range(max_iter):
        # E-step
        alpha, _ = forward(obs, A, B, pi)
        beta = backward(obs, A, B)
        
        # Compute gamma and xi
        gamma = alpha * beta
        gamma /= gamma.sum(axis=1, keepdims=True)
        
        xi = np.zeros((T-1, N, N))
        for t in range(T-1):
            for i in range(N):
                for j in range(N):
                    xi[t, i, j] = (alpha[t, i] * A[i, j] *
                                   B[j, obs[t+1]] * beta[t+1, j])
            xi[t] /= xi[t].sum()
        
        # M-step
        pi = gamma[0]
        A = xi.sum(axis=0) / gamma[:-1].sum(axis=0)[:, None]
        
        for k in range(M):
            mask = obs == k
            B[:, k] = gamma[mask].sum(axis=0) / gamma.sum(axis=0)
        
        # Log-likelihood
        ll = np.log(alpha[-1].sum())
        
        if abs(ll - prev_ll) < tol:
            break
        prev_ll = ll
    
    return A, B, pi
""")

print("\n=== Initialization ===")
print(""
Random:
  - Simple
  - Multiple restarts

Uniform:
  - Equal probabilities
  - Symmetric

K-means:
  - Cluster first
  - Initialize from clusters

Manual:
  - Domain knowledge
  - Expert

Multiple restarts:
  - Avoid local optima
  - Choose best
""")

print("\n=== Convergence ===")
print(""
Criteria:
  - Log-likelihood change < tol
  - Parameter change < tol
  - Max iterations

Local optima:
  - EM guarantees non-decrease
  - Not global optimum
  - Multiple restarts

Monitoring:
  - Log-likelihood
  - Parameters
  - Iterations
""")

print("\n=== Applications ===")
print(""
Baum-Welch:
  - Profile HMMs
  - Gene prediction
  - Speech recognition
  - NLP tagging
  - Clustering
  - Anomaly detection

Variants:
  - Viterbi training
  - Segmental k-means
  - Bayesian HMM
""")

print("\n=== Best Practices ===")
print(""
1. Multiple restarts
2. Good initialization
3. Monitor log-likelihood
4. Check convergence
5. Validate
6. Regularize
7. Report parameters
8. Compare models
""")

🧬 Profile HMMs

🎯 Protein Family Models

A profile HMM is an HMM built from a multiple sequence alignment, capturing position-specific conservation and indels. It's the standard for protein family modeling.

Profile HMM

Match: Conserved. Insert: Extra residues. Delete: Missing. Position: Specific. MSA: Input. HMMER: Tool.

🧬 Match🧬 Insert🧬 Delete🧬 HMMER

# Profile HMMs
print("=== Concept ===")
print(""
Profile HMM:
  - Built from MSA
  - Position-specific
  - Match/Insert/Delete
  - Protein families

Input:
  - Multiple sequence alignment
  - Conserved columns
  - Variable columns

Output:
  - Profile HMM
  - Position-specific scores
  - Gap penalties
""")

print("\n=== Architecture ===")
print(""
States per column:
  M_k: Match
  I_k: Insert
  D_k: Delete

Transitions:
  M_k → M_{k+1}, I_k, D_{k+1}
  I_k → M_{k+1}, I_k
  D_k → M_{k+1}, D_{k+1}

Begin → M_1, I_0, D_1
End → M_L, I_L

Emissions:
  M_k: amino acid distribution
  I_k: amino acid distribution
  D_k: no emission
""")

print("\n=== Building ===")
print(""
Steps:
  1. MSA
  2. Identify match columns
  3. Count emissions
  4. Count transitions
  5. Add pseudocounts
  6. Normalize

Tools:
  - hmmbuild (HMMER)
  - HMMER
  - SAM

Example:
  hmmbuild profile.hmm alignment.sto
""")

print("\n=== HMMER ===")
print(""
HMMER:
  - Standard tool
  - Profile HMMs
  - Database search

Commands:
  hmmbuild: build HMM
  hmmsearch: search HMM vs seqs
  hmmscan: search seqs vs HMMs
  hmmscan: scan
  hmmalign: align
  hmmemit: emit
  hmmfetch: fetch

Example:
  hmmbuild profile.hmm alignment.sto
  hmmsearch profile.hmm proteome.fa
""")

print("\n=== Pfam ===")
print(""
Pfam:
  - Protein families
  - Profile HMMs
  - Comprehensive

Content:
  - 20,000+ families
  - Curated
  - HMMs

Access:
  - Web
  - FTP
  - API

URL:
  pfam.xfam.org
""")

print("\n=== InterPro ===")
print(""
InterPro:
  - Integrative
  - Many databases
  - Pfam, PRINTS, PROSITE

Content:
  - Families
  - Domains
  - Sites

Access:
  - Web
  - API

URL:
  ebi.ac.uk/interpro
""")

print("\n=== Scoring ===")
print(""
Log-odds:
  score = log2( P(seq | HMM) / P(seq | background) )

Bit score:
  - Log-odds
  - Bits

E-value:
  - Expected hits
  - Database size
  - < 1e-5: significant

HMMER output:
  - Target
  - Query
  - E-value
  - Score
  - Bias
""")

print("\n=== Applications ===")
print(""
Profile HMMs:
  - Protein family classification
  - Remote homology detection
  - Domain annotation
  - Genome annotation
  - Metagenomics
  - Phylogenetics

Example:
  - Kinase domain
  - Zinc finger
  - Immunoglobulin
  - ABC transporter
""")

print("\n=== Building Custom HMM ===")
print(""
# 1. Prepare MSA
muscle -in seqs.fa -out aln.fa

# 2. Convert to Stockholm
# (or use HMMER's format)

# 3. Build HMM
hmmbuild profile.hmm aln.sto

# 4. Search
hmmsearch profile.hmm proteome.fa

# 5. Align hits
hmmalign profile.hmm hits.fa > hits.sto
""")

print("\n=== Best Practices ===")
print(""
1. Good MSA
2. Enough sequences
3. Curate alignment
4. Use HMMER
5. Check E-values
6. Validate hits
7. Compare to Pfam
8. Report
""")

🔬 HMMER

🎯 The Standard Tool

HMMER is the standard tool for profile HMM-based sequence analysis, used for protein family classification and remote homology detection.

HMMER

hmmbuild: Build. hmmsearch: Search. hmmscan: Scan. hmmalign: Align. Pfam: Database. E-value: Significance.

🔬 hmmbuild🔬 hmmsearch🔬 hmmscan🔬 hmmalign

# HMMER
print("=== Installation ===")
print(""
# Conda
conda install -c bioconda hmmer

# Ubuntu
sudo apt-get install hmmer

# macOS
brew install hmmer

# From source
wget http://eddylab.org/software/hmmer/hmmer.tar.gz
tar xzf hmmer.tar.gz
cd hmmer-3.4
./configure
make
make install
""")

print("\n=== hmmbuild ===")
print(""
Build profile HMM from MSA.

Command:
  hmmbuild profile.hmm alignment.sto

Options:
  --amino: Protein (default)
  --dna: DNA
  --rna: RNA
  --informat: Input format
  --symfrac: Symfrac
  --fragthresh: Fragment

Output:
  profile.hmm
""")

print("\n=== hmmsearch ===")
print(""
Search HMM vs sequence database.

Command:
  hmmsearch profile.hmm proteome.fa > results.txt

Options:
  -o: Output file
  -A: Alignment output
  --tblout: Table output
  --domtblout: Domain table
  -E: E-value threshold
  --domE: Domain E-value

Output:
  results.txt
  results.tblout
  results.domtblout
""")

print("\n=== hmmscan ===")
print(""
Scan sequence vs HMM database.

Command:
  hmmscan Pfam-A.hmm protein.fa > results.txt

Options:
  --tblout: Table
  --domtblout: Domain table
  -E: E-value
  --cut_ga: Gathering threshold

Output:
  results.txt
  results.tblout
""")

print("\n=== hmmalign ===")
print(""
Align sequences to profile HMM.

Command:
  hmmalign profile.hmm sequences.fa > aligned.sto

Options:
  -o: Output
  --trim: Trim
  --outformat: Format

Output:
  aligned.sto
""")

print("\n=== Other Commands ===")
print(""
hmmfetch:
  - Fetch HMM from database

hmmemit:
  - Emit sequences from HMM

hmmpress:
  - Index HMM database

hmmlogo:
  - Logo from HMM

hmmstat:
  - Statistics
""")

print("\n=== Pfam Database ===")
print(""
Pfam:
  - Protein families
  - 20,000+ HMMs
  - Curated

Download:
  wget http://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
  gunzip Pfam-A.hmm.gz

Index:
  hmmpress Pfam-A.hmm

Scan:
  hmmscan --cut_ga Pfam-A.hmm protein.fa > domains.txt
""")

print("\n=== Example Workflow ===")
print(""
# 1. Get sequences
# 2. Align (MAFFT)
mafft sequences.fa > aligned.fa

# 3. Convert to Stockholm
# (or use aligned.fa)

# 4. Build HMM
hmmbuild profile.hmm aligned.sto

# 5. Search proteome
hmmsearch --tblout hits.tbl profile.hmm proteome.fa

# 6. Filter by E-value
awk '$5 < 1e-5' hits.tbl > significant.tbl

# 7. Align hits
hmmalign profile.hmm hits.fa > hits.sto
""")

print("\n=== Python with pyhmmer ===")
print(""
import pyhmmer

# Load HMM
with pyhmmer.plan7.HMMFile("profile.hmm") as hmm_file:
    hmm = hmm_file.read()

# Load sequences
with pyhmmer.easel.SequenceFile("proteome.fa",
                                 digital=True) as seq_file:
    sequences = seq_file.read_block()

# Search
hits = pyhmmer.hmmsearch([hmm], sequences)

for hit in hits:
    for domain in hit.domains:
        print(f"Hit: {hit.name}")
        print(f"  E-value: {hit.evalue}")
        print(f"  Score: {hit.score}")
""")

print("\n=== Best Practices ===")
print(""
1. Use --cut_ga for Pfam
2. Filter by E-value
3. Check domain boundaries
4. Validate hits
5. Use multiple databases
6. Compare to BLAST
7. Report parameters
8. Cite HMMER
""")

🧬 Gene Prediction

🎯 Finding Genes with HMMs

Gene prediction uses HMMs to identify genes (exons, introns, splice sites) in genomic DNA.

Gene Prediction

States: Exon, intron. Emissions: Codons. Transitions: Splice. GENSCAN: Classic. AUGUSTUS: Modern. Glimmer: Prokaryotes.

🧬 Exon🧬 Intron🧬 Splice sites🧬 GENSCAN

# Gene Prediction
print("=== HMM Architecture ===")
print(""
Eukaryotic gene:
  - Exons
  - Introns
  - Splice sites
  - UTRs
  - Promoters

HMM states:
  - E: Exon
  - I: Intron
  - N: Intergenic
  - 5' UTR
  - 3' UTR
  - Promoter

Transitions:
  - E → I (donor site)
  - I → E (acceptor site)
  - N → E (start)
  - E → N (stop)
""")

print("\n=== GENSCAN ===")
print(""
GENSCAN (Burge & Karlin, 1997):
  - Classic
  - Eukaryotic
  - HMM-based

Features:
  - Multiple genes
  - Partial genes
  - Suboptimal exons
  - CpG islands
  - Promoters

States:
  - Exon
  - Intron
  - Intergenic
  - UTR

Web:
  genes.mit.edu/GENSCAN.html
""")

print("\n=== AUGUSTUS ===")
print(""
AUGUSTUS:
  - Modern
  - Eukaryotic
  - HMM-based

Features:
  - Multiple species
  - Training
  - RNA-Seq evidence
  - Homology

Command:
  augustus --species=human genome.fa > predictions.gff

Options:
  --species: Species model
  --gff3: Output format
  --protein: Protein output
  --codingseq: CDS output

Web:
  bioinf.uni-greifswald.de/augustus
""")

print("\n=== Glimmer ===")
print(""
Glimmer:
  - Prokaryotic
  - Interpolated Markov Models
  - Classic

Features:
  - Bacterial genomes
  - Archaea
  - Fast

Command:
  glimmer3 genome.fa output

Web:
  ccb.jhu.edu/software/glimmer
""")

print("\n=== GeneMark ===")
print(""
GeneMark:
  - Prokaryotic + eukaryotic
  - HMM-based
  - Self-training

Features:
  - Bacteria
  - Archaea
  - Fungi
  - Plants

Web:
  opal.biology.gatech.edu/GeneMark
""")

print("\n=== SNAP ===")
print(""
SNAP:
  - Eukaryotic
  - Semi-HMM
  - Training

Features:
  - Human
  - Drosophila
  - Arabidopsis

Web:
  home.snafu.de/korf/SNAP
""")

print("\n=== Evaluation ===")
print(""
Metrics:
  - Sensitivity (Sn)
  - Specificity (Sp)
  - Accuracy
  - CC (correlation)

At gene level:
  - Exact match
  - Partial match
  - Missing

At exon level:
  - Exon Sn
  - Exon Sp

Benchmarks:
  - GASP
  - ENCODE
  - EGASP
""")

print("\n=== Best Practices ===")
print(""
1. Use species-specific model
2. Train if possible
3. Use RNA-Seq evidence
4. Combine multiple tools
5. Validate manually
6. Check conserved genes
7. Report parameters
8. Compare to reference
""")

🐍 Python Integration

🎯 HMMs in Python

You can implement HMMs from scratch or use libraries like hmmlearn, pomegranate, and pyhmmer.

Python

hmmlearn: General. pomegranate: Flexible. pyhmmer: HMMER. NumPy: Scratch. Scratch: Learning. Pipeline: Automate.

🐍 hmmlearn🐍 pomegranate🐍 pyhmmer🐍 NumPy

# Python Integration
print("=== From Scratch: Viterbi ===")
print(""
import numpy as np

def viterbi(obs, A, B, pi):
    \"\"\"Viterbi algorithm.\"\"\"
    T = len(obs)
    N = len(pi)
    delta = np.zeros((T, N))
    psi = np.zeros((T, N), dtype=int)
    
    # Initialization
    delta[0] = pi * B[:, obs[0]]
    
    # Recursion
    for t in range(1, T):
        for j in range(N):
            probs = delta[t-1] * A[:, j]
            psi[t, j] = np.argmax(probs)
            delta[t, j] = np.max(probs) * B[j, obs[t]]
    
    # Termination
    path = np.zeros(T, dtype=int)
    path[-1] = np.argmax(delta[-1])
    prob = np.max(delta[-1])
    
    # Traceback
    for t in range(T-2, -1, -1):
        path[t] = psi[t+1, path[t+1]]
    
    return path, prob

# Example
obs = [0, 1, 0, 1]  # A B A B
A = np.array([[0.7, 0.3], [0.4, 0.6]])
B = np.array([[0.5, 0.5], [0.3, 0.7]])
pi = np.array([0.6, 0.4])

path, prob = viterbi(obs, A, B, pi)
print(f"Path: {path}")
print(f"Probability: {prob:.4f}")
""")

print("\n=== From Scratch: Forward ===")
print(""
def forward(obs, A, B, pi):
    \"\"\"Forward algorithm.\"\"\"
    T = len(obs)
    N = len(pi)
    alpha = np.zeros((T, N))
    
    # Initialization
    alpha[0] = pi * B[:, obs[0]]
    
    # Recursion
    for t in range(1, T):
        for j in range(N):
            alpha[t, j] = np.sum(
                alpha[t-1] * A[:, j]
            ) * B[j, obs[t]]
    
    # Termination
    return alpha, np.sum(alpha[-1])

# Use
alpha, likelihood = forward(obs, A, B, pi)
print(f"Likelihood: {likelihood:.4f}")
""")

print("\n=== From Scratch: Baum-Welch ===")
print(""
def baum_welch(obs, N, M, max_iter=100, tol=1e-6):
    \"\"\"Baum-Welch algorithm.\"\"\"
    T = len(obs)
    
    # Initialize
    A = np.random.rand(N, N)
    A /= A.sum(axis=1, keepdims=True)
    B = np.random.rand(N, M)
    B /= B.sum(axis=1, keepdims=True)
    pi = np.random.rand(N)
    pi /= pi.sum()
    
    prev_ll = -np.inf
    
    for _ in range(max_iter):
        # E-step
        alpha, _ = forward(obs, A, B, pi)
        beta = backward(obs, A, B)
        
        gamma = alpha * beta
        gamma /= gamma.sum(axis=1, keepdims=True)
        
        xi = np.zeros((T-1, N, N))
        for t in range(T-1):
            for i in range(N):
                for j in range(N):
                    xi[t, i, j] = (alpha[t, i] * A[i, j] *
                                   B[j, obs[t+1]] * beta[t+1, j])
            xi[t] /= xi[t].sum()
        
        # M-step
        pi = gamma[0]
        A = xi.sum(axis=0) / gamma[:-1].sum(axis=0)[:, None]
        
        for k in range(M):
            B[:, k] = gamma[obs == k].sum(axis=0) / gamma.sum(axis=0)
        
        # Log-likelihood
        ll = np.log(alpha[-1].sum())
        if abs(ll - prev_ll) < tol:
            break
        prev_ll = ll
    
    return A, B, pi

# Use
A, B, pi = baum_welch(obs, N=2, M=2)
""")

print("\n=== hmmlearn ===")
print(""
from hmmlearn import hmm
import numpy as np

# Discrete HMM
model = hmm.MultinomialHMM(n_components=2, n_iter=100)
model.fit(obs.reshape(-1, 1))

# Predict
hidden_states = model.predict(obs.reshape(-1, 1))

# Parameters
print(f"Transitions: {model.transmat_}")
print(f"Emissions: {model.emissionprob_}")
print(f"Start: {model.startprob_}")

# Score
log_prob = model.score(obs.reshape(-1, 1))
print(f"Log-likelihood: {log_prob}")

# Gaussian HMM
model = hmm.GaussianHMM(n_components=3)
model.fit(X)
""")

print("\n=== pomegranate ===")
print(""
from pomegranate import HiddenMarkovModel, DiscreteDistribution

# Define distributions
d1 = DiscreteDistribution({'A': 0.5, 'B': 0.5})
d2 = DiscreteDistribution({'A': 0.3, 'B': 0.7})

# Create model
model = HiddenMarkovModel()
model.add_states(d1, d2)

# Add transitions
model.add_transition(model.start, d1, 0.6)
model.add_transition(model.start, d2, 0.4)
model.add_transition(d1, d1, 0.7)
model.add_transition(d1, d2, 0.3)
model.add_transition(d2, d1, 0.4)
model.add_transition(d2, d2, 0.6)
model.bake()

# Predict
log_prob, path = model.viterbi(['A', 'B', 'A'])
""")

print("\n=== pyhmmer ===")
print(""
import pyhmmer

# Load HMM
with pyhmmer.plan7.HMMFile("profile.hmm") as f:
    hmm = f.read()

# Load sequences
with pyhmmer.easel.SequenceFile("proteome.fa",
                                 digital=True) as f:
    seqs = f.read_block()

# Search
for hits in pyhmmer.hmmsearch([hmm], seqs):
    for hit in hits:
        print(f"Hit: {hit.name}, E-value: {hit.evalue}")
""")

print("\n=== Complete Pipeline ===")
print(""
import numpy as np
from hmmlearn import hmm

def hmm_pipeline(sequences, n_states=3):
    \"\"\"Complete HMM pipeline.\"\"\"
    # Prepare data
    X = np.array(sequences).reshape(-1, 1)
    
    # Train
    model = hmm.MultinomialHMM(n_components=n_states,
                                n_iter=100)
    model.fit(X)
    
    # Predict
    states = model.predict(X)
    
    # Score
    score = model.score(X)
    
    print(f"States: {n_states}")
    print(f"Log-likelihood: {score:.2f}")
    print(f"Transitions:\\n{model.transmat_}")
    
    return model, states

# Use
# model, states = hmm_pipeline(sequences)
""")

⚖️ Tool Comparison

🎯 Which Tool?

Choosing the right HMM tool depends on your application: general modeling, profile HMMs, or gene prediction.

Comparison

hmmlearn: General. pomegranate: Flexible. pyhmmer: HMMER. HMMER: Profile. AUGUSTUS: Gene. GENSCAN: Classic.

⚖️ hmmlearn⚖️ pomegranate⚖️ HMMER⚖️ AUGUSTUS

Tool	Type	Language	Best For
hmmlearn	Library	Python	General HMMs
pomegranate	Library	Python	Flexible models
pyhmmer	Library	Python	HMMER bindings
HMMER	CLI	C	Profile HMMs
HMMER3	CLI	C	Protein search
Infernal	CLI	C	RNA families
AUGUSTUS	CLI	C++	Eukaryotic genes
GENSCAN	Web	C	Eukaryotic genes
Glimmer	CLI	C++	Prokaryotic genes
GeneMark	CLI	C	Prokaryotic + eukaryotic

# Decision Guide
print("=== Choose Tool by Application ===")
print(""
General HMM:
  → hmmlearn (Python)
  → pomegranate (Python)

Profile HMMs:
  → HMMER (CLI)
  → pyhmmer (Python)
  → Pfam database

RNA families:
  → Infernal
  → Rfam

Gene prediction:
  → AUGUSTUS (eukaryotic)
  → Glimmer (prokaryotic)
  → GeneMark

Speech:
  → HTK
  → Kaldi

NLP:
  → NLTK
  → hmmlearn
""")

print("\n=== Speed Comparison ===")
print(f"{'Tool':<20} {'100 seqs':<15} {'1000 seqs':<15}")
print("-" * 50)
print(f"{'hmmlearn':<20} {'1s':<15} {'10s':<15}")
print(f"{'pomegranate':<20} {'1s':<15} {'10s':<15}")
print(f"{'HMMER':<20} {'1s':<15} {'10s':<15}")
print(f"{'AUGUSTUS':<20} {'10s':<15} {'2 min':<15}")
print(f"{'Glimmer':<20} {'5s':<15} {'1 min':<15}")

print("\n=== When to Use Which ===")
print(""
USE hmmlearn:
  ✓ General HMMs
  ✓ Python
  ✓ scikit-learn style

USE pomegranate:
  ✓ Flexible
  ✓ Complex models
  ✓ Bayesian

USE HMMER:
  ✓ Protein families
  ✓ Remote homology
  ✓ Pfam

USE Infernal:
  ✓ RNA families
  ✓ Rfam
  ✓ Structure

USE AUGUSTUS:
  ✓ Eukaryotic genes
  ✓ RNA-Seq evidence
  ✓ Multiple species

USE Glimmer:
  ✓ Prokaryotic genes
  ✓ Fast
  ✓ Bacterial genomes
""")

print("\n=== Complexity ===")
print(f"{'Algorithm':<25} {'Time':<15} {'Space':<15}")
print("-" * 55)
print(f"{'Forward':<25} {'O(N²T)':<15} {'O(NT)':<15}")
print(f"{'Backward':<25} {'O(N²T)':<15} {'O(NT)':<15}")
print(f"{'Viterbi':<25} {'O(N²T)':<15} {'O(NT)':<15}")
print(f"{'Baum-Welch':<25} {'O(N²T)':<15} {'O(NT)':<15}")
print(f"{'Forward-Backward':<25} {'O(N²T)':<15} {'O(NT)':<15}")
""")

💼 Use Cases

🎯 Where HMMs Are Used

HMMs are used in virtually every area of sequence analysis: gene prediction, protein families, speech, and more.

Use Cases

Gene: Prediction. Protein: Families. Alignment: Profile. Speech: Recognition. NLP: Tagging. Finance: Regimes.

💼 Gene prediction💼 Protein families💼 Speech💼 NLP

# Use Cases
print("=== 1. Gene Prediction ===")
print(""
Use case:
  - Find genes
  - Exons, introns
  - Splice sites

Workflow:
  1. Genome
  2. HMM (AUGUSTUS)
  3. Predictions
  4. Validate

Examples:
  - Human
  - Mouse
  - Arabidopsis
  - Bacteria

Tools:
  - AUGUSTUS
  - GENSCAN
  - Glimmer
  - GeneMark
""")

print("\n=== 2. Protein Family Classification ===")
print(""
Use case:
  - Classify proteins
  - Domains
  - Families

Workflow:
  1. Protein sequence
  2. HMMER vs Pfam
  3. Domains
  4. Function

Examples:
  - Kinases
  - GPCRs
  - Zinc fingers
  - Immunoglobulins

Tools:
  - HMMER
  - Pfam
  - InterPro
""")

print("\n=== 3. Remote Homology Detection ===")
print(""
Use case:
  - Find distant homologs
  - No BLAST hits
  - Profile HMM

Workflow:
  1. Build HMM
  2. Search database
  3. E-values
  4. Validate

Examples:
  - Protein families
  - Metagenomics
  - Novel proteins

Tools:
  - HMMER
  - HHpred
  - PSI-BLAST
""")

print("\n=== 4. CpG Island Detection ===")
print(""
Use case:
  - Find CpG islands
  - Promoters
  - Methylation

Workflow:
  1. DNA sequence
  2. HMM
  3. CpG+ regions
  4. Validate

Tools:
  - CpGIE
  - CpGProD
  - Custom HMM
""")

print("\n=== 5. Splice Site Prediction ===")
print(""
Use case:
  - Find splice sites
  - Donor, acceptor
  - Gene structure

Workflow:
  1. DNA sequence
  2. HMM
  3. Splice sites
  4. Validate

Tools:
  - GENSCAN
  - AUGUSTUS
  - SplicePort
""")

print("\n=== 6. Speech Recognition ===")
print(""
Use case:
  - Acoustic → words
  - Phonemes
  - HMM states

Workflow:
  1. Audio
  2. Features
  3. HMM
  4. Decode

Tools:
  - HTK
  - Kaldi
  - CMU Sphinx
""")

print("\n=== 7. NLP Tagging ===")
print(""
Use case:
  - POS tagging
  - NER
  - Chunking

Workflow:
  1. Text
  2. HMM
  3. Tags
  4. Evaluate

Tools:
  - NLTK
  - hmmlearn
  - pomegranate
""")

print("\n=== 8. Financial Modeling ===")
print(""
Use case:
  - Regime switching
  - Volatility
  - Markets

Workflow:
  1. Time series
  2. HMM
  3. Regimes
  4. Predict

Tools:
  - hmmlearn
  - pomegranate
  - Statsmodels
""")

print("\n=== Comparison ===")
print(f"{'Use Case':<25} {'Tool':<20} {'Output':<20}")
print("-" * 65)
print(f"{'Gene prediction':<25} {'AUGUSTUS':<20} {'Genes':<20}")
print(f"{'Protein families':<25} {'HMMER':<20} {'Domains':<20}")
print(f"{'Remote homology':<25} {'HMMER':<20} {'Hits':<20}")
print(f"{'CpG islands':<25} {'CpGIE':<20} {'Islands':<20}")
print(f"{'Splice sites':<25} {'GENSCAN':<20} {'Sites':<20}")
print(f"{'Speech':<25} {'Kaldi':<20} {'Words':<20}")
print(f"{'NLP':<25} {'NLTK':<20} {'Tags':<20}")
print(f"{'Finance':<25} {'hmmlearn':<20} {'Regimes':<20}")

❓ FAQ

🎯 Frequently Asked Questions

Answers to common questions about Hidden Markov Models.

What is a Hidden Markov Model?

What is the difference between Viterbi, Forward-Backward, and Baum-Welch?

Viterbi finds the single most likely sequence of hidden states given observations. Forward-Backward computes the probability of being in each state at each position (posterior probabilities). Baum-Welch is an Expectation-Maximization algorithm that learns the model parameters (transition and emission probabilities) from unlabeled data.

What is a profile HMM and what is HMMER used for?

A profile HMM is an HMM built from a multiple sequence alignment of a protein family. It captures position-specific conservation and insertion/deletion patterns. HMMER is the standard tool for building profile HMMs and searching databases for remote homologs (e.g., Pfam, InterPro).

What are the three main problems of HMMs?

The three canonical HMM problems are: (1) Evaluation: compute the probability of an observation sequence given the model (Forward algorithm); (2) Decoding: find the most likely hidden state sequence (Viterbi algorithm); (3) Learning: estimate model parameters from data (Baum-Welch / EM algorithm).

Which Python library should I use for HMMs?

Use hmmlearn for general-purpose discrete and Gaussian HMMs (scikit-learn-style API). Use pomegranate for flexible probabilistic models. For profile HMMs and biological sequence search, use HMMER (pyhmmer for Python bindings). For teaching, implement Viterbi/Forward-Backward from scratch.

What is the difference between a Markov chain and an HMM?

A Markov chain has observable states; we see the state at each time step. An HMM has hidden states; we only see emissions generated by those states. The states in an HMM must be inferred from the observations, which is what Viterbi and Forward-Backward do.

How do I choose the number of states in an HMM?

There's no universal rule. Use domain knowledge (e.g., 2 states for CpG+/CpG-, 3 for exon/intron/intergenic), model selection (AIC/BIC, cross-validation), or Bayesian nonparametrics (infinite HMM). Start simple and increase complexity only if needed. Always validate with held-out data.

What is the difference between HMMER and BLAST?

BLAST uses pairwise sequence alignment and is fast for close homologs. HMMER uses profile HMMs built from MSAs, capturing position-specific conservation and indels. HMMER is more sensitive for remote homology detection — finding distant relatives that BLAST misses. They are complementary.

How do I interpret HMMER E-values?

E-value is the expected number of hits with a given score or better occurring by chance in a database of a given size. Lower is better: E-value < 1e-5 = significant, 1e-5 to 1e-3 = marginal, > 1e-3 = not significant. For Pfam, use the curated gathering threshold (--cut_ga).

What is the difference between a profile HMM and a PWM?

A PWM (Position Weight Matrix) is a fixed-length motif representation with no gap states. A profile HMM extends this with explicit insert and delete states, allowing variable-length matches with gaps. Profile HMMs are more powerful for protein families where insertions and deletions occur.

# Quick Reference
print("=== Quick Answers ===")
print(""
Q: What is an HMM?
A: Probabilistic model with hidden states.

Q: Three problems?
A: Evaluation, decoding, learning.

Q: Viterbi?
A: Best state sequence.

Q: Forward?
A: Likelihood P(O|λ).

Q: Baum-Welch?
A: EM for learning parameters.

Q: Profile HMM?
A: Protein family model.

Q: HMMER?
A: Profile HMM tool.

Q: Python library?
A: hmmlearn, pomegranate, pyhmmer.

Q: Gene prediction?
A: AUGUSTUS, GENSCAN, Glimmer.

Q: Free tools?
A: Yes — HMMER, hmmlearn, AUGUSTUS.
""")

🚀 Future

🎯 Where HMMs Are Going

The future is deep learning hybrids, single-cell, multi-omics, real-time, and interpretable AI.

Future Directions

Deep learning: Hybrids. Single-cell: scHMM. Multi-omics: Integration. Real-time: Streaming. Interpretable: AI. Foundation: Still relevant.

🚀 Deep learning🚀 Single-cell🚀 Multi-omics🚀 Real-time

# Future of HMMs
print("=== 1. Deep Learning Hybrids ===")
print(""
Trend: HMM + deep learning

Ideas:
  - Neural HMM
  - HMM + CNN
  - HMM + transformer
  - Differentiable HMM

Examples:
  - DeepHMM
  - Neural HMM
  - RNN-HMM

Future:
  - End-to-end
  - Learnable
  - Multi-modal
""")

print("\n=== 2. Single-Cell ===")
print(""
Trend: Single-cell

Current:
  - Bulk
  - Average

Single-cell:
  - scHMM
  - Per-cell states
  - Cell-type specific
  - Trajectories

Tools:
  - Monocle
  - Seurat
  - HMM-based

Future:
  - Cell-type specific
  - Rare populations
  - Developmental
""")

print("\n=== 3. Multi-Omics ===")
print(""
Trend: Integrate everything

Current:
  - Single modality
  - Separate analysis

Multi-omics:
  - HMM + RNA
  - HMM + ATAC
  - HMM + Hi-C
  - Single-cell multi-omics

Tools:
  - MOFA
  - Seurat
  - HMM-based

Future:
  - Unified models
  - Regulatory networks
  - Systems biology
""")

print("\n=== 4. Real-Time ===")
print(""
Trend: Real-time

Current:
  - Batch
  - Slow

Real-time:
  - Streaming
  - Online
  - Continuous
  - Adaptive

Examples:
  - Nanopore
  - Real-time surveillance
  - Wearables

Tools:
  - Online HMM
  - Streaming HMM
  - Adaptive HMM
""")

print("\n=== 5. Interpretable AI ===")
print(""
Trend: Interpretable

Current:
  - Black box
  - Deep learning

Future:
  - Interpretable
  - HMM + attention
  - Explainable
  - Trustworthy

Examples:
  - Attention-HMM
  - Explainable HMM
  - Causal HMM

Future:
  - Clinical
  - Regulatory
  - Trust
""")

print("\n=== 6. Foundation Models ===")
print(""
Trend: Foundation models

Current:
  - Task-specific
  - Small

Future:
  - Pretrained
  - Transfer
  - Multi-task
  - HMM + LLM

Examples:
  - Protein language models
  - DNA language models
  - HMM hybrids

Future:
  - Unified
  - Scalable
  - General
""")

# TIMELINE
print("\n=== Timeline ===")
print(""
Past:
  1960s: Baum-Welch
  1966: Viterbi
  1970s: Speech
  1989: Profile HMM
  1994: HMMER
  1998: GENSCAN
  2000s: Pfam
  2010s: Deep learning
  2020s: Hybrids

Present:
  2025: Single-cell
  2026: Multi-omics

Future:
  2030: Deep learning hybrids?
  2035: Real-time?
  2040: Foundation models?
""")

print("\n=== Final Thoughts ===")
print(""
HMMs are:
  - Foundational
  - Probabilistic
  - Interpretable
  - Still relevant

The future is:
  - Deep learning hybrids
  - Single-cell
  - Multi-omics
  - Real-time
  - Interpretable AI
  - Foundation models

"HMMs gave us a probabilistic
 language for sequences.
 That language is still spoken."

Keep exploring. Keep building.
""")