| 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 |
| 6 | using namespace seqan;
| 7 |
| 8 | int main() {
|
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 |
|
23 | TVertexDescriptor begState = addVertex(hmm);
| 24 | assignBeginState(hmm, begState);
| 25 |
|
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 |
|
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 |
|
55 | TVertexDescriptor eState = addVertex(hmm);
| 56 | assignEndState(hmm, eState);
| 57 |
|
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 |
|
75 | ::std::cout << hmm << ::std::endl;
| 76 |
|
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 |
|
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 |
|
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 | }
|
|