xPairwiseScreen

 xPairwiseScreen returns a sparse matrix with non-zero entries corresponding to pairs of bases which satisfy all given constraints
 Codes are the base codes from File

0001 % xPairwiseScreen returns a sparse matrix with non-zero entries corresponding to pairs of bases which satisfy all given constraints
0002 % Codes are the base codes from File
0003 %
0004 
0005 function [Screen] = xPairwiseScreen(File,Codes,Query,p,q,PC);
0006 
0007 if Query.Geometric == 0,
0008   D = File.Distance .* (File.Distance < Query.Diameter);
0009                                         % cap distance for non-geometric search
0010 else
0011   D = File.Distance;
0012 end
0013 
0014 % --------- Screen according to interaction between nucleotides
0015 
0016 if isfield(Query,'EdgeNums'),           % if screening by edges, incl stacking
0017   if length(Query.EdgeNums{p,q} > 0),
0018     E = sparse(zeros(size(D)));
0019     for i=1:length(Query.EdgeNums{p,q}),
0020       E = E + (fix(File.Edge) == Query.EdgeNums{p,q}(i));
0021     end
0022     D = D .* (E > 0);                         % include only those that match
0023   end
0024 end
0025     
0026 if isfield(Query,'ExcludeEdges'),                 % if excluding by edges
0027   if length(Query.ExcludeEdges{p,q} > 0),
0028     E = sparse(zeros(size(D)));
0029     for i=1:length(Query.ExcludeEdges{p,q}),
0030       E = E + (fix(File.Edge) == Query.ExcludeEdges{p,q}(i));
0031     end
0032     D = D .* (E == 0);
0033   end
0034 end
0035     
0036 if isfield(Query,'OKPairs'),                 % if screening by paircode
0037   if length(Query.OKPairs{p,q} > 0),
0038     E = sparse(zeros(size(D)));
0039     for i=1:length(Query.OKPairs{p,q}),
0040       E = E + sparse(PC == Query.OKPairs{p,q}(i));
0041     end
0042     D = D .* (E > 0);
0043   end
0044 end
0045     
0046 if isfield(Query,'ExPairs'),                 % if screening by paircode
0047   if length(Query.ExPairs{p,q} > 0),
0048     E = sparse(zeros(size(D)));
0049     for i=1:length(Query.ExPairs{p,q}),
0050       E = E + sparse(PC == Query.ExPairs{p,q}(i));
0051     end
0052     D = D .* (E == 0);
0053   end
0054 end
0055 
0056 if isfield(Query,'OKBB'),           % if screening by backbone
0057   if length(Query.OKBB{p,q} > 0),
0058     E = sparse(zeros(size(D)));
0059     for i=1:length(Query.OKBB{p,q}),
0060       E = E + (fix(File.Backbone) + fix(File.Backbone') == Query.OKBB{p,q}(i));
0061     end
0062     D = D .* (E > 0);                         % include only those that match
0063   end
0064 end
0065     
0066 if isfield(Query,'ExBB'),                 % if excluding by backbone
0067   if length(Query.ExBB{p,q} > 0),
0068     E = sparse(zeros(size(D)));
0069     for i=1:length(Query.ExBB{p,q}),
0070       E = E + (fix(File.Backbone) + fix(File.Backbone') == Query.ExBB{p,q}(i));
0071     end
0072     D = D .* (E == 0);
0073   end
0074 end
0075 
0076 if isfield(Query,'BasePhos'),                 % if screening by base-phosphate
0077   if length(Query.BasePhos{p,q} > 0),
0078     E = sparse(zeros(size(D)));
0079     for i=1:length(Query.BasePhos{p,q}),
0080       E = E + (fix(File.BasePhosphate) == Query.BasePhos{p,q}(i));
0081     end
0082     D = D .* (E > 0);
0083   end
0084   if length(Query.BasePhos{q,p} > 0),
0085     E = sparse(zeros(size(D)));
0086     for i=1:length(Query.BasePhos{q,p}),
0087       E = E + (fix(File.BasePhosphate') == Query.BasePhos{q,p}(i));
0088     end
0089     D = D .* (E > 0);
0090   end
0091 end
0092     
0093 if isfield(Query,'ExcludeBasePhos'),                 % if excluding by edges
0094   if length(Query.ExcludeBasePhos{p,q} > 0),
0095     E = sparse(zeros(size(D)));
0096     for i=1:length(Query.ExcludeBasePhos{p,q}),
0097       E = E + (fix(File.BasePhosphate) == Query.ExcludeBasePhos{p,q}(i));
0098     end
0099     D = D .* (E == 0);
0100   end
0101   if length(Query.ExcludeBasePhos{q,p} > 0),
0102     E = sparse(zeros(size(D)));
0103     for i=1:length(Query.ExcludeBasePhos{q,p}),
0104       E = E + (fix(File.BasePhosphate') == Query.ExcludeBasePhos{q,p}(i));
0105     end
0106     D = D .* (E == 0);
0107   end
0108 end
0109 
0110 if isfield(Query,'Flank'),
0111   if ~isempty(Query.Flank{p,q}),
0112     D = D .* File.Flank;                 % only keep pairs between nested cWW
0113   end
0114 end
0115 
0116 if isfield(Query,'Range'),               % screen by range restriction
0117   if ~isempty(Query.Range{p,q}),
0118     R = Query.Range{p,q};
0119     D = D .* (File.Range >= R(1)) .* (File.Range <= R(2));
0120   end
0121 end
0122 
0123 if isfield(Query,'Coplanar'),            % screen by coplanarity
0124  if ~isempty(Query.Coplanar{p,q}),
0125 
0126   switch Query.Coplanar{p,q}
0127   case 0,
0128     D = D .* (File.Coplanar < 0.5);      % ~cp
0129   case 1,
0130     D = D .* (File.Coplanar >= 0.5);     % cp
0131   case 2,
0132     D = D .* (File.Coplanar > 0);        % ncp cp
0133   case 3,
0134     D = D .* (File.Coplanar > 0) .* (File.Coplanar < 0.5);  % ncp
0135   case 4,
0136     D = D .* ((File.Coplanar == 0) + (File.Coplanar >= 0.5)); % ~ncp
0137   case 5,
0138     D = D .* (File.Coplanar == 0);       % ~ncp ~cp
0139   end
0140 
0141 %figure(3)
0142 %hist(nonzeros(D),30)
0143 
0144  end
0145 end
0146 
0147     
0148 [i,j] = find(D);         % nucleotide pairs with OK distances and interactions
0149 d = nonzeros(D);         % distances below the large cutoff
0150 
0151 % --------- Screen according to maximum difference in nucleotide numbers
0152 
0153 if isfield(Query,'MaxDiffMat'),
0154   if Query.MaxDiffMat(p,q) < Inf,
0155     k = find(abs(i-j) <= Query.MaxDiffMat(p,q)); %retain pairs close enough together
0156     i = i(k);
0157     j = j(k);
0158     d = d(k);
0159   end
0160 end
0161 
0162 % --------- Screen according to minimum difference in nucleotide numbers
0163 
0164 if isfield(Query,'MinDiffMat'),
0165   if Query.MinDiffMat(p,q) > 1,
0166     k = find(abs(i-j) >= Query.MinDiffMat(p,q)); %retain pairs far enough apart
0167     i = i(k);
0168     j = j(k);
0169     d = d(k);
0170   end
0171 end
0172 
0173 % --------- Screen according to sign of nucleotide number difference
0174 
0175 if isfield(Query,'DifferenceSignMat'),
0176   if Query.DifferenceSignMat(p,q) < 0,
0177     k = find(j > i);                          % retain pairs with num2>num1
0178     i = i(k);
0179     j = j(k);
0180     d = d(k);
0181   elseif Query.DifferenceSignMat(p,q) > 0,
0182     k = find(j < i);                          % retain pairs with num2<num1
0183     i = i(k);
0184     j = j(k);
0185     d = d(k);
0186   end
0187 end
0188 
0189 % --------- Screen according to the nucleotide mask
0190 
0191 if (min(Query.OKCodes{p}) == 0) | (min(Query.OKCodes{q}) == 0),
0192   if (min(Query.OKCodes{p}) == 1) & (min(Query.OKCodes{q}) == 0),
0193     k = find(Query.OKCodes{q}(Codes(j)));
0194   elseif (min(Query.OKCodes{p}) == 0) & (min(Query.OKCodes{q}) == 1),
0195     k = find(Query.OKCodes{p}(Codes(i)));
0196   else
0197     k = find(Query.OKCodes{p}(Codes(i)) .* Query.OKCodes{q}(Codes(j)));
0198   end
0199 
0200   i = i(k); 
0201   j = j(k);
0202   d = d(k);
0203 end
0204 
0205 % --------- Screen according to configuration (syn or anti)
0206 
0207 if length(Query.Config{p}) > 0,
0208   switch Query.Config{p}
0209     case 'syn'
0210       k = find(cat(1,File.NT(i).Syn) == 1);
0211     case 'anti'
0212       k = find(cat(1,File.NT(i).Syn) == 0);
0213     otherwise
0214       k = 1:length(i);
0215   end
0216 
0217   i = i(k);
0218   j = j(k);
0219   d = d(k);
0220 end
0221 
0222 if length(Query.Config{q}) > 0,
0223   switch Query.Config{q}
0224     case 'syn'
0225       k = find(cat(1,File.NT(j).Syn) == 1);
0226     case 'anti'
0227       k = find(cat(1,File.NT(j).Syn) == 0);
0228     otherwise
0229       k = 1:length(j);
0230   end
0231 
0232   i = i(k);
0233   j = j(k);
0234   d = d(k);
0235 end
0236 
0237 % --------- Screen according to pairwise distance in model
0238 
0239 if (Query.Geometric > 0),
0240 
0241   % --------- Compute square of distance difference from model
0242 
0243   d = (d - Query.Distance(p,q)).^2;   % squared difference in distances
0244 
0245   d = d + 0.00000001 * (d == 0);       % avoid rejecting model; make d nonzero
0246 
0247   % --------- Impose upper limit on distance differences; 2-nucleotide cutoff
0248 
0249 
0250   if Query.NumNT > 2,
0251     Wp = Query.LocWeight(p);
0252     Wq = Query.LocWeight(q);
0253     MaxD = (Wp + Wq) * (Query.NumNT * Query.DiscCutoff)^2 / (Wp * Wq);
0254   else
0255     Wp = 1;
0256     Wq = 1;
0257     MaxD = (Query.NumNT * Query.DiscCutoff)^2;
0258   end
0259 
0260   if isfield(Query,'Flex'),
0261     MaxD = max(MaxD,Query.Flex(p,q)^2);  % allow larger distance if desired
0262   end
0263 
0264   k = find(d <= MaxD);            % keep ones with small difference from model
0265 
0266   i = i(k);
0267   j = j(k);
0268   d = d(k) * Wp * Wq;
0269 
0270 end
0271 
0272 % --------- Construct sparse matrix of retained distance difference squares
0273 
0274 % [length(i) length(j) length(d) max(k)];
0275 
0276 Screen = sparse(i,j,d,File.NumNT,File.NumNT,length(i));
0277 
0278 
0279