% p_23S_1420_1444 is designed for a particular loop

Filename = 'p_23S_1420_1444.mat';

if exist(Filename,'file') > 0,            % load previous work, if it exists
  load(Filename);
  k = k + 1;                              % go to next sequence
  fprintf('Loaded previous computations\n');
else                                      % define model, load sequences

% ----------------------------------------- Basic parameters

% 1AA 2CA 3GA 4UA 5AC 6CC 7GC 8UC 9AG 10CG 11GG 12UG 13AU 14CU 15GU 16UU

% --------------------------------------- initialize model nodes ---------

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,2,3,0.02);    % slightly higher del prob
  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(1,2,3,0.012);
  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(1,4,1,0.011);
  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;
  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=6;
  Node(n).type      = 'Basepair';                % node type
  Node(n).nextnode  = n+1;                      % index of next node in tree
  Score             = pIsoScore(1,1,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=7;                                            % 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 4 5]; %3 is ignored because it may be bulged               % nucleotides to use on left
  Node(n).Left(2,:) = [1 2 3 4];                % nucleotides to use on left
  Node(n).LIP       = [0.9 0.1];                % probs for insertion possibs
  Node(n).Right(1,:)= [3 1];     %numbers in the opposite way, meaning the same       % nucleotides to use on right
  Node(n).Right(2,:)= [2 1];                    % nucleotides to use on right
  Node(n).RIP       = [0.9 0.1];                % probs for insertion possibs
%  Node(n).RIP        = 1;
  Node(n).IBases(1,:)  = [1 6];
  Node(n).Score(:,:,1) = pIsoScore(1,3,2);
  Node(n).IBases(2,:)  = [2 5];
  Node(n).Score(:,:,2) = pIsoScore(-4,2,1);
  Node(n).IBases(3,:)  = [3 6];
  Node(n).Score(:,:,3) = pIsoScore(5,2,2);
  Node(n).IBases(4,:)  = [4 5];
  Node(n).Score(:,:,4) = pIsoScore(1,4,1);
%%%%%%%%%%%%%%
n=8;
  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(1,1)*Score;    %1,1 because there a HP after it      % 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;                                            % current node
  Node(n).type      = 'Hairpin';                % node type
  Node(n).subtype   = '';    %put 'GNRA' if you have that loop, only this is coded now       % nothing special
  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

% ----------------------------------------- Read sequences

% KeyOrg = 'Haloarcula_marismortui';
% KeySeq = 'CCUCAGCACUGCUGAUCAAGCUGAGG';  % as in pdb file

%Sequence = pReadFASTA('2004_23sA_UNIQUE',KeyOrg,KeySeq);
  % use this form if you know the key organism and what sequence to look for
FastaFilename='C:\AliRes\_MatlabWork\MATLAB programs\Alignment_Analysis\FASTA_alignments\23S_ABE_2004_UNIQUE_Archaea';
Sequence = pReadFASTA(FastaFilename,4130,4158);
  % use this form if you know which columns of the FASTA file to use

% I cannot determine the appropriate columns to use in this FASTA file!
% So I am just using some columns and we'll see what we get

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

k = 1;
end

numtodo = 25;

for j=k:min(length(Sequence),(k+numtodo-1)),
  [Node,Sequence] = pParseSequences(Node,Sequence,j,j);
  k = j;
  drawnow
  if mod(j,10) == 0,
    pDisplayMultipleAlignment(Node,Sequence);
    save(Filename,'Node','Sequence','k');
  end
end

save(Filename,'Node','Sequence','k');

return

for j=1:length(Sequence),
  fprintf('%s %s\n', Sequence(j).Fasta(4130:4158),Sequence(j).Organism);
end

%close %Ali (not needed anymore because I removed clf from pParseSequences)