% p_5S_LoopE_Motif sets up a model for the Loop E motif found in the 5S,
% plus closing canonical basepairs on either side
% It is for demonstration purposes.

clear Node

n=1;                                            % current node
  Node(n).type      = 'Initial';                % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Node(n).lpar      = 0.05;                     % left insertion parameter
  Node(n).rpar      = 0.05;                     % right insertion parameter

n=2;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(1,3,2,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=3;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(10,3,1,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=4;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(8,1,1,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=5;                                            % current node
  Node(n).type      = 'Motif';                  % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Node(n).P         = [0.05*ones(17,1) 0.95*ones(17,1)];
                                                % state to state transitions
  Node(n).PIns	    = [0.05 0.95];              % when no previous state
  Node(n).Left(1,:) = [1 2];                    % nucleotides to use on left
  Node(n).LIP       = [1];                      % probs for insertion possibs
  Node(n).Right(1,:)= [1];                      % nucleotides to use on right
  Node(n).RIP       = [1];                      % probs for insertion possibs
  Node(n).IBases(1,:)  = [1 2];
  Node(n).Score(:,:,1) = pIsoScore(-9,3,4);
  Node(n).IBases(2,:)  = [2 3];
  Node(n).Score(:,:,2) = pIsoScore(-4,4,1);

n=6;
  Node(n).type      = 'Initial';                % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Node(n).lpar      = 0.05;                     % left insertion parameter
  Node(n).rpar      = 0.05;                     % right insertion parameter

n=7;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(1,2,3,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=8;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(1,4,3,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=9;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(1,3,4,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=10;
  Node(n).type      = 'Basepair';               % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(1,3,2,0.01);
  Node(n).P         = ones(17,1)*Score;         % state to state transitions
  Node(n).PIns	    = Score;                    % when no previous state
  Node(n).lpar      = [0.01*ones(16,1); 0];     % left insertion parameters
  Node(n).rpar      = [0.01*ones(16,1); 0];     % right insertion parameters

n=11;                                           % current node
  Node(n).type      = 'Hairpin';                % node type
  Node(n).subtype   = '....';                   % pretend 'AAAA' loop
  Node(n).nextnode  = Inf;                      % index of next node in tree
  Node(n).P         = ones(17,1);               % state to state transitions
  Node(n).PIns      = 1;                        % when no previous state

% ----------------------------------------- End of user-specified part

% ----------------------------------------- Specify sequences to parse

Sequence(1).Organism = 'Viroid 1';
Sequence(1).X        = 'GGAGUACUGG....CUGGGAAAC';
Sequence(1).FastaNum = num2str(1);

Sequence(2).Organism = 'Viroid 2';
Sequence(2).X        = 'GGACUACUGG....CUGGGAAAC';
Sequence(2).FastaNum = num2str(2);

% ----------------------------------------- Parse sequences, display alignment

[Node,Sequence] = pParseSequences(Node,Sequence,1,length(Sequence));

pDisplayMultipleAlignment(Node,Sequence);

