diff options
| author | ziejd2 | 2017-09-28 15:04:40 -0500 |
|---|---|---|
| committer | ziejd2 | 2017-09-28 15:04:40 -0500 |
| commit | 8070dc963753142bb86c4ed698d91fd623ed28e7 (patch) | |
| tree | d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/BNT/examples/dynamic/HHMM/Motif | |
| parent | 7cc31810d53176e805532b2789955f4eedbce6bb (diff) | |
| download | BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz | |
BNW using Octave instead of Matlab.
This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning. I am calling this BNW_1.02. It can be accessed at: compbio.uthsc.edu/BNW_1.02
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 |
