Example Program
HMM
Hidden Markov Model code example
File "graph_hmm.cpp"
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() {

HMM creation

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 |

Begin state

23 | TVertexDescriptor begState = addVertex(hmm);
24 | assignBeginState(hmm, begState);
25 |

Add exon state

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 |

Add 5' splice site

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 |

Add intron state

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 |

End state

44 | TVertexDescriptor eState = addVertex(hmm);
45 | assignEndState(hmm, eState);
46 |

Transitions

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 |

Print the whole model

54 | ::std::cout << hmm << ::std::endl;
55 |

Viterbi algorithm

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 |

Forward algorithm

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 |

Backward algorithm

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 | }

See
Graph
SeqAn - Sequence Analysis Library - www.seqan.de

Page built @2012/03/14 16:34:45