|
wolffd@0
|
1 function bnet = fixed_args_mk_motif_hhmm(motif_length, motif_pattern, background_char)
|
|
wolffd@0
|
2 %
|
|
wolffd@0
|
3 % BNET = MK_MOTIF_HHMM(MOTIF_LENGTH)
|
|
wolffd@0
|
4 % Make the following HHMM
|
|
wolffd@0
|
5 %
|
|
wolffd@0
|
6 % S2 <----------------------> S1
|
|
wolffd@0
|
7 % | |
|
|
wolffd@0
|
8 % | |
|
|
wolffd@0
|
9 % M1 -> M2 -> M3 -> end B1 -> end
|
|
wolffd@0
|
10 %
|
|
wolffd@0
|
11 % where Mi represents the i'th letter in the motif
|
|
wolffd@0
|
12 % and B is the background state.
|
|
wolffd@0
|
13 % Si chooses between running the motif or the background.
|
|
wolffd@0
|
14 % The Si and B states have self loops (not shown).
|
|
wolffd@0
|
15 %
|
|
wolffd@0
|
16 % The transition params are defined to respect the above topology.
|
|
wolffd@0
|
17 % The background is uniform; each motif state has a random obs. distribution.
|
|
wolffd@0
|
18 %
|
|
wolffd@0
|
19 % BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN)
|
|
wolffd@0
|
20 % In this case, we make the motif submodel deterministically
|
|
wolffd@0
|
21 % emit the motif pattern.
|
|
wolffd@0
|
22 %
|
|
wolffd@0
|
23 % BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN, BACKGROUND_CHAR)
|
|
wolffd@0
|
24 % In this case, we make the background submodel
|
|
wolffd@0
|
25 % deterministically emit the specified character (to make the pattern
|
|
wolffd@0
|
26 % easier to see).
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28 if nargin < 2, motif_pattern = []; end
|
|
wolffd@0
|
29 if nargin < 3, background_char = []; end
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 chars = ['a', 'c', 'g', 't'];
|
|
wolffd@0
|
32 Osize = length(chars);
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 motif_length = length(motif_pattern);
|
|
wolffd@0
|
35 Qsize = [2 motif_length];
|
|
wolffd@0
|
36 Qnodes = 1:2;
|
|
wolffd@0
|
37 D = 2;
|
|
wolffd@0
|
38 transprob = cell(1,D);
|
|
wolffd@0
|
39 termprob = cell(1,D);
|
|
wolffd@0
|
40 startprob = cell(1,D);
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 % startprob{d}(k,j), startprob{1}(1,j)
|
|
wolffd@0
|
43 % transprob{d}(i,k,j), transprob{1}(i,j)
|
|
wolffd@0
|
44 % termprob{d}(k,j)
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 % LEVEL 1
|
|
wolffd@0
|
48
|
|
wolffd@0
|
49 startprob{1} = zeros(1, 2);
|
|
wolffd@0
|
50 startprob{1} = [1 0]; % always start in the background model
|
|
wolffd@0
|
51
|
|
wolffd@0
|
52 % When in the background state, we stay there with high prob
|
|
wolffd@0
|
53 % When in the motif state, we immediately return to the background state.
|
|
wolffd@0
|
54 transprob{1} = [0.8 0.2;
|
|
wolffd@0
|
55 1.0 0.0];
|
|
wolffd@0
|
56
|
|
wolffd@0
|
57
|
|
wolffd@0
|
58 % LEVEL 2
|
|
wolffd@0
|
59 startprob{2} = 'leftstart'; % both submodels start in substate 1
|
|
wolffd@0
|
60 transprob{2} = zeros(motif_length, 2, motif_length);
|
|
wolffd@0
|
61 termprob{2} = zeros(2, motif_length);
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 % In the background model, we only use state 1.
|
|
wolffd@0
|
64 transprob{2}(1,1,1) = 1; % self loop
|
|
wolffd@0
|
65 termprob{2}(1,1) = 0.2; % prob transition to end state
|
|
wolffd@0
|
66
|
|
wolffd@0
|
67 % Motif model
|
|
wolffd@0
|
68 transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); % no self loops
|
|
wolffd@0
|
69 termprob{2}(2,end) = 1.0; % last state immediately terminates
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72 % OBS LEVEl
|
|
wolffd@0
|
73
|
|
wolffd@0
|
74 obsprob = zeros([Qsize Osize]);
|
|
wolffd@0
|
75 if isempty(background_char)
|
|
wolffd@0
|
76 % uniform background model
|
|
wolffd@0
|
77 obsprob(1,1,:) = normalise(ones(Osize,1));
|
|
wolffd@0
|
78 else
|
|
wolffd@0
|
79 % deterministic background model (easy to see!)
|
|
wolffd@0
|
80 m = find(chars==background_char);
|
|
wolffd@0
|
81 obsprob(1,1,m) = 1.0;
|
|
wolffd@0
|
82 end
|
|
wolffd@0
|
83
|
|
wolffd@0
|
84 if gen_motif
|
|
wolffd@0
|
85 % initialise with true motif (cheating)
|
|
wolffd@0
|
86 for i=1:motif_length
|
|
wolffd@0
|
87 m = find(chars == motif_pattern(i));
|
|
wolffd@0
|
88 obsprob(2,i,m) = 1.0;
|
|
wolffd@0
|
89 end
|
|
wolffd@0
|
90 else
|
|
wolffd@0
|
91 obsprob(2,:,:) = mk_stochastic(ones(motif_length, Osize));
|
|
wolffd@0
|
92 end
|
|
wolffd@0
|
93
|
|
wolffd@0
|
94 Oargs = {'CPT', obsprob};
|
|
wolffd@0
|
95
|
|
wolffd@0
|
96 [bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
|
|
wolffd@0
|
97 'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
|
|
wolffd@0
|
98 'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
|
|
wolffd@0
|
99
|
