zClassifyPairs

 zClassifyPairs(File) calculates the rotation matrix, axis, angle, and shift
 between bases in File that are close enough to possibly be interacting, then
 classifies the interaction

0001 % zClassifyPairs(File) calculates the rotation matrix, axis, angle, and shift
0002 % between bases in File that are close enough to possibly be interacting, then
0003 % classifies the interaction
0004 
0005 function [File] = zClassifyPairs(File,Verbose)
0006 
0007 if isfield(File,'Pair'),
0008   File = rmfield(File,'Pair');                  % remove previous pair info
0009 end
0010 
0011 if nargin < 2,
0012   Verbose = 1;
0013 end
0014 
0015 if File.NumNT > 0,
0016 
0017 t = cputime;
0018 
0019 CL = zClassLimits;                              % read ClassLimits matrix
0020 
0021 if exist('PairExemplars.mat','file') > 0,
0022   load('PairExemplars','Exemplar');
0023 else
0024   Exemplar = [];
0025 end
0026 
0027 % -------- First screening of base pairs ------------------------------------
0028 
0029 DistCutoff = 10.5;                              % max distance for interaction
0030 [i,j] = find((File.Distance < DistCutoff).*(File.Distance > 0)); 
0031                                                 % screen by C-C distance
0032 k = find(i<j);                                  % look at each pair only once
0033 i = i(k);                                       % reduce list of indices
0034 j = j(k);                                       % reduce list of indices
0035 
0036 if Verbose > 0,
0037   fprintf('Classifying %5d pairs of bases for interactions ...', length(i));
0038 end
0039 
0040 % -------- Screen and analyze base pairs ------------------------------------
0041 % 1-AA  2-CA  3-GA  4-UA  5-AC  6-CC  7-GC  8-UC
0042 % 9-AG 10-CG 11-GG 12-UG 13-AU 14-CU 15-GU 16-UU
0043 
0044 pc = 1;                                         % index for pairs
0045 
0046 for k = 1:length(i),                            % loop through possible pairs
0047 
0048   Ni = File.NT(i(k));                           % nucleotide i information
0049   Nj = File.NT(j(k));                           % nucleotide j information
0050 
0051   [Pair,s,coplanar] = zClassifyPair(Ni,Nj,CL,Exemplar,0,Verbose);
0052 
0053   File.Coplanar(i(k),j(k)) = coplanar;
0054   File.Coplanar(j(k),i(k)) = coplanar;
0055 
0056   if ~isempty(Pair),
0057 
0058     if (s == 1),
0059       Pair.Base1Index = i(k);                       % bases in original order
0060       Pair.Base2Index = j(k);
0061       File.Edge(i(k),j(k)) =  Pair.Edge;
0062       File.Edge(j(k),i(k)) = -Pair.Edge;
0063     else
0064       Pair.Base1Index = j(k);                       % bases in reversed order
0065       Pair.Base2Index = i(k);
0066       File.Edge(i(k),j(k)) = -Pair.Edge;
0067       File.Edge(j(k),i(k)) =  Pair.Edge;
0068     end
0069 
0070     % --------------------------- code class to distinguish AA, CC, ... cases
0071 
0072     if (Ni.Code == Nj.Code) & (abs(Pair.Class) < 15),
0073       if i(k) > j(k),
0074         Pair.Class = -Pair.Class;                  
0075         % negative indicates that the lower-indexed base uses the dominant edge
0076       end
0077     end
0078 
0079     pc = pc + 1;                                    % increment pair counter
0080 
0081   end
0082 end   % loop over pairs
0083 
0084 if Verbose > 1,
0085   fprintf('Found %5d pairs that are possibly interacting\n', pc-1);
0086   fprintf('Classification took %4.2f minutes, or %4.0f classifications per minute\n', (cputime-t)/60, 60*(length(i))/(cputime-t));
0087 end
0088 
0089 end