| Example Program HMM Hidden Markov Model code example A tutorial about HMMs.
| 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, Default()> > 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 exonState = addVertex(hmm);
| | 27 | emissionProbability(hmm, exonState, dnaA) = 0.25;
| | 28 | emissionProbability(hmm, exonState, dnaC) = 0.25;
| | 29 | emissionProbability(hmm, exonState, dnaG) = 0.25;
| | 30 | emissionProbability(hmm, exonState, dnaT) = 0.25;
| | 31 |
|
| 32 | TVertexDescriptor spliceState = addVertex(hmm);
| | 33 | emissionProbability(hmm, spliceState, dnaA) = 0.05;
| | 34 | emissionProbability(hmm, spliceState, dnaC) = 0.0;
| | 35 | emissionProbability(hmm, spliceState, dnaG) = 0.95;
| | 36 | emissionProbability(hmm, spliceState, dnaT) = 0.0;
| | 37 |
|
| 38 | TVertexDescriptor intronState = addVertex(hmm);
| | 39 | emissionProbability(hmm, intronState, dnaA) = 0.4;
| | 40 | emissionProbability(hmm, intronState, dnaC) = 0.1;
| | 41 | emissionProbability(hmm, intronState, dnaG) = 0.1;
| | 42 | emissionProbability(hmm, intronState, dnaT) = 0.4;
| | 43 |
|
| 44 | TVertexDescriptor eState = addVertex(hmm);
| | 45 | assignEndState(hmm, eState);
| | 46 |
|
| 47 | addEdge(hmm, exonState, exonState, 0.9);
| | 48 | addEdge(hmm, exonState, spliceState, 0.1);
| | 49 | addEdge(hmm, spliceState, intronState, 1.0);
| | 50 | addEdge(hmm, begState, exonState, 1.0);
| | 51 | addEdge(hmm, intronState, intronState, 0.9);
| | 52 | addEdge(hmm, intronState, eState, 0.1);
| | 53 |
|
| 54 | ::std::cout << hmm << ::std::endl;
| | 55 |
|
| 56 | String<Dna> sequence = "CTTCATGTGAAAGCAGACGTAAGTCA";
| | 57 | String<TVertexDescriptor> path;
| | 58 | TProbability p = viterbiAlgorithm(hmm, sequence, path);
| | 59 | ::std::cout << "Viterbi algorithm" << ::std::endl;
| | 60 | ::std::cout << "Probability of best path: " << p << ::std::endl;
| | 61 | ::std::cout << "Sequence: " << ::std::endl;
| | 62 | for(TSize i = 0; i<length(sequence); ++i) ::std::cout << sequence[i] << ',';
| | 63 | ::std::cout << ::std::endl;
| | 64 | ::std::cout << "State path: " << ::std::endl;
| | 65 | for(TSize i = 0; i<length(path); ++i) {
| | 66 | ::std::cout << path[i];
| | 67 | if (isSilent(hmm, path[i])) ::std::cout << " (Silent)";
| | 68 | if (i < length(path) - 1) ::std::cout << ',';
| | 69 | }
| | 70 | ::std::cout << ::std::endl;
| | 71 |
|
| 72 | ::std::cout << "Forward algorithm" << ::std::endl;
| | 73 | p = forwardAlgorithm(hmm, sequence);
| | 74 | ::std::cout << "Probability that the HMM generated the sequence: " << p << ::std::endl;
| | 75 |
|
| 76 | ::std::cout << "Backward algorithm" << ::std::endl;
| | 77 | p = backwardAlgorithm(hmm, sequence);
| | 78 | ::std::cout << "Probability that the HMM generated the sequence: " << p << ::std::endl;
| | 79 |
| | 80 | return 0;
| | 81 | }
|
|