9: Hidden Markov Models (HMMs) in Biological Sequence Analysis
Chapter 9: Hidden Markov Models (HMMs) in Biological Sequence Analysis
Biological sequences are generated by underlying, unobserved biological states: an exon state, an intron state, a splice junction, a transmembrane domain, or a CpG island. Hidden Markov Models (HMMs) provide an elegant probabilistic framework for modeling sequential dependencies and decoding hidden biological features directly from primary sequence data.
An HMM is defined by five mathematical components: \(\lambda = (S, V, A, B, \pi)\):
- Hidden States (\(S\)): The unobserved biological labels (e.g., \(S = \{ ext{Exon}, ext{Intron}, ext{5'SS}, ext{3'SS}\}\)).
- Observation Alphabet (\(V\)): The emitted characters (e.g., \(V = \{A, C, G, T\}\)).
- Transition Probabilities (\(A\)): Matrix \(A_{ij} = P(q_{t+1} = S_j \mid q_t = S_i)\), representing the probability of transitioning from hidden state \(i\) to state \(j\).
- Emission Probabilities (\(B\)): Matrix \(B_j(k) = P(O_t = v_k \mid q_t = S_j)\), the probability of emitting character \(k\) while in state \(j\).
- Initial State Distribution (\(\pi\)): Vector \(\pi_i = P(q_1 = S_i)\).
1. The Viterbi Decoding Algorithm
Given an observed DNA sequence \(O = (O_1, O_2, \dots, O_T)\), the Viterbi algorithm uses dynamic programming to find the single most probable hidden state path \(Q^* = (q_1^*, q_2^*, \dots, q_T^*)\):
vt(j) = maxi [ vt-1(i) × Aij ] × Bj(Ot)
Viterbi guarantees optimal global state path decoding in \(O(K^2 T)\) time, where \(K\) is the number of states and \(T\) is sequence length.
2. Protein Domain Profiling with HMMER3
Profile HMMs capture position-specific insertion and deletion probabilities across protein families, dramatically outperforming standard BLAST for remote homology detection:
# 1. Build a profile HMM from a multiple sequence alignment
hmmbuild kinase_domain.hmm aligned_kinases.fasta
# 2. Search a translated proteome for members containing the kinase domain
hmmsearch \
--domE 1e-5 \
--cpu 8 \
kinase_domain.hmm \
human_proteome.fasta > kinase_hits.txt
# 3. Scan a single query protein against the entire Pfam domain database
hmmscan \
--domE 1e-5 \
/databases/Pfam-A.hmm \
query_protein.fasta > domain_architecture.txtCommon Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.