diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif')
7 files changed, 328 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries new file mode 100644 index 00000000..93258c2a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Entries @@ -0,0 +1,5 @@ +/fixed_args_mk_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/learn_motif_hhmm.m/1.1.1.1/Tue Jul 2 22:56:14 2002// +/mk_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +/sample_motif_hhmm.m/1.1.1.1/Wed May 29 15:59:54 2002// +D diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository new file mode 100644 index 00000000..5062cfd4 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/examples/dynamic/HHMM/Motif diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m new file mode 100644 index 00000000..ce4e288b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/fixed_args_mk_motif_hhmm.m @@ -0,0 +1,99 @@ +function bnet = fixed_args_mk_motif_hhmm(motif_length, motif_pattern, background_char) +% +% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH) +% 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. +% +% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN) +% In this case, we make the motif submodel deterministically +% emit the motif pattern. +% +% BNET = MK_MOTIF_HHMM(MOTIF_LENGTH, MOTIF_PATTERN, BACKGROUND_CHAR) +% In this case, we make the background submodel +% deterministically emit the specified character (to make the pattern +% easier to see). + +if nargin < 2, motif_pattern = []; end +if nargin < 3, background_char = []; end + +chars = ['a', 'c', 'g', 't']; +Osize = length(chars); + +motif_length = length(motif_pattern); +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)); +else + % deterministic background model (easy to see!) + m = find(chars==background_char); + obsprob(1,1,m) = 1.0; +end + +if gen_motif + % 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(ones(motif_length, Osize)); +end + +Oargs = {'CPT', obsprob}; + +[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); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m new file mode 100644 index 00000000..54553460 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/learn_motif_hhmm.m @@ -0,0 +1,75 @@ + +seed = 0; +rand('state', seed); +randn('state', seed); + +chars = ['a', 'c', 'g', 't']; +motif = 'accca'; +motif_length = length(motif); +motif_code = zeros(1, motif_length); +for i=1:motif_length + motif_code(i) = find(chars == motif(i)); +end + +[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_length', length(motif)); +%[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_pattern', motif); +ss = bnet_init.nnodes_per_slice; + + + +% We generate a training set by creating uniform sequences, +% and inserting a single motif at a random location. +ntrain = 100; +T = 20; +cases = cell(1, ntrain); + +if 1 + % uniform background + background_dist = normalise(ones(1, length(chars))); +end +if 0 + % use a constant background + background_dist = zeros(1, length(chars)); + m = find(chars=='t'); + background_dist(m) = 1.0; +end +if 0 + % use a background skewed away from the motif + p = 0.01; q = (1-(2*p))/2; + background_dist = [p p q q]; +end + +unif_pos = normalise(ones(1, T-length(motif))); +cases = cell(1, ntrain); +data = zeros(1,T); +for i=1:ntrain + data = sample_discrete(background_dist, 1, T); + L = sample_discrete(unif_pos, 1, 1); + data(L:L+length(motif)-1) = motif_code; + cases{i} = cell(ss, T); + cases{i}(Onode,:) = num2cell(data); +end +disp('sample training cases') +for i=1:5 + chars(cell2num(cases{i}(Onode,:))) +end + +engine_init = hmm_inf_engine(bnet_init); + +[bnet_learned, LL, engine_learned] = ... + learn_params_dbn_em(engine_init, cases, 'max_iter', 100, 'thresh', 1e-2); +% 'anneal', 1, 'anneal_rate', 0.7); + +% extract the learned motif profile +eclass = bnet_learned.equiv_class; +CPDO=struct(bnet_learned.CPD{eclass(Onode,1)}); +fprintf('columns = chars, rows = states\n'); +profile_learned = squeeze(CPDO.CPT(2,:,:)) +[m,ndx] = max(profile_learned, [], 2); +map_motif_learned = chars(ndx) +back_learned = squeeze(CPDO.CPT(1,1,:))' +%map_back_learned = chars(argmax(back_learned)) + +CPDO_init = struct(bnet_init.CPD{eclass(Onode,1)}); +profile_init = squeeze(CPDO_init.CPT(2,:,:)); +back_init = squeeze(CPDO_init.CPT(1,1,:))'; diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m new file mode 100644 index 00000000..32980397 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m @@ -0,0 +1,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); + diff --git a/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m new file mode 100644 index 00000000..b5822e10 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif/sample_motif_hhmm.m @@ -0,0 +1,10 @@ +%bnet = mk_motif_hhmm('motif_pattern', 'acca', 'background', 't'); +bnet = mk_motif_hhmm('motif_pattern', 'accaggggga', 'background', []); + +chars = ['a', 'c', 'g', 't']; +Tmax = 100; + +for seqi=1:5 + evidence = cell2num(sample_dbn(bnet, 'length', Tmax)); + chars(evidence(end,:)) +end |