xScatterPairs

 xScatterPairs displays multiple 3d scatter plots of pair parameters, with a menu that allows one to color points in different orders, or display different nucleotides from the candidates

0001 % xScatterPairs displays multiple 3d scatter plots of pair parameters, with a menu that allows one to color points in different orders, or display different nucleotides from the candidates
0002 
0003 function [ppp] = xScatterPairs(Search,N1,N2,ViewParam,Param)
0004 
0005 warning off
0006 
0007 if nargin < 3,
0008   N1 = 1;
0009   N2 = 2;
0010 end
0011 
0012 % ---------------------------------------------
0013 % --------------------------------------------- Collect basic data
0014 
0015 Candidates = Search.Candidates;
0016 [L,t] = size(Candidates);
0017 N     = t - 1;                                     % number of nucleotides
0018 File  = Search.File;
0019 CL    = zClassLimits;                              % read ClassLimits matrix
0020 
0021 ppp = 1:L;                                         % default ordering
0022 
0023 if exist('PairExemplars.mat','file') > 0,
0024   load('PairExemplars','Exemplar');
0025 else
0026   Exemplar = [];
0027 end
0028 
0029 % --------------------------------------------- Provide memory space
0030 
0031 for i = 1:N,
0032   for j = 1:N,
0033     SavedPairs(i,j).Classified = 0;
0034   end
0035 end
0036 
0037 % --------------------------------------------- Set defaults
0038 
0039 if nargin < 4,
0040   ViewParam.Color  = 6;
0041   ViewParam.FigNum = 1;
0042   ViewParam.Normal = 0;
0043   ViewParam.ClassLimits = 1;
0044 end
0045 
0046 if nargin < 5,
0047   Param.Decimal    = 1;
0048 end
0049 
0050 %---------------------------------------------- Loop for the menu
0051 
0052 Stop       = 0;
0053 Recolor    = 1;
0054 Reclassify = 1;
0055 Replot     = 1;
0056 
0057 while Stop == 0,
0058 
0059 
0060 % --------------------------------------------- Analyze pairs of nucleotides
0061 
0062 if Reclassify > 0,
0063 
0064  if SavedPairs(N1,N2).Classified == 0,
0065 
0066   pc = 1;                                          % pair counter
0067   FirstBaseCodes = [];
0068   Paircode       = [];
0069   fprintf('Nucleotide %3d at the origin, %3d away\n', N1, N2);
0070 
0071   for k = 1:L,                                     % Loop through candidates
0072     f  = Candidates(k,N+1);                        % file number
0073     i1 = Candidates(k,N1);                         % first nucleotide index
0074     i2 = Candidates(k,N2);                         % second nucleotide index
0075 
0076     Ni = File(f).NT(i1);
0077     Nj = File(f).NT(i2);
0078 
0079     [p,s] = zClassifyPair(Ni,Nj,CL,Exemplar,1);
0080 
0081     p.NT1 = Ni;
0082     p.NT2 = Nj;
0083     p.Filename = File(f).Filename;
0084     p.C1pC1p   = zDistance(Ni.Sugar(1,:),Nj.Sugar(1,:));  % C1' - C1' distance
0085     p.Resol    = File(f).Info.Resolution;
0086     if isempty(p.Resol)
0087       p.Resol = NaN;
0088     end
0089 
0090     if ~isempty(p),
0091       Pair(pc) = p;                                % store this pair
0092       pc = pc + 1;                                 % increment pair counter
0093       FirstBaseCodes(Ni.Code) = 1;                 % this base was seen
0094       Paircode(p.Paircode)    = 1;                 % this paircode was seen
0095     end
0096 
0097   end
0098 
0099   modcode = find(FirstBaseCodes);                  % codes of first bases seen
0100   Paircodes = find(Paircode);
0101 
0102   SavedPairs(N1,N2).Pair = Pair;
0103   SavedPairs(N1,N2).modcode = modcode;
0104   SavedPairs(N1,N2).Paircodes = Paircodes;
0105 
0106  else
0107 
0108   Pair      = SavedPairs(N1,N2).Pair;
0109   modcode   = SavedPairs(N1,N2).modcode;
0110   Paircodes = SavedPairs(N1,N2).Paircodes;
0111   
0112  end
0113 
0114 end
0115 
0116 %---------------------------------------------- Set the colors for the points
0117 
0118 if Recolor > 0,
0119 
0120  for k = 1:length(Pair),                                   % Loop through pairs
0121   p = Pair(k);
0122 
0123   % -------- Round categories to nearest integer, if desired -------------
0124 
0125   if Param.Decimal == 1,                  % select with decimal places in categ
0126     pClass = p.Edge;
0127     eClass = p.Classes(1);
0128   else
0129     pClass = fix(p.Edge);
0130     eClass = fix(p.Classes(1));
0131   end
0132 
0133   switch ViewParam.Color,
0134     case  1, Color(k) = k;
0135     case {2, 3, 4}, Color(k) = p.Edge;
0136     case  5, Color(k) = p.Paircode;
0137     case  6, Color(k) = p.Ang;
0138     case  7, Color(k) = p.Gap;
0139     case  8, Color(k) = p.Normal(3);
0140     case  9, Color(k) = p.Distances(1);
0141     case 11, Color(k) = p.StackingOverlap;
0142     case 12, Color(k) = p.Edge;
0143     case 13, Color(k) = p.Displ(1);
0144     case 14, Color(k) = p.Displ(2);
0145     case 15, Color(k) = p.Displ(3);
0146     otherwise, Color(k) = 1;
0147   end
0148 
0149  end
0150 
0151   [y,i] = sort(Color);
0152   ppp   = i;                                    % permutation
0153   Pair  = Pair(i);
0154   Color = Color(i);
0155 
0156 end
0157 
0158 % --------------------------------------------- Set color axis
0159 
0160 switch ViewParam.Color,
0161   case 1, ColorAxis =  [min(Color) max(Color)];
0162   case 2, ColorAxis =  [1 12];
0163   case 3, ColorAxis =  [15 19];
0164   case 4, ColorAxis =  [-12 30];
0165   case 5, ColorAxis =  [1 16];
0166   case 6, ColorAxis =  [min(Color) max(Color)];      % color by angle
0167   case 7, ColorAxis =  [min(Color) max(Color)];
0168   case 8, ColorAxis =  [-1 3];
0169   case 9, ColorAxis =  [0 4];
0170   case 10, ColorAxis = [-12 12];
0171   case 11, ColorAxis = [0 10];
0172   case 12, ColorAxis = [min(Color) max(Color)];
0173   case 13, ColorAxis =  [min(Color) max(Color)];
0174   case 14, ColorAxis =  [min(Color) max(Color)];
0175   case 15, ColorAxis =  [min(Color) max(Color)];
0176   otherwise, ColorAxis =  [0 10];
0177 end
0178 
0179 %---------------------------------------------- Plot displacements
0180 
0181 if Replot > 0,
0182 
0183 figure(ViewParam.FigNum)
0184 clf
0185 figure(ViewParam.FigNum)
0186 clf
0187 
0188 %set(gcf,'Renderer','OpenGL');     % fast rotation
0189 %set(gcf,'Renderer','zbuffer')
0190 %set(gcf,'Renderer','painters')
0191 
0192 for k = 1:length(Pair),                              % Loop through pairs
0193   p      = Pair(k);
0194   c(k)   = Color(k);
0195   e(k,:) = p.Displ;
0196 
0197   if ViewParam.Normal == 1,
0198     v = p.Normal/3;                                      % add normal vector
0199     plot3([e(k,1) e(k,1)+v(1)], [e(k,2) e(k,2)+v(2)], [e(k,3) e(k,3)+v(3)], 'b');
0200     hold on
0201   end
0202 end
0203 
0204 scatter3(e(:,1),e(:,2),e(:,3),18*ones(size(c)),c,'filled')
0205 
0206 if max(modcode) == min(modcode),             % all have same first base
0207   zPlotStandardBase(modcode);                % plot base at the origin
0208   Lett = 'ACGU';  
0209   Title =[Lett(modcode) ' shown at the origin, N1/N9 atom of second base shown by dots'];
0210 else
0211   Title = ['N1/N9 atom of second base shown by dots'];
0212 end
0213 
0214 % ------------------ Display class limits, if desired
0215 
0216 if isfield(ViewParam,'ClassLimits'),
0217  if ViewParam.ClassLimits == 1,                  % show all boxes
0218    B = CL(:,:,Paircodes(1));                     % use limits for this paircode
0219 
0220    for row = 1:length(B(:,1)),
0221      hold on
0222      if (abs(B(row,1)) < 14) & (abs(B(row,1)) > 0),
0223        zSquare([B(row,[2 4]) 0],[B(row,[3 5]) 0],'k');
0224        tex = zEdgeText(B(row,1),1,Paircodes(1));
0225        tex = tex(1:4);
0226        text(B(row,2)+0.1,B(row,4)+0.2,0,tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0227      end
0228    end
0229  elseif ViewParam.ClassLimits == 2,              % just show one box
0230    ClassLimits = zClassLimits;                   % load class limits
0231    B = CL(:,:,Paircodes(1));       % use limits for this paircode
0232 
0233    for row = 1:length(B(:,1)),
0234      hold on
0235      if any(fix(B(row,1)) == fix(Param.Category)),
0236        zSquare([B(row,[2 4]) 0],[B(row,[3 5]) 0],'k');
0237        tex = zEdgeText(B(row,1),1,Paircodes(1));
0238        tex = tex(1:4);
0239        text(B(row,2)+0.1,B(row,4)+0.2,0,tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0240      end
0241    end
0242  end
0243 end
0244 
0245 zoom on
0246 
0247 caxis(ColorAxis);
0248 view(2)
0249 axis equal
0250 
0251 if ViewParam.Normal == 1,
0252  Title = [Title ', lines are the normal vector of second base'];
0253 end
0254 title(Title);
0255 xlabel('Perpendicular to glycosidic bond');
0256 ylabel('Parallel to glycosidic bond');
0257 zlabel('Vertical with respect to first base');
0258 
0259 %----------------------------------------------------- Plot angle and normal
0260 figure(ViewParam.FigNum + 1)
0261 clf
0262 figure(ViewParam.FigNum + 1)
0263 clf
0264 
0265 %set(gcf,'Renderer','OpenGL');
0266 %set(gcf,'Renderer','zbuffer')
0267 
0268 cut = 90;                                             % cut at angle -cut
0269 
0270 for k = 1:length(Pair),                               % Loop through pairs
0271   p = Pair(k);
0272   e = p.Displ;
0273 
0274   c = Color(k);
0275 
0276   scatter3(mod(p.Ang+cut,360)-cut,p.Normal(3),p.Gap,18,c,'filled')
0277   hold on
0278 end
0279 
0280 xlabel('Angle of rotation (in degrees)');
0281 ylabel('Vertical component of normal');
0282 zlabel('Gap');
0283 title('');
0284 caxis(ColorAxis);
0285 grid on
0286 
0287 view(2)
0288 
0289 if isfield(ViewParam,'ClassLimits'),
0290   if ViewParam.ClassLimits == 1,
0291     B = CL(:,:,Paircodes(1));     % use limits for this paircode
0292 
0293     B(:,10:11) = mod(B(:,10:11)+cut,360)-cut;
0294     B(:,8:9)   = B(:,8:9) + 0.1*(rand(size(B(:,8:9)))-0.5).*(abs(B(:,8:9))>1);
0295 
0296     for row = 1:length(B(:,1)),
0297       hold on
0298       if (abs(B(row,1)) < 14) & (abs(B(row,1)) > 0),
0299         tex = zEdgeText(B(row,1),1,Paircodes(1));
0300         tex = tex(1:4);
0301         if (B(row,10) < B(row,11)),
0302           zSquare([B(row,[10 8]) 0],[B(row,[11 9]) 0],'k');
0303           text(B(row,10)+4,B(row,8)+0.08,B(row,11),tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0304         else
0305           zSquare([B(row,[10 8]) 0],[270 B(row,9) 0],'k');
0306           text(B(row,10)+4,B(row,8)+0.08,B(row,11),tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0307           zSquare([-90 B(row,8) 0],[B(row,[11 9]) 0],'k');
0308           text(-86,B(row,8)+0.08,B(row,11),tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0309         end
0310       end
0311     end
0312   elseif ViewParam.ClassLimits == 2,            % show only one box
0313     B = CL(:,:,Paircodes(1));     % use limits for this paircode
0314     B(:,10:11) = mod(B(:,10:11)+cut,360)-cut;
0315     B(:,8:9)   = B(:,8:9) + 0.1*(rand(size(B(:,8:9)))-0.5).*(abs(B(:,8:9))>1);
0316 
0317     for row = 1:length(B(:,1)),
0318       hold on
0319       if any(fix(B(row,1)) == fix(Param.Category)),
0320         tex = zEdgeText(B(row,1),1,Paircodes(1));
0321         tex = tex(1:4);
0322         if (B(row,10) < B(row,11)),
0323           zSquare([B(row,[10 8]) 0],[B(row,[11 9]) 0],'k');
0324           text(B(row,10)+4,B(row,8)+0.08,B(row,11),tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0325         else
0326           zSquare([B(row,[10 8]) 0],[270 B(row,9) 0],'k');
0327           text(B(row,10)+4,B(row,8)+0.08,B(row,11),tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0328           zSquare([-90 B(row,8) 0],[B(row,[11 9]) 0],'k');
0329           text(-86,B(row,8)+0.08,B(row,11),tex,'horizontalalignment','left','fontweight','bold','FontSize',8);
0330         end
0331       end
0332     end
0333     v = axis;
0334     axis([-cut 360-cut -1.1 1.1 -2 2]);
0335   end
0336 end
0337 
0338 zoom on
0339 
0340 %----------------------------------------------------- Plot angle and normal
0341 
0342 if ViewParam.Color > 1,
0343   figure(ViewParam.FigNum + 2)
0344   clf
0345   figure(ViewParam.FigNum + 2)
0346   clf
0347 
0348   hist(Color,30)
0349   switch ViewParam.Color,
0350     case {2, 3, 4}, title('Histogram of interaction code');
0351     case  5, title('Histogram of paircode'); % Color(k) = p.Paircode;
0352     case  6, title('Histogram of angle of rotation'); % Color(k) = p.Ang;
0353     case  7, title('Histogram of gap'); % Color(k) = p.Gap;
0354     case  8, title('Histogram of third component of normal vector');
0355     case  9, title('Histogram of distance'); % Color(k) = p.Distances(1);
0356     case 11, title('Histogram of stacking overlap'); % Color(k) = p.StackingOverlap;
0357     case 12, title('Histogram of interaction code'); % Color(k) = p.Edge;
0358     case 13, title('Histogram of perpendicual displacement'); % Color(k) = p.Displ(1);
0359     case 14, title('Histogram of parallel displacement'); % Color(k) = p.Displ(2);
0360     case 15, title('Histogram of vertical displacement'); % Color(k) = p.Displ(3);
0361   end
0362 end
0363 
0364 
0365 end
0366 % ---------------------------------------------------- Display the menu
0367 
0368   k=menu('Scatterplot controls',...
0369          'Increment first nucleotide',...             % 1
0370          'Decrement first nucleotide',...             % 2
0371          'Increment second nucleotide',...            % 3
0372          'Decrement second nucleotide',...            % 4
0373          'Color by angle',...                         % 5
0374          'Color by pair code',...                     % 6
0375          'Color by interaction',...                   % 7
0376          'Color by subcategory',...                   % 8
0377          'Color by displacement perp to bond',...     % 9
0378          'Color by displacement parallel to bond',... % 10
0379          'Color by vertical displacement',...         % 11
0380          'Color by position in list',...              % 12
0381          'Color by gap',...                           % 13
0382          'Color by third component of normal',...     % 14
0383          'Toggle normal vector',...                   % 15
0384          'Toggle category limits',...                 % 16
0385          'List pair parameters',...                   % 17
0386          'Quit');                                     % 18
0387 
0388   switch k,
0389   case 1, N1 = N1 + 1;
0390           if N1 > N,
0391             N1 = 1;
0392           end
0393           if N1 == N2,
0394             N1 = N2 + 1;
0395           end
0396           if N1 > N,
0397             N1 = 1;
0398           end
0399   case 2, N1 = N1 - 1;
0400           if N1 < 1,
0401             N1 = N;
0402           end
0403           if N1 == N2,
0404             N1 = N2 - 1;
0405           end
0406           if N1 < 1,
0407             N1 = N;
0408           end
0409   case 3, N2 = N2 + 1;
0410           if N2 > N,
0411             N2 = 1;
0412           end
0413           if N2 == N1,
0414             N2 = N2 + 1;
0415           end
0416           if N2 > N,
0417             N2 = 1;
0418           end
0419   case 4, N2 = N2 - 1;
0420           if N2 < 1,
0421             N2 = N;
0422           end
0423           if N2 == N1,
0424             N2 = N2 - 1;
0425           end
0426           if N2 < 1,
0427             N2 = N;
0428           end
0429   case  5, ViewParam.Color = 6;
0430   case  6, ViewParam.Color = 5;
0431   case  7, ViewParam.Color = 2;
0432   case  8, ViewParam.Color = 12;
0433   case  9, ViewParam.Color = 13;
0434   case 10, ViewParam.Color = 14;
0435   case 11, ViewParam.Color = 15;
0436   case 12, ViewParam.Color = 1;
0437   case 13, ViewParam.Color = 7;
0438   case 14, ViewParam.Color = 8;
0439   case 15, ViewParam.Normal = 1 - ViewParam.Normal;
0440            Recolor = 0;
0441   case 16, ViewParam.ClassLimits = 1 - ViewParam.ClassLimits;
0442            Recolor = 0;
0443   case 17, xListPairs(Pair,ViewParam);
0444            Replot = 0;
0445   case 18, Stop = 1;
0446   end
0447 
0448   if (k >= 1) && (k <= 4) && (N == 2),     % switch order
0449     N3 = N1;
0450     N1 = N2;
0451     N2 = N3;
0452   end
0453 
0454   if k < 17,
0455     Replot  = 1;
0456   end
0457 
0458   if k < 13,
0459     Recolor  = 1;
0460   end
0461 
0462   if k < 5,
0463     Reclassify = 1;
0464   end
0465 
0466 end
0467 
0468 FigsDone = 2;
0469 
0470 warning on