Example Program
HMM Silent States
Hidden Markov Model with silent states code example
A tutorial about HMMs with silent states
1#include <iostream>
2#include <fstream>
3#include <seqan/graph_algorithms.h>
4#include <seqan/basic/basic_logvalue.h>
5
6using namespace seqan;
7
8int main() {
HMM creation
9    typedef LogProb<> TProbability;
10    typedef Dna TAlphabet;
11    typedef Size<TAlphabet>::Type TSize;
12    typedef Graph<Hmm<TAlphabet, TProbability> > THmm;
13    typedef VertexDescriptor<THmm>::Type TVertexDescriptor;
14    typedef EdgeDescriptor<THmm>::Type TEdgeDescriptor;
15    
16    Dna dnaA = Dna('A');
17    Dna dnaC = Dna('C');
18    Dna dnaG = Dna('G');
19    Dna dnaT = Dna('T');
20
21    THmm hmm;
22
Begin state
23    TVertexDescriptor begState = addVertex(hmm);
24    assignBeginState(hmm, begState);
25
Emission states
26    TVertexDescriptor emitState1 = addVertex(hmm);
27    emissionProbability(hmm, emitState1, dnaA) = 0.0;
28    emissionProbability(hmm, emitState1, dnaC) = 0.8;
29    emissionProbability(hmm, emitState1, dnaG) = 0.1;
30    emissionProbability(hmm, emitState1, dnaT) = 0.1;
31
32    TVertexDescriptor emitState2 = addVertex(hmm);
33    emissionProbability(hmm, emitState2, dnaA) = 0.0;
34    emissionProbability(hmm, emitState2, dnaC) = 0.2;
35    emissionProbability(hmm, emitState2, dnaG) = 0.2;
36    emissionProbability(hmm, emitState2, dnaT) = 0.6;
37
38    TVertexDescriptor emitState3 = addVertex(hmm);
39    emissionProbability(hmm, emitState3, dnaA) = 0.7;
40    emissionProbability(hmm, emitState3, dnaC) = 0.1;
41    emissionProbability(hmm, emitState3, dnaG) = 0.1;
42    emissionProbability(hmm, emitState3, dnaT) = 0.1;
43
44    TVertexDescriptor emitState4 = addVertex(hmm);
45    emissionProbability(hmm, emitState4, dnaA) = 0.25;
46    emissionProbability(hmm, emitState4, dnaC) = 0.25;
47    emissionProbability(hmm, emitState4, dnaG) = 0.25;
48    emissionProbability(hmm, emitState4, dnaT) = 0.25;
49
Silent states (deletion states)
50    TVertexDescriptor delState1 = addVertex(hmm, true);
51    TVertexDescriptor delState2 = addVertex(hmm, true);
52    TVertexDescriptor delState3 = addVertex(hmm, true);
53    TVertexDescriptor delState4 = addVertex(hmm, true);
54
End state
55    TVertexDescriptor eState = addVertex(hmm);
56    assignEndState(hmm, eState);
57
Transitions
58    addEdge(hmm, begState, emitState1, 0.5);
59    addEdge(hmm, begState, delState1, 0.5);
60    addEdge(hmm, emitState1, emitState2, 0.5);
61    addEdge(hmm, emitState1, delState2, 0.5);
62    addEdge(hmm, delState1, emitState2, 0.5);
63    addEdge(hmm, delState1, delState2, 0.5);
64    addEdge(hmm, emitState2, emitState3, 0.5);
65    addEdge(hmm, emitState2, delState3, 0.5);
66    addEdge(hmm, delState2, emitState3, 0.5);
67    addEdge(hmm, delState2, delState3, 0.5);
68    addEdge(hmm, emitState3, emitState4, 0.5);
69    addEdge(hmm, emitState3, delState4, 0.5);
70    addEdge(hmm, delState3, emitState4, 0.5);
71    addEdge(hmm, delState3, delState4, 0.5);
72    addEdge(hmm, emitState4, eState, 1.0);
73    addEdge(hmm, delState4, eState, 1.0);
74
Print the whole model
75    ::std::cout << hmm << ::std::endl;
76
Viterbi algorithm
77    String<Dna> sequence = "CA";
78    String<TVertexDescriptor> path;
79    TProbability p = viterbiAlgorithm(hmm, sequence, path);
80    ::std::cout << "Viterbi algorithm" << ::std::endl;
81    ::std::cout << "Probability of the best path: " << p << ::std::endl;
82    ::std::cout << "Sequence: " << ::std::endl;
83    for(TSize i = 0; i<length(sequence); ++i) ::std::cout << sequence[i] << ',';
84    ::std::cout << ::std::endl;
85    ::std::cout << "State path: " << ::std::endl;
86    for(TSize i = 0; i<length(path); ++i) {
87        ::std::cout << path[i];
88        if (isSilent(hmm, path[i])) ::std::cout << " (Silent)";
89        if (i < length(path) - 1) ::std::cout << ',';
90    }
91    ::std::cout << ::std::endl;
92
Forward algorithm
93    ::std::cout << "Forward algorithm" << ::std::endl;
94    p = forwardAlgorithm(hmm, sequence);
95    ::std::cout << "Probability that the HMM generated the sequence: " << p << ::std::endl;
96
Backward algorithm
97    ::std::cout << "Backward algorithm" << ::std::endl;
98    p = backwardAlgorithm(hmm, sequence);
99    ::std::cout << "Probability that the HMM generated the sequence: " << p << ::std::endl;
100
101    return 0;
102}
SeqAn - Sequence Analysis Library - www.seqan.de
 

Page built @2011/02/08 21:37:01