zAnalyzePair

 zAnalyzePair(N1,N2,CL) computes distances, angles, and classification
 codes.

0001 % zAnalyzePair(N1,N2,CL) computes distances, angles, and classification
0002 % codes.
0003 
0004 function [Pair] = zAnalyzePair(N1,N2,CL,Exemplar,Displ,Verbose)
0005 
0006   if nargin < 3,
0007     CL = zClassLimits;                              % read ClassLimits matrix
0008 
0009     if exist('PairExemplars.mat','file') > 0,
0010       load('PairExemplars','Exemplar');
0011     else
0012       Exemplar = [];
0013     end
0014   end
0015 
0016   if nargin < 5,
0017     Displ = (N2.Fit(1,:)-N1.Fit(1,:)) * N1.Rot;   % vector shift from 1 to 2
0018   end
0019 
0020   Pair.Paircode = 4*(N2.Code-1) + N1.Code;     % AA is 1, CA is 2, etc.
0021 
0022   Pair.Displ = Displ;                        % vector shift (from 1 to 2)
0023 
0024   ro = N1.Rot'*N2.Rot;                       % rotation matrix from 1 to 2
0025   Pair.Normal = ro(:,3)';                    % normal to second plane
0026 
0027   if ro(3,3) > 0,                            % depending on orientation of 2,
0028     [ax,ang] = zAxisAngle(ro);               % rotation angle without a flip
0029   else
0030     [ax,ang] = zAxisAngle(ro*diag([-1 1 -1])); % flip base 2 first
0031   end
0032 
0033   Pair.Rot      = ro;
0034   Pair.RotAx    = ax';
0035   Pair.Ang      = ang;
0036 
0037   Pair.PlaneAng = acos(abs(N1.Rot(:,3)'*N2.Rot(:,3)))*57.29577951308232; 
0038                                              % angle between planes
0039 
0040   Lim(2,:) = [15 13 16 12];     % total number of atoms, including hydrogen
0041 
0042   d = zDistance(N2.Fit(1:Lim(2,N2.Code),:), N1.Center); 
0043                                            % distances to base 1 center
0044   [y,m] = min(d);                          % identify the closest atom
0045   m = m(1);                                % in case of a tie, use the first
0046   Pair.Gap = N1.Rot(:,3)'*(N2.Fit(m,:)-N1.Center)';% height above plane of 1
0047 
0048   Pair.MinDist = min(min(zDistance(N1.Fit,N2.Fit)));
0049 
0050   a = zCheckCutoffs(Pair.Displ,Pair.Normal,Pair.Ang,Pair.Gap,CL(:,:,Pair.Paircode));
0051                                            % find possible classifications
0052 
0053   % ------------------------ check for coplanarity
0054 
0055   Pair.Coplanar = 0;                      % default value
0056   if (abs(Pair.Gap) < 2) && (Pair.MinDist < 5),
0057     v  = N1.Center - N2.Center;           % vector from center to center
0058     v  = v / norm(v);                     % normalize
0059 
0060     dot1 = abs(v * N1.Rot(:,3));          % to calculate angle: v and normal
0061     dot2 = abs(v * N2.Rot(:,3));
0062 
0063     yy = 0.5;
0064 
0065     if (dot1 < yy) && (dot2 < yy),         % angle > acos(yy) = 60 degrees
0066 
0067       d = zDistance(N1.Fit(1:Lim(2,N1.Code),:), N2.Center); 
0068                                            % distances to base 2 center
0069       [y,m] = min(d);                      % identify the closest atom
0070       m = m(1);                            % in case of a tie, use the first
0071       Gap2 = N2.Rot(:,3)'*(N1.Fit(m,:)-N2.Center)';% height above plane of 1
0072 
0073       if abs(Gap2) < 2,
0074         Pair.Coplanar = min([(2-abs(Pair.Gap))/2 (2-abs(Gap2))/2 (yy-dot1)/yy (yy-dot2)/yy min(1,5-Pair.MinDist)]);
0075 
0076       end
0077     end
0078   end
0079 
0080   % ---------- Notify and remove multiple classifications
0081 
0082   if length(a) > 1,
0083     if max(fix(a)) > min(fix(a)),              % different integer parts
0084       if Verbose > 1,
0085         fprintf('Bases %1s%5s(%1s) and %1s%5s(%1s) fall into categories ', N1.Base, N1.Number, N1.Chain, N2.Base, N2.Number, N2.Chain);
0086         for k=1:length(a),
0087           fprintf('%6.2f ',a(k));
0088         end
0089         fprintf('\n');
0090       end
0091       a = a(1);
0092     else
0093       a = sign(a(1))*min(abs(a));               % use primary version of class
0094     end
0095   end
0096 
0097   % ---------- Calculate hydrogen bonds for base pairing interactions
0098 
0099   if (abs(a) < 14),                           % if it appears to be a basepair
0100     Pair.Hydrogen = zCheckHydrogen(N1,N2,a);
0101   else
0102     Pair.Hydrogen = [];
0103   end
0104 
0105   % ---------- Eliminate out of plane interactions
0106 
0107   if (abs(a) < 11) || ((abs(a) >= 13) && (abs(a) < 14)),   % avoid cSS, tSS
0108     if length(Pair.Hydrogen) > 0,
0109       if abs(Pair.Gap) > 2.0,
0110         a = 30.1;
0111       end
0112     elseif abs(Pair.Gap) > 1.6,
0113       a = 30.2;
0114     end
0115   end
0116 
0117   % ---------- Eliminate bad hydrogen bonds
0118 
0119   if (abs(a)<14) && (abs(a) - fix(abs(a)) < 0.5),% standard planar interaction
0120     if (length(Pair.Hydrogen) > 0),             % with hydrogens to check
0121       goodhydrogens = 0;
0122       for h = 1:length(Pair.Hydrogen),
0123         if isempty(Pair.Hydrogen(h).Angle),     % no third atom available
0124           if Pair.Hydrogen(h).Distance <= 4.5,  %
0125             goodhydrogens = goodhydrogens + 1;
0126           end
0127         else 
0128           if (Pair.Hydrogen(h).Angle >= 110)&&(Pair.Hydrogen(h).Distance <= 4),
0129             goodhydrogens = goodhydrogens + 1;
0130           end
0131         end
0132       end
0133 
0134       if goodhydrogens < length(Pair.Hydrogen),  % missing hydrogen bond
0135         gh = length(Pair.Hydrogen);              % required number
0136 
0137         if (length(Pair.Hydrogen) == 4),
0138           gh = 3;
0139         elseif (length(Pair.Hydrogen) == 3),
0140           gh = 2;
0141         end
0142 
0143         % Implement a few exceptions to hydrogen bonding rules
0144 
0145         if (fix(abs(a)) == 13) && (Pair.Paircode == 1)
0146           gh = 0;
0147         elseif (fix(a) == 13) && (Pair.Paircode == 5)
0148           gh = 0;
0149         elseif (fix(a) == -13) && (Pair.Paircode == 5)
0150           gh = 1;
0151         elseif (fix(abs(a)) == 13) && (Pair.Paircode == 6)
0152           gh = 1;
0153         elseif (fix(a) == 11) && (Pair.Paircode == 6)
0154           gh = 0;
0155         end
0156      
0157         if goodhydrogens < gh,
0158           a = 30.3;                              % reject this pair
0159         end
0160       end
0161     else
0162 
0163     end
0164   end
0165 
0166   % ---------- Measure stacking overlap
0167 
0168   SO1 = zStackingOverlap(N1,N2);
0169   SO2 = zStackingOverlap(N2,N1);
0170 
0171   if (SO1 > 0) && (SO2 > 0),
0172     Pair.StackingOverlap = (SO1+SO2)/2;
0173   else
0174     Pair.StackingOverlap = 0;
0175   end
0176   
0177   % ----------- Is an unclassified pair really stacked?
0178 
0179   if (fix(a) == 30) && (Pair.StackingOverlap > 0) && (Pair.MinDist < 4),
0180     if Pair.Displ(3) > 0,
0181       if Pair.Normal(3) > 0,
0182         a = 21;                       % second base above, pointing up
0183       else
0184         a = 22;                       % second base above, pointing down
0185       end
0186     else
0187       if Pair.Normal(3) > 0,
0188         a = -21;                      % second base below, pointing up
0189       else
0190         a = 23;                       % second base below, pointing down
0191       end
0192     end
0193   end
0194 
0195   % ----------- If classified as stacking, is there really overlap?
0196 
0197   if (fix(a) >= 21) && (fix(a) < 24) && (Pair.StackingOverlap == 0),
0198     a = 30;
0199 %    if N1.Loc(1,1) < N2.Loc(1,1),     % only print each pair once
0200       % fprintf('Bases %1s%5s(%1s) and %1s%5s(%1s) have no stacking overlap\n', N1.Base, N1.Number, N1.Chain, N2.Base, N2.Number, N2.Chain);
0201 %    end
0202   end
0203 
0204   % -------------------------- find distance to nearest exemplar
0205 
0206   if ~isempty(Exemplar),
0207     [c,d,ff,gg,h] = zDistanceToExemplars(Exemplar,N1,N2);
0208     Pair.Classes   = c(1:3);
0209     Pair.Distances = d(1:3);
0210   else
0211     Pair.Classes   = 99 * ones(1,3);
0212     Pair.Distances = 99999999 * ones(1,3);
0213   end
0214 
0215   % --------------------------- check hydrogen bonds if nearest exemplar
0216   % --------- is a basepair, even if this pair is far from that exemplar
0217 
0218   if (abs(a) >= 14) && (abs(Pair.Classes(1)) < 14),
0219     Pair.Hydrogen = zCheckHydrogen(N1,N2,Pair.Classes(1));
0220   end
0221 
0222   % ------------------------ record nearest exemplar if no classification yet
0223 
0224   if (fix(a) == 30),
0225     if (Pair.Distances(1) < 0.8) && (Pair.Normal(3) * Exemplar(ff(1),gg(1)).R(3,3) > 0), % same flip
0226       b = a-fix(a);                          % extract decimal code for reason
0227       c = Pair.Classes(1);
0228       a = sign(c) * (100 + abs(c) + b/1000);
0229     elseif (Pair.Distances(2) < 0.8) && (Pair.Normal(3) * Exemplar(ff(2),gg(2)).R(3,3) > 0), % same flip
0230       b = a-fix(a);                          % extract decimal code for reason
0231       c = Pair.Classes(2);
0232       a = sign(c) * (100 + abs(c) + b/1000);
0233     elseif (Pair.Distances(3) < 0.8) && (Pair.Normal(3) * Exemplar(ff(3),gg(3)).R(3,3) > 0), % same flip
0234       b = a-fix(a);                          % extract decimal code for reason
0235       c = Pair.Classes(3);
0236       a = sign(c) * (100 + abs(c) + b/1000);
0237     end
0238   end
0239 
0240   % ------------------------ store the classification
0241 
0242   Pair.Class = a;
0243 
0244   % ------------------------ store the edge information
0245 
0246   % reverse classification for GC and CG pairs, but why, exactly????
0247 
0248   if ((Pair.Paircode == 7) || (Pair.Paircode == 10)) && (abs(Pair.Class) < 14),
0249     Pair.Edge = -Pair.Class;
0250   else
0251     Pair.Edge = Pair.Class;
0252   end
0253 
0254   Pair.EdgeText = zEdgeText(Pair.Edge,1,Pair.Paircode);  % ever used?
0255 
0256