zPhosphateInteractions

 zPhosphateInteractions checks all nearby pairs of bases for base-phosphate
 interactions, and stores them in a sparse matrix field BasePhosphate

0001 % zPhosphateInteractions checks all nearby pairs of bases for base-phosphate
0002 % interactions, and stores them in a sparse matrix field BasePhosphate
0003 
0004 function [File,D] = zPhosphateInteractions(File,Verbose,Self)
0005 
0006 if nargin < 2,
0007   Verbose = 0;
0008 end
0009 
0010 if nargin < 3,
0011   Self = 1;
0012 end
0013 
0014 t = cputime;
0015 
0016 % ------------------------------  Meanings of the classification codes:
0017     %  1  A C2-H2  interacts with oxygen of phosphate, called 2BPh
0018     %  2  A N6-1H6 interacts with oxygen of phosphate, called 6BPh
0019     %  3  A N6-2H6 interacts with oxygen of phosphate, called 7BPh
0020     %  4  A C8-H8  interacts with oxygen of phosphate, called 0BPh
0021     %  5  C N4-2H4 interacts with oxygen of phosphate, called 6BPh
0022     %  6  C N4-1H4 interacts with oxygen of phosphate, called 7BPh
0023     %  7  C N4-1H4 and C5-H5 interact with 2 oxygens of phosphate, called 8BPh
0024     % 18  C N4-1H4 and C5-H5 interact with just one oxygen, called 7BPh
0025     %  8  C C5-H5  interacts with oxygen of phosphate, called 9BPh
0026     %  9  C C6-H6  interacts with oxygen of phosphate, called 0BPh
0027     % 10  G N2-1H2 interacts with oxygen of phosphate, called 1BPh
0028     % 11  G N2-2H2 interacts with oxygen of phosphate, called 3BPh
0029     % 12  G N2-2H2 and N1-H1 interacts with 2 oxygens of phosphate, called 4BPh
0030     % 19  G N2-2H2 and N1-H2 interact with just one oxygen, called 4BPh
0031     % 13  G N1-H1  interacts with oxygen of phosphate, called 5BPh
0032     % 14  G C8-H8  interacts with oxygen of phosphate, called 0BPh
0033     % 15  U N3-H3  interacts with oxygen of phosphate, called 5BPh
0034     % 16  U C5-H5  interacts with oxygen of phosphate, called 9BPh
0035     % 17  U C6-H6  interacts with oxygen of phosphate, called 0BPh
0036 
0037     % Categories 18 and 19 are internal to FR3D.  They are displayed 8BPh, 4BPh
0038 
0039 % ------------------------------ temporary cutoffs!
0040 
0041 CarbonDist    = 4.0;                           % max massive - oxygen distance
0042 nCarbonDist   = 5;                           % near category
0043 
0044 NitrogenDist  = 3.5;                           % max massive - oxygen distance
0045 nNitrogenDist = 5.0;                           % near category
0046 
0047 AL = 130;                                      % angle limit for BPh
0048 nAL = 90;                                     % angle limit for nBPh
0049 
0050 DL([4 6 11]) = NitrogenDist;
0051 DL([7 8 9])  = CarbonDist;
0052 
0053 nDL([4 6 11]) = nNitrogenDist;
0054 nDL([7 8 9])  = nCarbonDist;
0055 
0056 % ------------------------------ temporary cutoffs!
0057 
0058 CarbonDist    = 4.0;                           % max massive - oxygen distance
0059 nCarbonDist   = 5;                           % near category
0060 
0061 NitrogenDist  = 3.5;                           % max massive - oxygen distance
0062 nNitrogenDist = 5.0;                           % near category
0063 
0064 AL = 130;                                      % angle limit for BPh
0065 nAL = 90;                                     % angle limit for nBPh
0066 
0067 DL([4 6 11]) = NitrogenDist;
0068 DL([7 8 9])  = CarbonDist;
0069 
0070 nDL([4 6 11]) = nNitrogenDist;
0071 nDL([7 8 9])  = nCarbonDist;
0072 
0073 % ------------------------------ specify cutoffs for classification
0074 
0075 CarbonDist    = 4.0;                           % max massive - oxygen distance
0076 nCarbonDist   = 4.5;                           % near category
0077 
0078 NitrogenDist  = 3.5;                           % max massive - oxygen distance
0079 nNitrogenDist = 4.0;                           % near category
0080 
0081 AL = 130;                                      % angle limit for BPh
0082 nAL = 110;                                     % angle limit for nBPh
0083 
0084 DL([4 6 11]) = NitrogenDist;
0085 DL([7 8 9])  = CarbonDist;
0086 
0087 nDL([4 6 11]) = nNitrogenDist;
0088 nDL([7 8 9])  = nCarbonDist;
0089 
0090 % ------------------------------- define basic data
0091 
0092 D = [];                          % where to save data if Verbose
0093 T = [];                          % data on which hydrogen with which oxygen(s)
0094 
0095 zStandardBases
0096 Sugar = {'C1*','C2*','O2*','C3*','O3*','C4*','O4*','C5*','O5*','P','O1P','O2P','O3 of next'};
0097 
0098 p   = [9 13 11 12];                           % rows of the phosphate oxygens
0099 pn  = {'O5*','O3*','O1P','O2P'};              % names of phosphate oxygens
0100 
0101 % ------------------------------ loop through files and classify
0102 
0103 for f = 1:length(File),
0104  if File(f).NumNT > 1,
0105   if isempty(File(f).Distance),
0106     c = cat(1,File(f).NT(1:File(f).NumNT).Center); % nucleotide centers
0107     File(f).Distance = zMutualDistance(c,16); % compute distances < 16 Angs
0108   end
0109 
0110   File(f).BasePhosphate = sparse(zeros(File(f).NumNT));
0111 
0112   % -------- First screening of base pairs -----------------------------------
0113 
0114   DistCutoff = 16;                              % max distance for interaction
0115   [i,j] = find((File(f).Distance < DistCutoff).*(File(f).Distance > 0)); 
0116                                                 % screen by C-C distance
0117 
0118   if Self > 0,
0119     i = [i; (1:length(File(f).NT))'];             % allow self interactions
0120     j = [j; (1:length(File(f).NT))'];             % allow self interactions
0121 %    i = [(1:length(File(f).NT))'; i];             % allow self interactions
0122 %    j = [(1:length(File(f).NT))'; i];             % allow self interactions
0123   end
0124 
0125   % -------- Screen and analyze pairs ----------------------------------------
0126 
0127   for k = 1:length(i),                          % loop through possible pairs
0128    N1 = File(f).NT(i(k));                       % nucleotide i information
0129    N2 = File(f).NT(j(k));                       % nucleotide j information
0130 
0131    DT = [];                                     % accumulate data here
0132 
0133    ph = (N2.Sugar(10,:)-N1.Fit(1,:)) * N1.Rot;  % phosphorus displacement
0134    if abs(ph(3)) < 4.5,                         % phosphorus close to plane
0135 
0136     switch N1.Code
0137       case 1,                         % Base A
0138               h   = [11 12 14 15];    % rows of the base hydrogens
0139               hn  = {'H2','H8','1H6','2H6'}; % names of the base hydrogens
0140               m   = [ 9  7  6  6];    % rows of the corresponding massive atoms
0141               e   = [ 1  4  2  3];    % code for location of the interaction
0142       case 2,                         % Base C
0143               h   = [10 11 12 13];
0144               hn  = {'H6','H5','1H4','2H4'}; % names of the base hydrogens
0145               m   = [ 7  8  6  6];
0146               e   = [ 9  8  6  5];
0147       case 3,                         % Base G
0148               h   = [12 13 15 16];
0149               hn  = {'H1','H8','1H2','2H2'}; % names of the base hydrogens
0150               m   = [ 4  7 11 11];
0151               e   = [13 14 10 11];
0152       case 4,                         % Base U
0153               h   = [ 9 11 12];
0154               hn  = {'H5','H3','H6'}; % names of the base hydrogens
0155               m   = [ 8  4  7];
0156               e   = [16 15 17];
0157     end
0158 
0159     dis = zDistance(N1.Fit(m,:), N2.Sugar(p,:)); % distances between mass & O's
0160     nearDL = nDL(m)' * ones(1,4);     % limits to compare to
0161     dis = dis .* (dis < nearDL);      % massive-oxygen pairs close enough
0162 
0163     g = [];                           % internal classification number
0164     w = [];                           % which oxygen interacts
0165     Angle = [];
0166     Dist  = [];
0167 
0168     for mm = 1:length(m),             % massive atom to consider
0169      pp = find(dis(mm,:));            % oxygens close enough to consider
0170      for n = 1:length(pp),            % loop through potential oxygens
0171       Angle(n)=zAngle(N1.Fit(m(mm),:),N1.Fit(h(mm),:),N2.Sugar(p(pp(n)),:));
0172                                       % base massive - hydrogen - oxygen angle
0173       Dist(n) = dis(mm,pp(n));        % distance
0174      end
0175 
0176      PAngle = zAngle(N1.Fit(m(mm),:),N1.Fit(h(mm),:),N2.Sugar(10,:));
0177                                   % base massive - hydrogen - phosphorus angle
0178      PDist  = zDistance(N1.Fit(m(mm),:),N2.Sugar(10,:));        % distance
0179 
0180      [u,v] = min(-Angle+60*Dist);     % order by quality of potential bond
0181 
0182      for n = 1:length(pp),            % loop through potential oxygens
0183       if Angle(n) > nAL,              % good enough to be "near" base-phosph
0184 
0185         if ((Angle(n) > AL) && (Dist(n) < DL(m(mm)))) % true BPh pair
0186           g = [g e(mm)];              % assign a non-near class.
0187           w = [w pp(n)];              % record which oxygen it is
0188           T = [T; [f i(k) j(k) e(mm) pp(n)]];
0189         else
0190           g = [g e(mm) + 100];        % > 100 means "near"
0191           w = [w pp(n)];              % record which oxygen it is
0192         end
0193 
0194         if Verbose > 1,
0195           % store information for later display
0196           ox = (N2.Sugar(p(pp(n)),:)-N1.Fit(1,:)) * N1.Rot; % oxygen displ
0197 
0198           a = [f i(k) j(k) N1.Code g(end) mm pp(n) Angle(n) Dist(n) ox ph File(f).Distance(i(k),j(k)) (v(1)==n) -u(1) PAngle PDist str2num(N1.Number) str2num(N2.Number)];
0199 
0200 % [ g(end) v(1) == n]
0201 
0202           % Columns:
0203           % 1  file number
0204           % 2  index of base
0205           % 3  index of nucleotide using phosphate
0206           % 4  code of base
0207           % 5  classification number for this massive-oxygen pair
0208           % 6  which massive atom is interacting
0209           % 7  which oxygen is interacting
0210           % 8  angle of interaction, in degrees
0211           % 9  distance from massive atom to oxygen, in Angstroms
0212           %10  displacement of oxygen atom relative to C1' of base
0213           %13  displacement of phophorus atom relative to C1' of base
0214           %16  distance between centers of the two bases
0215           %17  1 if this is the best oxygen for this hydrogen, 0 otherwise
0216           %18  approximate quality of the hydrogen bond
0217           %19  angle made by massive, hydrogen, phosphorus
0218           %20  distance from massive to phosphorus
0219           %21  nucleotide number of base
0220           %22  nucleotide number of phosphate donor
0221           %23  (to be added below) classification of this interaction
0222 
0223           if length(a(1,:)) == 22,
0224             DT = [DT; a];                  % append data to data matrix
0225           end
0226 
0227 %          zOtherPhosphateInteractions
0228 
0229           if Verbose > 3,
0230             fprintf('%6s base %s%5s %3s BPcode %3d %4s phosphate donor %s%5s %13s length %6.2f angle %6.2f interaction %s', File(f).Filename, N1.Base, N1.Number, AtomNames{h(mm),N1.Code}, g(end), zBasePhosphateText(g(end)), N2.Base, N2.Number, Sugar{p(pp(n))}, dis(mm,pp(n)), Angle(n), zEdgeText(File(f).Edge(i(k),j(k))));
0231 
0232 
0233             if a(17) == 1,
0234               fprintf('Best\n');
0235             else
0236               fprintf('\n');
0237             end
0238           end
0239         end
0240 
0241       end
0242 
0243      end  % loop over potential oxygens
0244     end   % loop over massive atoms
0245 
0246      if length(g) > 0,
0247        if (min(g) < 100) && (max(g) > 100),
0248          w = w(find(g < 100));
0249          g = g(find(g < 100));               % remove near classifications
0250        end
0251        if length(g) > 1,                     % multiple bonds
0252          g = sort(g);
0253          if (g(1) == 6) && (g(end) == 8) && (min(w) < max(w)),
0254            g = 7;                            % two different oxygens
0255          elseif (g(1) == 6) && (g(end) == 8) && (min(w) == max(w)),
0256            g = 18;
0257          elseif (g(1) == 11) && (g(end) == 13) && (min(w) < max(w)),
0258            g = 12;                           % two bonds formed
0259          elseif (g(1) == 11) && (g(end) == 13) && (min(w) == max(w)),
0260            g = 19;                           % two bonds formed
0261          elseif (min(g) == max(g)),
0262            g = g(1);
0263          elseif (g(end) < 100) && (Verbose > 0),
0264            fprintf('zPhosphateInteractions: Another case to consider: ');
0265            fprintf('File %s Base index %4d Phosphate index %4d Classes ', File(f).Filename, i(k), j(k));
0266            for m = 1:length(g),
0267              fprintf('%d ', g(m));
0268            end
0269            fprintf('\n');
0270          end
0271        end
0272 
0273        File(f).BasePhosphate(i(k),j(k)) =   g(1);  % record classification
0274 
0275        if Verbose > 1,
0276          if ~isempty(DT),
0277            D = [D; [DT g(1)*ones(length(DT(:,1)),1)]];
0278          end 
0279        end
0280 
0281      end
0282 
0283    end    % if vertical displacement of phosphate is less than 6 Angstroms
0284   end     % loop over nucleotide pairs
0285  end      % if File.NumNT > 1
0286 end       % loop over files
0287 
0288 
0289 if Verbose > 1,
0290   fprintf('Classifying base-phosphate interactions took %8.2f minutes\n', (cputime-t)/60);
0291   zPhosDisplay(D,0);                     % display parameters graphically
0292 end
0293 
0294 
0295 
0296 if Verbose > 2,
0297   for i = 1:(length(T(:,1))-1),
0298     if (T(i,2) == T(i+1,2)) && (T(i,3) == T(i+1,3)),
0299       if min(T(i,4),T(i+1,4)) == 6 && max(T(i,4),T(i+1,4)) == 8,
0300         T(i,4)   = 7;
0301         T(i+1,4) = 7;
0302         T(i,5)   = 5+4*min(T(i,5),T(i+1,5))+max(T(i,5),T(i+1,5));
0303         T(i+1,5) = 5+4*min(T(i,5),T(i+1,5))+max(T(i,5),T(i+1,5));
0304       elseif min(T(i,4),T(i+1,4)) == 11 && max(T(i,4),T(i+1,4)) == 13,
0305         T(i,4)   = 12;
0306         T(i+1,4) = 12;
0307         T(i,5)   = 5+4*min(T(i,5),T(i+1,5))+max(T(i,5),T(i+1,5));
0308         T(i+1,5) = 5+4*min(T(i,5),T(i+1,5))+max(T(i,5),T(i+1,5));
0309       end
0310     end
0311   end
0312 
0313   [Table,Chi1,P,Labels] = crosstab(T(:,4),T(:,5));
0314 end
0315 
0316 return
0317 
0318 % File = zAddNTData({'1s72','1j5e','2avy','2aw4','2j01'});
0319 % zPhosphateInteractions(File,3);
0320 
0321 Res = [];
0322 for f = 1:length(File);
0323   if isempty(File(f).Info.Resolution),
0324     Res(f) = 10;
0325   else
0326     Res(f) = File(f).Info.Resolution;
0327   end
0328 end
0329 File = File(find(Res <= 3.0));
0330 
0331 
0332 
0333 
0334 i = find(T(:,4) == 12);
0335 Search = [T(i,2) T(i,3) T(i,1)];
0336 xDisplayCandidates(File,Search)