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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
function [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(varargin)
% [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(...)
%
% Make the following HHMM
%
% S2 <----------------------> S1
% | |
% | |
% M1 -> M2 -> M3 -> end B1 -> end
%
% where Mi represents the i'th letter in the motif
% and B is the background state.
% Si chooses between running the motif or the background.
% The Si and B states have self loops (not shown).
%
% The transition params are defined to respect the above topology.
% The background is uniform; each motif state has a random obs. distribution.
%
% Optional params:
% motif_length - required, unless we specify motif_pattern
% motif_pattern - if specified, we make the motif submodel deterministically
% emit this pattern
% background - if specified, we make the background submodel
% deterministically emit this (makes the motif easier to see!)
args = varargin;
nargs = length(args);
% extract pattern, if any
motif_pattern = [];
for i=1:2:nargs
switch args{i},
case 'motif_pattern', motif_pattern = args{i+1};
end
end
% set defaults
motif_length = length(motif_pattern);
background_char = [];
% get params
for i=1:2:nargs
switch args{i},
case 'motif_length', motif_length = args{i+1};
case 'background', background_char = args{i+1};
end
end
chars = ['a', 'c', 'g', 't'];
Osize = length(chars);
Qsize = [2 motif_length];
Qnodes = 1:2;
D = 2;
transprob = cell(1,D);
termprob = cell(1,D);
startprob = cell(1,D);
% startprob{d}(k,j), startprob{1}(1,j)
% transprob{d}(i,k,j), transprob{1}(i,j)
% termprob{d}(k,j)
% LEVEL 1
startprob{1} = zeros(1, 2);
startprob{1} = [1 0]; % always start in the background model
% When in the background state, we stay there with high prob
% When in the motif state, we immediately return to the background state.
transprob{1} = [0.8 0.2;
1.0 0.0];
% LEVEL 2
startprob{2} = 'leftstart'; % both submodels start in substate 1
transprob{2} = zeros(motif_length, 2, motif_length);
termprob{2} = zeros(2, motif_length);
% In the background model, we only use state 1.
transprob{2}(1,1,1) = 1; % self loop
termprob{2}(1,1) = 0.2; % prob transition to end state
% Motif model
transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); % no self loops
termprob{2}(2,end) = 1.0; % last state immediately terminates
% OBS LEVEl
obsprob = zeros([Qsize Osize]);
if isempty(background_char)
% uniform background model
%obsprob(1,1,:) = normalise(ones(Osize,1));
obsprob(1,1,:) = normalise(rand(Osize,1));
else
% deterministic background model (easy to see!)
m = find(chars==background_char);
obsprob(1,1,m) = 1.0;
end
if ~isempty(motif_pattern)
% initialise with true motif (cheating)
for i=1:motif_length
m = find(chars == motif_pattern(i));
obsprob(2,i,m) = 1.0;
end
else
obsprob(2,:,:) = mk_stochastic(rand(motif_length, Osize));
end
if 0
Oargs = {'CPT', obsprob};
else
% We use a minent prior for the emission distribution for the states in the motif model
% (but not the background model). This encourages nearly deterministic distributions.
% We create an index matrix (where M = motif length)
% [2 1
% 2 2
% ...
% 2 M]
% and then convert this to a list of integers, which
% specifies when to use the minent prior (Q1=2 specifies motif model).
M = motif_length;
ndx = [2*ones(M,1) (1:M)'];
pcases = subv2ind([2 motif_length], ndx);
Oargs = {'CPT', obsprob, 'prior_type', 'entropic', 'entropic_pcases', pcases};
end
[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
'startprob', startprob, 'transprob', transprob, 'termprob', termprob);
|
