1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
//![includes]
#include <iostream>
#include <seqan/graph_algorithms.h>
using namespace seqan2;
//![includes]
//![typedefs]
int main()
{
typedef LogProb<> TProbability;
typedef Dna TAlphabet;
typedef Size<TAlphabet>::Type TSize;
typedef Graph<Hmm<TAlphabet, TProbability, Default> > THmm;
typedef VertexDescriptor<THmm>::Type TVertexDescriptor;
//![typedefs]
//![variables]
Dna dnaA = Dna('A');
Dna dnaC = Dna('C');
Dna dnaG = Dna('G');
Dna dnaT = Dna('T');
THmm hmm;
//![variables]
//![begin-state]
TVertexDescriptor begState = addVertex(hmm);
assignBeginState(hmm, begState);
//![begin-state]
//![main-states-emissions]
TVertexDescriptor exonState = addVertex(hmm);
emissionProbability(hmm, exonState, dnaA) = 0.25;
emissionProbability(hmm, exonState, dnaC) = 0.25;
emissionProbability(hmm, exonState, dnaG) = 0.25;
emissionProbability(hmm, exonState, dnaT) = 0.25;
TVertexDescriptor spliceState = addVertex(hmm);
emissionProbability(hmm, spliceState, dnaA) = 0.05;
emissionProbability(hmm, spliceState, dnaC) = 0.0;
emissionProbability(hmm, spliceState, dnaG) = 0.95;
emissionProbability(hmm, spliceState, dnaT) = 0.0;
TVertexDescriptor intronState = addVertex(hmm);
emissionProbability(hmm, intronState, dnaA) = 0.4;
emissionProbability(hmm, intronState, dnaC) = 0.1;
emissionProbability(hmm, intronState, dnaG) = 0.1;
emissionProbability(hmm, intronState, dnaT) = 0.4;
//![main-states-emissions]
//![end-state]
TVertexDescriptor endState = addVertex(hmm);
assignEndState(hmm, endState);
//![end-state]
//![transitions]
addEdge(hmm, begState, exonState, 1.0);
addEdge(hmm, exonState, exonState, 0.9);
addEdge(hmm, exonState, spliceState, 0.1);
addEdge(hmm, spliceState, intronState, 1.0);
addEdge(hmm, intronState, intronState, 0.9);
addEdge(hmm, intronState, endState, 0.1);
//![transitions]
std::cout << "//![print-model]" << std::endl;
//![print-model]
std::cout << hmm << std::endl;
//![print-model]
std::cout << "//![print-model]" << std::endl;
std::cout << "//![viterbi]" << std::endl;
//![viterbi]
String<Dna> sequence = "CTTCATGTGAAAGCAGACGTAAGTCA";
String<TVertexDescriptor> path;
TProbability p = viterbiAlgorithm(path, hmm, sequence);
std::cout << "Viterbi algorithm" << std::endl;
std::cout << "Probability of best path: " << p << std::endl;
std::cout << "Sequence: " << std::endl;
for (TSize i = 0; i < length(sequence); ++i)
std::cout << sequence[i] << ',';
std::cout << std::endl;
std::cout << "State path: " << std::endl;
for (TSize i = 0; i < length(path); ++i)
{
std::cout << path[i];
if (isSilent(hmm, path[i]))
std::cout << " (Silent)";
if (i < length(path) - 1)
std::cout << ',';
}
std::cout << std::endl;
//![viterbi]
std::cout << "//![viterbi]" << std::endl;
std::cout << "//![forward-backward]" << std::endl;
//![forward-algorithm]
std::cout << "Forward algorithm" << std::endl;
p = forwardAlgorithm(hmm, sequence);
std::cout << "Probability that the HMM generated the sequence: " << p << std::endl;
//![forward-algorithm]
//![backward-algorithm]
std::cout << "Backward algorithm" << std::endl;
p = backwardAlgorithm(hmm, sequence);
std::cout << "Probability that the HMM generated the sequence: " << p << std::endl;
//![backward-algorithm]
std::cout << "//![forward-backward]" << std::endl;
//![backward-algorithm]
return 0;
}
//![backward-algorithm]
|
