xAddNucleotide

 xAddNucleotide adds a nucleotide according to the screen matrices S

0001 % xAddNucleotide adds a nucleotide according to the screen matrices S
0002 
0003 function [NewList,NewSS] = xAddNucleotide(Model,List,S,SS,r)
0004 
0005 L = 100000;
0006 
0007 [s,t] = size(List);
0008 
0009 NewList = uint16(zeros(L,t+1));
0010 NewSS   = zeros(L,1);
0011 
0012 A = S{r,1};                % for some reason, this is MUCH faster
0013 B = S{r,2};
0014 C = S{r,3};
0015 D = S{r,4};
0016 
0017 if r >  5,  E = S{r,5}; end
0018 if r >  6,  F = S{r,6}; end
0019 if r >  7,  G = S{r,7}; end
0020 if r >  8,  H = S{r,8}; end
0021 if r >  9,  I = S{r,9}; end
0022 
0023 c = 0;                                     % number of candidates
0024 
0025 for n=1:s,
0026   a = List(n,1:(r-1));                     % a is the ith candidate
0027   switch r                                 % search for one more nucleotide
0028     case 5,
0029       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4));
0030     case 6,
0031       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4)) .* E(:,a(5));
0032     case 7,
0033       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4)) .* E(:,a(5)) .*  ...
0034           F(:,a(6));
0035     case 8,
0036       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4)) .* E(:,a(5)) .*  ...
0037           F(:,a(6)) .* G(:,a(7));
0038     case 9,
0039       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4)) .* E(:,a(5)) .*  ...
0040           F(:,a(6)) .* G(:,a(7)) .* H(:,a(8));
0041     case 10,
0042       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4)) .* E(:,a(5)) .*  ...
0043           F(:,a(6)) .* G(:,a(7)) .* H(:,a(8)) .* I(:,a(9));
0044     otherwise
0045       h = A(:,a(1)) .* B(:,a(2)) .* C(:,a(3)) .* D(:,a(4)) .* E(:,a(5)) .*  ...
0046           F(:,a(6)) .* G(:,a(7)) .* H(:,a(8)) .* I(:,a(9));
0047       for q = 10:(r-1), 
0048         h = h .* S{r,q}(:,a(q));           % slower than hard-coding
0049       end
0050   end
0051 
0052   m = find(h);                             % indices of new nucleotides
0053 
0054   if c + length(m) > L,                    % memory management; faster
0055     NewList = [NewList; uint16(zeros(L,t+1))];
0056     NewSS   = [NewSS; zeros(L,1)];
0057     L = 2*L;
0058   end
0059 
0060   for i=1:length(m),                       % append these new candidates
0061     c = c + 1;
0062     NewList(c,:) = [a m(i)];
0063     NewSS(c,1)   = SS(n,1);
0064   end
0065 end
0066 
0067 NewList = NewList(1:c,:);
0068 NewSS   = NewSS(1:c,:);
0069 
0070 if Model.Geometric > 0,
0071   Acceptable = uint16(zeros(size(NewList(:,1))));
0072   for n = 1:length(NewList(:,1)),
0073     a = NewList(n,:);
0074     b = a(r);
0075     y = NewSS(n,1) + A(b,a(1)) + B(b,a(2)) + C(b,a(3)) + D(b,a(4));
0076     switch r,
0077       case  6, y = y + E(b,a(5));
0078       case  7, y = y + E(b,a(5)) + F(b,a(6));
0079       case  8, y = y + E(b,a(5)) + F(b,a(6)) + G(b,a(7));
0080       case  9, y = y + E(b,a(5)) + F(b,a(6)) + G(b,a(7)) + H(b,a(8));
0081       case 10, y = y + E(b,a(5)) + F(b,a(6)) + G(b,a(7)) + H(b,a(8)) + ...
0082                        I(b,a(9));
0083       case 11, 
0084         y = y + E(b,a(5)) + F(b,a(6)) + G(b,a(7)) + H(b,a(8)) + I(b,a(9));
0085         for q = 10:(r-1),
0086           y = y + S{r,q}(b,a(q));
0087         end
0088     end
0089 
0090     if y < Model.SSCutoff(r),
0091       Acceptable(n) = 1;
0092     end
0093 
0094   end
0095 
0096   OK      = find(Acceptable);  
0097   NewList = NewList(OK,:);
0098   NewSS   = NewSS(OK,:);
0099 
0100 %  fprintf('Retained %6d, rejected %5d candidates on nucleotide %1d\n', length(OK), length(Acceptable)-length(OK), r);
0101 end
0102