001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021package org.biojava.nbio.structure.domain.pdp; 022 023import org.biojava.nbio.structure.Atom; 024 025import java.util.List; 026 027public class Cut { 028 029 static boolean verbose = CutDomain.verbose; 030 031 public int cut( Atom[] ca, Domain dom, CutValues val, int[][] dist, PDPDistanceMatrix pdpMatrix) { 032 033 int nclose = pdpMatrix.getNclose(); 034 035 int[] iclose = pdpMatrix.getIclose(); 036 int[] jclose = pdpMatrix.getJclose(); 037 038 int[] contacts = new int[PDPParameters.MAXLEN]; 039 double[] max_contacts = new double [PDPParameters.MAXLEN]; 040 double[] contact_density = new double[PDPParameters.MAXLEN]; 041 double average_density,x,y; 042 043 int endsf,endst; 044 int k,l,nc; 045 int size1,size2; 046 int size1t,size2t; 047 int size11,size22,size0; 048 int contactsd; 049 int iseg,jseg,kseg; 050 int from,to,from1,to1,from2,to2,lseg; 051 052 053 int site_min = -1; 054 055 056 // AP add sort here.. 057 //qsort(dom.segment,dom.nseg,sizeof(struct Segment),segcmp); 058 // what is going on with the segments?? 059 060 List<Segment> segments = dom.getSegments(); 061 062 java.util.Collections.sort(segments, new SegmentComparator()); 063 064 if ( verbose) 065 System.out.println(" --- Cut.cut " + dom + " "); 066 average_density = 0.0; 067 size0=0; 068 for(iseg=0;iseg<dom.nseg;iseg++) { 069 contactsd=1; 070 size1t=0; 071 size2t=0; 072 for(jseg=0;jseg<iseg;jseg++) 073 size1t+=(dom.getSegmentAtPos(jseg).getFrom() - dom.getSegmentAtPos(jseg).getFrom() + 1); 074 for(jseg=iseg+1;jseg<dom.nseg;jseg++) 075 size2t+=(dom.getSegmentAtPos(jseg).getTo() - dom.getSegmentAtPos(jseg).getFrom() + 1); 076 for(jseg=0;jseg<iseg;jseg++) { 077 from1 = dom.getSegmentAtPos(jseg).getFrom(); 078 to1 = dom.getSegmentAtPos(jseg).getTo(); 079 for(int i=from1;i<to1;i++) { 080 for(kseg=iseg+1;kseg<dom.nseg;kseg++) { 081 from2 = dom.getSegmentAtPos(kseg).getFrom(); 082 to2 = dom.getSegmentAtPos(kseg).getFrom(); 083 for(int j=from2;j<to2;j++) 084 if(Math.abs(i-j)>4) contactsd+=(dist[i][j]); 085 } 086 } 087 } 088 from = dom.getSegmentAtPos(iseg).getFrom(); 089 to = dom.getSegmentAtPos(iseg).getTo(); 090 for(k=from;k<to;k++) { 091 contacts[k] = contactsd; 092 /* 093 if(k==392) printf("init contacts = %d\n",contacts[k]); 094 */ 095 size11=size1t+(k-from+1); 096 size22=size2t+(to-k); 097 for(int i=from;i<=k;i++) { 098 for(kseg=iseg+1;kseg<dom.nseg;kseg++) { 099 from2 = dom.getSegmentAtPos(kseg).getFrom(); 100 to2 = dom.getSegmentAtPos(kseg).getTo(); 101 for(int j=from2;j<=to2;j++) 102 if(Math.abs(i-j)>4) contacts[k]+=(dist[i][j]); 103 } 104 } 105 /* 106 if(k==392) printf("[from,k]x]iseg,nseg[ = %d\n",contacts[k]); 107 */ 108 for(int i=from;i<=k;i++) { 109 for (int j=k+1;j<=to;j++) 110 if(Math.abs(i-j)>4) contacts[k]+=(dist[i][j]); 111 } 112 /* 113 if(k==392) printf("[from,k]x]k,to[ = %d\n",contacts[k]); 114 */ 115 for(int i=k+1;i<=to;i++) { 116 for(kseg=0;kseg<iseg;kseg++) { 117 from2 = dom.getSegmentAtPos(kseg).getFrom(); 118 to2 = dom.getSegmentAtPos(kseg).getTo(); 119 for(int j=from2;j<to2;j++) 120 if(Math.abs(i-j)>4) contacts[k]+=(dist[j][i]); 121 } 122 } 123 /* 124 if(k==392) printf("]k,to]x]0,iseg[ = %d\n",contacts[k]); 125 */ 126 size1=Math.min(size11,size22); 127 size2=Math.max(size11,size22); 128 x=Math.min(PDPParameters.MAXSIZE,size1); 129 y=Math.min(PDPParameters.MAXSIZE,size2); 130 if(x>150&&y>1.5*x) y=1.5*x; 131 else if(y>2*x) y=2*x; 132 /* 133 */ 134 x=Math.min(Math.pow(x,1.3/3)+PDPParameters.RG,Math.pow(x,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG); 135 y=Math.min(Math.pow(y,1.3/3)+PDPParameters.RG,Math.pow(y,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG); 136 /* max_ contacts depend on the size of domains */ 137 /* stella wanted comments at this point */ 138 /* 139 max_contacts[k] = min(MAXCONT,x*y); 140 */ 141 max_contacts[k] = 10*x*y; 142 if(size1>150) max_contacts[k] = 9*x*y; 143 contact_density[k]=contacts[k]/max_contacts[k]; 144 /* 145 if(contact_density[k]>2.5) contact_density[k]=2.5; 146 */ 147 /* 148 if(first_cut) 149 if(k==277 && !first_cut) 150 printf("k=%d s1=%d s2=%d x=%f y=%f mc=%d c=%d cd=%f\n",k,size1,size2,x,y,max_contacts[k],contacts[k],contact_density[k]); 151 */ 152 //if(verbose) System.out.println(String.format("%d %d %d %f %f %f %d %f",k,size1,size2,x,y,max_contacts[k],contacts[k],contact_density[k])); 153 /* 154 */ 155 if(from==0) endsf = PDPParameters.ENDSEND; 156 else endsf = PDPParameters.ENDS; 157 if(to==ca.length-1) endst = PDPParameters.ENDSEND; 158 else endst = PDPParameters.ENDS; 159 if((contact_density[k])<val.s_min&&k>from+endsf&&k<to-endst) { 160 val.s_min = (contact_density[k]); 161 site_min=k+1; 162 } 163 if(k>from+endsf&&k<to-endst) { 164 average_density+=contact_density[k]; 165 size0++; 166 } 167 } 168 } 169 average_density/=size0; 170 if(verbose) System.out.printf(" --- Trying to cut domain of size %d having %d segments and average cont_density %f%n%n",dom.size,dom.nseg,average_density); 171 172 if ( verbose ) 173 for(kseg=0;kseg<dom.nseg;kseg++) 174 System.out.printf(" --- segment %d from %d to %d av density: %f%n",kseg,dom.getSegmentAtPos(kseg).getFrom(),dom.getSegmentAtPos(kseg).getTo(), average_density); 175 176 177 if(val.first_cut) { 178 /* 179 for(k=from+ENDS;k<to-ENDS;k++) 180 printf("%d %s %d %d %d %d %f\n",k,protein.res[k].type,k-from,to-k,contacts[k],max_contacts[k],contact_density[k]); 181 */ 182 val.AD = average_density; 183 if(verbose) System.out.printf(" --- AD=%f%n", average_density); 184 /* 185 */ 186 } 187 val.AD = average_density; 188 189 val.s_min/=val.AD; 190 191 if(verbose) System.out.printf(" --- after single cut: s_min = %f site_min = %d%n",val.s_min,site_min); 192 193 /// 194 k=0; 195 196 /* check double cuts */ 197 if ( verbose ) 198 System.out.println(" --- checking double cuts up to: " + nclose); 199 nc=0; 200 for(l=0;l<nclose;l++) { 201 202 /************ find iseg, jseg ****************/ 203 iseg=jseg=-1; 204 for(kseg=0;kseg<dom.nseg;kseg++) { 205 from=dom.getSegmentAtPos(kseg).getFrom(); 206 to=dom.getSegmentAtPos(kseg).getTo(); 207 if(from==0) endsf = PDPParameters.ENDSEND; 208 else endsf = PDPParameters.ENDS; 209 if(to==ca.length-1) endst = PDPParameters.ENDSEND; 210 else endst = PDPParameters.ENDS; 211 if(iclose[l]>from+endsf&&iclose[l]<to-endst) 212 iseg=kseg; 213 if(jclose[l]>from+endsf&&jclose[l]<to-endst) 214 jseg=kseg; 215 } 216 if(iseg<0||jseg<0) continue; 217 /* 218 printf("l = %d iclose[l] = %d jclose[l] = %d\n",l,iclose[l],jclose[l]); 219 */ 220 /*********************************************/ 221 222 from=dom.getSegmentAtPos(iseg).getFrom(); 223 to=dom.getSegmentAtPos(iseg).getTo(); 224 from1=dom.getSegmentAtPos(jseg).getFrom(); 225 to1=dom.getSegmentAtPos(jseg).getTo(); 226 227 /************ count contacts *****************/ 228 contacts[nc] = 1; 229 //no=0; 230 231 /******* contacts between [0,iseg[ and ]iseg,jseg[ ********/ 232 for(kseg=0;kseg<iseg;kseg++) 233 for(lseg=iseg+1;lseg<jseg;lseg++) 234 for( int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++) 235 for(int j=dom.getSegmentAtPos(lseg).getFrom();j<dom.getSegmentAtPos(lseg).getTo();j++) { 236 contacts[nc]+=(dist[i][j]); 237 } 238 239 //System.out.println(String.format("[0,iseg[ - ]iseg,jseg[ : %d\n",contacts[nc]-no)); 240 241 //no=contacts[nc]; 242 /**********************************************************/ 243 /**********************************************************/ 244 /******* contacts between ]jseg,nseg[ and ]iseg,jseg[ ********/ 245 for(kseg=jseg+1;kseg<dom.nseg;kseg++) 246 for(lseg=iseg+1;lseg<jseg;lseg++) 247 for(int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++) 248 for(int j=dom.getSegmentAtPos(lseg).getFrom();j<dom.getSegmentAtPos(lseg).getTo();j++) { 249 contacts[nc]+=(dist[j][i]); 250 } 251 /* 252 printf("]jseg,nseg] - ]iseg,jseg[ : %d\n",contacts[nc]-no); 253 */ 254 //no=contacts[nc]; 255 /*************************************************************/ 256 /*************************************************************/ 257 /**** contacts between [from,iclose] in iseg and ]iseg,jseg[ ****/ 258 if(iseg==jseg) { 259 //System.out.println(" CONTACT: " + from + " " + iclose[l] + " " + iseg + " " + jseg); 260 for(int i=from;i<=iclose[l];i++) { 261 for (int j=iclose[l]+1;j<=jclose[l];j++) { 262 contacts[nc]+=(dist[i][j]); 263 } 264 } 265 for (int j=iclose[l]+1;j<jclose[l];j++) { 266 for(kseg=0;kseg<iseg;kseg++) 267 for(int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++) { 268 contacts[nc]+=(dist[i][j]); 269 } 270 for(int i=jclose[l];i<to;i++) { 271 contacts[nc]+=(dist[j][i]); 272 } 273 for(kseg=iseg+1;kseg<dom.nseg;kseg++) 274 for(int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++) { 275 contacts[nc]+=(dist[j][i]); 276 } 277 } 278 /* 279 printf("iclose==jclose : %d\n",contacts[nc]-no); 280 */ 281 //no=contacts[nc]; 282 } 283 else { 284 //System.out.println(" ISEG!=JSEG " + " " + from + " " + iclose[l]); 285 for(int i=from;i<=iclose[l];i++) { 286 for(kseg=iseg+1;kseg<jseg;kseg++) 287 for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) { 288 contacts[nc]+=(dist[i][j]); 289 } 290 for(int j=from1;j<jclose[l];j++) { 291 contacts[nc]+=(dist[i][j]); 292 } 293 for(int j=iclose[l]+1;j<to;j++) { 294 contacts[nc]+=(dist[i][j]); 295 } 296 } 297 for(int i=iclose[l]+1;i<to;i++) { 298 for(kseg=0;kseg<iseg;kseg++) 299 for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) { 300 contacts[nc]+=(dist[j][i]); 301 } 302 for(kseg=jseg+1;kseg<dom.nseg;kseg++) 303 for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) { 304 contacts[nc]+=(dist[i][j]); 305 } 306 for(int j=jclose[l];j<=to1;j++) { 307 contacts[nc]+=(dist[i][j]); 308 } 309 } 310 for (int i=from1;i<jclose[l];i++) { 311 for(kseg=0;kseg<iseg;kseg++) 312 for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) { 313 contacts[nc]+=(dist[j][i]); 314 } 315 for(kseg=jseg+1;kseg<dom.nseg;kseg++) { 316 for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) 317 contacts[nc]+=(dist[i][j]); 318 } 319 for(int j=jclose[l];j<to1;j++) { 320 contacts[nc]+=(dist[i][j]); 321 } 322 } 323 for(int i=jclose[l];i<to1;i++) 324 for(kseg=iseg+1;kseg<jseg;kseg++) 325 for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) { 326 /* 327 if(iclose[l]==33&&jclose[l]==69&&dist[i][j]) printf("%d %s %d %s %d\n",i,protein.res[i].type,j,protein.res[j].type,dist[i][j]); 328 */ 329 contacts[nc]+=(dist[j][i]); 330 } 331 } 332 /*******************************************************************/ 333 /*******************************************************************/ 334 /*******************************************************************/ 335 size11=0; 336 size22=0; 337 for(kseg=0;kseg<iseg;kseg++) 338 size11+=(dom.getSegmentAtPos(kseg).getTo()-dom.getSegmentAtPos(kseg).getFrom()+1); 339 for(kseg=jseg+1;kseg<dom.nseg;kseg++) 340 size11+=(dom.getSegmentAtPos(kseg).getTo()-dom.getSegmentAtPos(kseg).getFrom()+1); 341 size11+=(iclose[l]-from+1); 342 size11+=(to1-jclose[l]+1); 343 /* 344 printf("size11 = %d from = %d to1 = %d \n",size11,from,to1); 345 */ 346 for(kseg=iseg+1;kseg<jseg;kseg++) 347 size22+=(dom.getSegmentAtPos(kseg).getTo()-dom.getSegmentAtPos(kseg).getFrom()+1); 348 if(iseg==jseg) 349 size22+=(jclose[l]-iclose[l]); 350 else { 351 size22+=(jclose[l]-from1); 352 size22+=(to-iclose[l]); 353 } 354 /* 355 printf("size22 = %d from1 = %d jclose %d iclose %d \n",size22,from1,jclose[l],iclose[l]); 356 */ 357 size1=Math.min(size11,size22); 358 size2=Math.max(size11,size22); 359 x=Math.min(PDPParameters.MAXSIZE,size1); 360 y=Math.max(PDPParameters.MAXSIZE,size2); 361 if(y>2*x) y=2*x; 362 /* 363 */ 364 x=Math.min(Math.pow(x,1.3/3)+PDPParameters.RG,Math.pow(x,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG); 365 y=Math.min(Math.pow(y,1.3/3)+PDPParameters.RG,Math.pow(y,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG); 366 /* 367 max_contacts[nc] = min(MAXCONT,x*y); 368 */ 369 max_contacts[nc] = x*y*10; 370 if(size1>150) max_contacts[k] = 9*x*y; 371 contact_density[nc]=contacts[nc]/max_contacts[nc]; 372 /* 373 if(first_cut) 374 */ 375 //if ( verbose) 376 // System.out.println("double cut" + l + " " + ca[iclose[l]].getGroup().getType() + " " + iclose[l] + " " + jclose[l] + " " + contacts[nc]+ 377 // " " + max_contacts[nc] + " " + x+ " " + y + " " + size11 + " " + size22 + " " + contact_density[nc] + " " + contact_density[nc]/val.AD ); 378 //if(verbose) System.out.println(String.format(" double cut: %d %s %d %d c=%d mc=%d x=%f y=%f s1=%d s2=%d cd=%f cd/ad=%f\n",l,ca[iclose[l]].getGroup().getType(),iclose[l],jclose[l],contacts[nc],max_contacts[nc],x,y,size11,size22,contact_density[nc],contact_density[nc]/val.AD)); 379 if((contact_density[nc]/val.AD+PDPParameters.DBL)<val.s_min&&contact_density[nc]/val.AD+PDPParameters.DBL<PDPParameters.CUT_OFF_VALUE2) { 380 /* 381 printf("________________\n"); 382 */ 383 val.s_min = (contact_density[nc]/val.AD)+PDPParameters.DBL; 384 site_min=iclose[l]; 385 val.site2=jclose[l]; 386 } 387 388 nc++; 389 if ( nc >= PDPParameters.MAXSIZE) 390 nc = PDPParameters.MAXSIZE-1; 391 } 392 val.first_cut=false; 393 if(verbose) 394 System.out.printf(" --- E ... at the end of cut: s_min %f CUTOFF %f site_min %d *site2 %d%n",val.s_min,PDPParameters.CUT_OFF_VALUE,site_min,val.site2); 395 if(val.s_min> PDPParameters.CUT_OFF_VALUE) return -1; 396 397 return(site_min); 398 } 399}