| Example Program HMM Silent States Hidden Markov Model with silent states code example File "graph_hmm_silent.cpp" 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 |
6 | using namespace seqan;
7 |
8 | int 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 | }
See Graph SeqAn - Sequence Analysis Library - www.seqan.de