zInteractionRange

 zInteractionRange(File) returns a sparse matrix S which indicates
 how far a pairwise interaction is from being a local interaction

0001 % zInteractionRange(File) returns a sparse matrix S which indicates
0002 % how far a pairwise interaction is from being a local interaction
0003 
0004 function [File] = zInteractionRange(File, Verbose)
0005 
0006 if nargin < 2,
0007   Verbose = 1;
0008 end
0009 
0010 for f = 1:length(File),
0011 
0012 File(f).Crossing = sparse(zeros(length(File.NT)));
0013 File(f).Range    = sparse(zeros(length(File.NT)));
0014 
0015 if length(File(f).NT) > 1,
0016 
0017 E = fix(abs(File(f).Edge)); % matrix of pairwise interactions w/o BPh
0018 
0019 C = fix(abs(File(f).Edge))+0.001*(fix(abs(File(f).BasePhosphate+File(f).BasePhosphate')));                  % matrix of interactions, including base-phosphate
0020 
0021 for i = 1:length(C(:,1)),
0022   C(i,i) = 0;               % remove base-phosphate self interactions
0023 end
0024 
0025 C = triu(C);                % sparse matrix of pairwise interactions
0026 
0027 [i,j,c] = find(C);          % locations of all pairwise interactions
0028 
0029 c = fix(c);                 % remove decimal from base-phosphate
0030 
0031 k = find((j > i));          % ignore self interactions
0032 i = i(k);
0033 j = j(k);
0034 c = c(k);
0035 
0036 [y,k] = sort(j-i);          % order by distance from diagonal
0037 i = i(k);
0038 j = j(k);
0039 c = c(k);
0040 
0041 d = zeros(size(c));                          % current interaction range
0042 cww = zeros(size(c));                        % number of cWW's crossed
0043 
0044 for m = 1:length(i),                         % go through all pairs once
0045   if (c(m) == 1) && (d(m) == 0),             % cWW pair, not pseudoknotted
0046 
0047     p =     (i <= i(m)) .* (j >= i(m)) .* (j <= j(m));
0048     p = p + (i >= i(m)) .* (i <= j(m)) .* (j >= j(m));
0049     q = find(p);                      % indices of pairs that overlap i(m),j(m)
0050 
0051 
0052 
0053 %    d(q) = max(d(q), abs(i(m)-i(q))+abs(j(m)-j(q)));
0054                             % increase interaction range for these pairs
0055                             % to account for how far they cross i(m),j(m)
0056 
0057     d(q) = max([d(q) min([abs(i(m)-i(q)) abs(j(m)-j(q))],[],2)],[],2); 
0058                             % for all nested cWW pairs crossed, keep track
0059                             % of how many nucleotides this pair would need to
0060                             % move over to get uncrossed, keep the maximum
0061 
0062     p =     (i < i(m)) .* (j > i(m)) .* (j < j(m));
0063     p = p + (i > i(m)) .* (i < j(m)) .* (j > j(m));
0064     q = find(p);                      % indices of pairs that overlap i(m),j(m)
0065 
0066 
0067     cww(q) = cww(q) + 1;        % count number of nested cWW's crossed
0068 
0069     d(m) = 0;               % don't inadvertently make this non-nested!
0070 
0071     if Verbose > 2,
0072       figure(2)
0073       clf
0074       plot(j(q),i(q),'.');
0075       hold on
0076       plot(j(m),i(m),'r.');
0077       axis([1 max(i) 1 max(i)]);
0078       axis ij
0079       drawnow
0080       if Verbose > 3,
0081         pause
0082       end
0083     end
0084 
0085   end
0086 end
0087 
0088 d = min(d,abs(i-j)-1);              % adjust for local in-strand interactions
0089 
0090 i = [i; length(File(f).NT)];        % append a blank interaction for size
0091 j = [j; length(File(f).NT)];
0092 d = [d; 0];
0093 c = [c; 0];
0094 cww = [cww; 0];
0095 
0096 S = sparse(i,j,d);
0097 S = S + S';
0098 File(f).Range = S;
0099 
0100 S = sparse(i,j,cww);
0101 S = S + S';
0102 File(f).Crossing = S;
0103 
0104 if Verbose > 1,
0105   fprintf('%s has %d basepairs, of which %d are local.\n', File(f).Filename, full(sum(sum((C > 0) .* (C < 15)))), full(sum(sum(((S<=10).*C > 0) .* ((S<=10).*C < 15)))));
0106   for m = 1:14,
0107     all = length(find(c==m));
0108     local = length(find((c==m).*(d<=10)));
0109     longrange = all - local;
0110     fprintf('%4s %3d total, %3d local, %3d long-range\n', zEdgeText(m), all, local, longrange);
0111   end
0112 
0113   fprintf('Long-range cWW interactions:\n');
0114   q = find((c==1).*(d>10));
0115   [y,h] = sort(i(q));
0116   q = q(h);
0117   for m = 1:length(q),
0118     fprintf('%s%s with %s%s\n', File(f).NT(i(q(m))).Base,File(f).NT(i(q(m))).Number, File(f).NT(j(q(m))).Base, File(f).NT(j(q(m))).Number);
0119   end
0120 end
0121 
0122 else
0123 
0124   File(f).Range = [];
0125 
0126 end
0127 
0128 end