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 * Created on Jan 28, 2006 021 * 022 */ 023package org.biojava.nbio.structure.align.pairwise; 024 025import org.biojava.nbio.structure.Atom; 026import org.biojava.nbio.structure.Calc; 027import org.biojava.nbio.structure.StructureException; 028import org.biojava.nbio.structure.align.StrucAligParameters; 029import org.biojava.nbio.structure.align.helper.AlignUtils; 030import org.biojava.nbio.structure.align.helper.JointFragments; 031import org.biojava.nbio.structure.jama.Matrix; 032import org.slf4j.Logger; 033import org.slf4j.LoggerFactory; 034 035import java.io.Serializable; 036import java.util.ArrayList; 037import java.util.Collections; 038import java.util.Comparator; 039import java.util.List; 040 041 042/** 043 * Joins the initial Fragments together to larger Fragments 044 * 045 * @author Andreas Prlic 046 * @author Peter Lackner 047 * @since 1.5 048 * @version %I% %G% 049 */ 050public class FragmentJoiner { 051 052 public static final Logger logger = LoggerFactory.getLogger(FragmentJoiner.class); 053 054 public FragmentJoiner() { 055 super(); 056 057 } 058 059 /** 060 * Reallocates an array with a new size, and copies the contents 061 * of the old array to the new array. 062 * @param oldArray the old array, to be reallocated. 063 * @param newSize the new array size. 064 * @return A new array with the same contents. 065 */ 066 @SuppressWarnings("rawtypes") 067 public static Object resizeArray (Object oldArray, int newSize) { 068 int oldSize = java.lang.reflect.Array.getLength(oldArray); 069 Class elementType = oldArray.getClass().getComponentType(); 070 Object newArray = java.lang.reflect.Array.newInstance( 071 elementType,newSize); 072 int preserveLength = Math.min(oldSize,newSize); 073 if (preserveLength > 0) 074 System.arraycopy (oldArray,0,newArray,0,preserveLength); 075 return newArray; 076 } 077 078 079 /** 080 * In helices often many similar fragments can be found. To reduce these to a few 081 * representative ones this check can be used. It does a distance check between 082 * all known Fragments and a new one. If this one is on a similar diagonal and it 083 * has a lower rms, this one is a better representation. Note: shifts of one are 084 * not allowed. 085 * 086 * @param fragments 087 * @param f 088 * @param rmsmat 089 * @return true - if this is a better representant for a group of locala fragments. 090 */ 091 public static boolean reduceFragments(List<FragmentPair> fragments, FragmentPair f, Matrix rmsmat){ 092 boolean doNotAdd =false; 093 int i = f.getPos1(); 094 int j = f.getPos2(); 095 096 for ( int p =0; p < fragments.size(); p++){ 097 FragmentPair tmp = fragments.get(p); 098 099 int di1 = Math.abs(f.getPos1() - tmp.getPos1()); 100 int di2 = Math.abs(f.getPos2() - tmp.getPos2()); 101 if (( Math.abs(di1-di2) == 2)) { 102 double rms1 = rmsmat.get(tmp.getPos1(),tmp.getPos2()); 103 double rms2 = rmsmat.get(i,j); 104 doNotAdd = true; 105 if ( rms2 < rms1){ 106 fragments.remove(p); 107 fragments.add(f); 108 break; 109 } 110 p = fragments.size(); 111 } 112 } 113 return doNotAdd; 114 } 115 116 117 public JointFragments[] approach_ap3(Atom[] ca1, Atom[] ca2, FragmentPair[] fraglst, 118 StrucAligParameters params) throws StructureException { 119 120 //the final list of joined fragments stores as apairs 121 List<JointFragments> fll = new ArrayList<JointFragments>(); 122 123 FragmentPair[] tmpfidx = new FragmentPair[fraglst.length]; 124 for ( int i=0 ; i < fraglst.length; i++){ 125 tmpfidx[i] = (FragmentPair)fraglst[i].clone(); 126 } 127 128 int n = tmpfidx.length; 129 130 // if two fragments after having been joint have rms < X 131 // keep the new fragment. 132 133 for (int i=0;i< fraglst.length;i++){ 134 135 boolean[] used = new boolean[n]; 136 137 int p1i = tmpfidx[i].getPos1(); 138 int p1j = tmpfidx[i].getPos2(); 139 int l1 = tmpfidx[i].getLength(); 140 141 JointFragments f = new JointFragments(); 142 143 int maxi=p1i+l1-1; 144 145 f.add(p1i,p1j,0,l1); 146 used[i] = true; 147 148 //n = tmpfidx.length; 149 150 for (int j=(i+1);j<n;j++){ 151 152 if ( used[j]) 153 continue; 154 155 int p2i = tmpfidx[j].getPos1(); 156 int p2j = tmpfidx[j].getPos2(); 157 int l2 = tmpfidx[j].getLength(); 158 159 if ( p2i > maxi) { 160 161 162 // TODO: replace this with plo angle calculation 163 if ( params.isDoAngleCheck()){ 164 165 // was 0.174 166 if (! angleCheckOk(tmpfidx[i], tmpfidx[j], 0.4f)) 167 continue; 168 } 169 170 if ( params.isDoDistanceCheck()) { 171 172 if (! distanceCheckOk(tmpfidx[i],tmpfidx[j], params.getFragmentMiniDistance())) 173 continue; 174 } 175 176 if ( params.isDoDensityCheck()) { 177 178 if ( ! densityCheckOk(ca1,ca2, f.getIdxlist(), p2i,p2j, l2, params.getDensityCutoff())) 179 continue; 180 } 181 182 183 if ( params.isDoRMSCheck()) { 184 185 double rms = rmsCheck(ca1,ca2, f.getIdxlist(), p2i, p2j, l2); 186 if ( rms > params.getJoinRMSCutoff()) 187 continue; 188 f.setRms(rms); 189 } 190 191 f.add(p2i,p2j,0,l2); 192 used[j] = true; 193 maxi = p2i+l2-1; 194 195 196 } 197 } 198 fll.add(f); 199 } 200 201 202 Comparator<JointFragments> comp = new JointFragmentsComparator(); 203 Collections.sort(fll,comp); 204 Collections.reverse(fll); 205 int m = Math.min(params.getMaxrefine(),fll.size()); 206 List<JointFragments> retlst = new ArrayList<JointFragments>(); 207 for ( int i = 0 ; i < m ; i++){ 208 JointFragments jf = fll.get(i); 209 retlst.add(jf); 210 } 211 212 return retlst.toArray(new JointFragments[retlst.size()]); 213 214 } 215 216 private boolean densityCheckOk(Atom[] aa1, Atom[] aa2, List<int[]> idxlist, 217 int p2i, int p2j, int l2, 218 double densityCutoff) throws StructureException { 219 JointFragments ftmp = new JointFragments(); 220 ftmp.setIdxlist(idxlist); 221 ftmp.add(p2i,p2j,0,l2); 222 AlternativeAlignment ali = new AlternativeAlignment(); 223 ali.apairs_from_idxlst(ftmp); 224 225 Atom[] aa3 = aa2.clone(); 226 227 int[] idx1 = ali.getIdx1(); 228 int[] idx2 = ali.getIdx2(); 229 230 Atom[] ca1subset = AlignUtils.getFragmentFromIdxList(aa1, idx1); 231 232 Atom[] ca2subset = AlignUtils.getFragmentFromIdxList(aa3,idx2); 233 234 double density = getDensity(ca1subset, ca2subset); 235 236 return density <= densityCutoff; 237 238 239 } 240 241 242 /** this is probably useless 243 * 244 * @param ca1subset 245 * @param ca2subset 246 * @return a double 247 * @throws StructureException 248 */ 249 private double getDensity(Atom[] ca1subset, Atom[] ca2subset ) throws StructureException{ 250 251 Atom centroid1 = Calc.getCentroid(ca1subset); 252 Atom centroid2 = Calc.getCentroid(ca2subset); 253 254 // get Average distance to centroid ... 255 256 double d1 = 0; 257 double d2 = 0; 258 259 for ( int i = 0 ; i < ca1subset.length;i++){ 260 double dd1 = Calc.getDistance(centroid1, ca1subset[i]); 261 double dd2 = Calc.getDistance(centroid2, ca2subset[i]); 262 263 d1 += dd1; 264 d2 += dd2; 265 266 } 267 268 double avd1 = d1 / ca1subset.length; 269 double avd2 = d2 / ca2subset.length; 270 271 return Math.min(avd1,avd2); 272 } 273 274 private double rmsCheck(Atom[] a1, Atom[] a2,List<int[]> idxlist, int p2i, int p2j, int l2) throws StructureException { 275 276 //System.out.println(dd); 277 // check if a joint fragment has low rms ... 278 JointFragments ftmp = new JointFragments(); 279 ftmp.setIdxlist(idxlist); 280 ftmp.add(p2i,p2j,0,l2); 281 Atom[] a3 = new Atom[a2.length]; 282 for (int i=0;i < a2.length;i++){ 283 a3[i] = (Atom)a2[i].clone(); 284 } 285 return getRMS(a1,a3,ftmp); 286 } 287 288 /** 289 * Get the RMS of the JointFragments pair frag 290 * 291 * @param ca1 the array of all atoms of structure1 292 * @param ca2 the array of all atoms of structure1 293 * @param frag the JointFragments object that contains the list of identical positions 294 * @return the rms 295 */ 296 public static double getRMS(Atom[] ca1, Atom[]ca2,JointFragments frag) throws StructureException { 297 // now svd ftmp and check if the rms is < X ... 298 AlternativeAlignment ali = new AlternativeAlignment(); 299 ali.apairs_from_idxlst(frag); 300 double rms = 999; 301 302 int[] idx1 = ali.getIdx1(); 303 int[] idx2 = ali.getIdx2(); 304 305 Atom[] ca1subset = AlignUtils.getFragmentFromIdxList(ca1, idx1); 306 307 Atom[] ca2subset = AlignUtils.getFragmentFromIdxList(ca2,idx2); 308 309 ali.calculateSuperpositionByIdx(ca1,ca2); 310 311 Matrix rot = ali.getRotationMatrix(); 312 Atom atom = ali.getShift(); 313 314 for (Atom a : ca2subset) { 315 Calc.rotate(a, rot); 316 Calc.shift(a, atom); 317 } 318 319 rms = Calc.rmsd(ca1subset,ca2subset); 320 321 return rms; 322 } 323 324 public boolean angleCheckOk(FragmentPair a, FragmentPair b, float distcutoff){ 325 326 double dist = -999; 327 328 Atom v1 = a.getUnitv(); 329 Atom v2 = b.getUnitv(); 330 dist = Calc.getDistance(v1,v2); 331 332 return dist <= distcutoff; 333 } 334 335 private boolean distanceCheckOk(FragmentPair a, FragmentPair b, float fragCompatDist){ 336 337 double dd ; 338 339 Atom c1i = a.getCenter1(); 340 Atom c1j = b.getCenter1(); 341 Atom c2i = a.getCenter2(); 342 Atom c2j = b.getCenter2(); 343 dd = Calc.getDistance(c1i,c1j) - Calc.getDistance(c2i,c2j); 344 345 346 if ( dd < 0) dd = -dd; 347 return dd <= fragCompatDist; 348 349 } 350 351 /** 352 * Calculate the pairwise compatibility of fpairs. 353 * Iterates through a list of fpairs and joins them if 354 * they have compatible rotation and translation parameters. 355 * @param fraglst FragmentPair[] array 356 * @param angleDiff angle cutoff 357 * @param fragCompatDist distance cutoff 358 * @param maxRefine max number of solutions to keep 359 * @return JointFragments[] 360 */ 361 public JointFragments[] frag_pairwise_compat(FragmentPair[] fraglst, int angleDiff, float fragCompatDist, int maxRefine ){ 362 363 364 FragmentPair[] tmpfidx = new FragmentPair[fraglst.length]; 365 for ( int i=0 ; i < fraglst.length; i++){ 366 tmpfidx[i] = (FragmentPair)fraglst[i].clone(); 367 } 368 369 int n = tmpfidx.length; 370 371 //indicator if a fragment was already used 372 int[] used = new int[n]; 373 374 //the final list of joined fragments stores as apairs 375 List<JointFragments> fll = new ArrayList<JointFragments>(); 376 377 double adiff = angleDiff * Math.PI / 180d; 378 logger.debug("addiff" + adiff); 379 //distance between two unit vectors with angle adiff 380 double ddiff = Math.sqrt(2.0-2.0*Math.cos(adiff)); 381 logger.debug("ddiff" + ddiff); 382 383 // the fpairs in the flist have to be sorted with respect to their positions 384 385 while (tmpfidx.length > 0){ 386 int i = 0; 387 int j = 1; 388 used[i]=1; 389 JointFragments f = new JointFragments(); 390 391 int p1i = tmpfidx[i].getPos1(); 392 int p1j = tmpfidx[i].getPos2(); 393 394 int maxi = p1i+tmpfidx[i].getLength(); 395 396 f.add(p1i,p1j,0,tmpfidx[i].getLength()); 397 398 n = tmpfidx.length; 399 400 while ( j < n) { 401 int p2i = tmpfidx[j].getPos1(); 402 int p2j = tmpfidx[j].getPos2(); 403 int l2 = tmpfidx[j].getLength(); 404 if ( p2i > maxi) { 405 406 double dist = Calc.getDistance(tmpfidx[i].getUnitv(), tmpfidx[j].getUnitv()); 407 if ( dist < ddiff) { 408 409 // centroids have to be closer than fragCompatDist 410 double dd = Calc.getDistance(tmpfidx[i].getCenter1(),tmpfidx[j].getCenter1()) - 411 Calc.getDistance(tmpfidx[i].getCenter2(),tmpfidx[j].getCenter2()); 412 if ( dd < 0) 413 dd = -dd; 414 if ( dd < fragCompatDist){ 415 maxi = p2i+l2; 416 used[j]=1; 417 f.add(p2i,p2j,0,tmpfidx[j].getLength()); 418 } 419 } 420 421 422 } 423 j+=1; 424 } 425 426 int red = 0; 427 for (int k = 0 ; k < n ; k ++) { 428 if (used[k] == 0) { 429 tmpfidx[red] = tmpfidx[k]; 430 red+=1; 431 } 432 } 433 434 435 used = new int[n]; 436 tmpfidx = (FragmentPair[])resizeArray(tmpfidx,red); 437 438 fll.add(f); 439 } 440 441 442 Comparator<JointFragments> comp = new JointFragmentsComparator(); 443 Collections.sort(fll,comp); 444 Collections.reverse(fll); 445 int m = Math.min(maxRefine,fll.size()); 446 List<JointFragments> retlst = new ArrayList<JointFragments>(); 447 for ( int i = 0 ; i < m ; i++){ 448 JointFragments jf = fll.get(i); 449 retlst.add(jf); 450 } 451 452 return retlst.toArray(new JointFragments[retlst.size()]); 453 454 } 455 456 457 public void extendFragments(Atom[] ca1, Atom[] ca2 ,JointFragments[] fragments, StrucAligParameters params) throws StructureException { 458 459 for(JointFragments p : fragments){ 460 extendFragments(ca1, ca2, p, params); 461 } 462 463 } 464 465 public void extendFragments(Atom[] ca1, Atom[] ca2 , JointFragments fragments, StrucAligParameters params) throws StructureException { 466 467 List<int[]> pos = fragments.getIdxlist(); 468 469 int[] firstP = pos.get(0); 470 int pstart1 = firstP[0]; 471 int pstart2 = firstP[1]; 472 473 int[] lastP = pos.get(pos.size()-1); 474 int pend1 = lastP[0]; 475 int pend2 = lastP[1]; 476 477 boolean startReached = false; 478 boolean endReached = false; 479 480 while (! (startReached && endReached)){ 481 482 if ( ! (startReached) && ((pstart1 <= 0) || (pstart2 <= 0))) { 483 startReached = true; 484 485 } else { 486 pstart1--; 487 pstart2--; 488 } 489 490 if ( ! (endReached) && (( pend1 >= (ca1.length -1) ) || ( pend2 >= ca2.length -1 )) ){ 491 endReached = true; 492 } else { 493 pend1++; 494 pend2++; 495 } 496 497 if ( ! startReached){ 498 double newRms1 = testAdd(ca1, ca2, fragments,pstart1,pstart2); 499 500 if (newRms1 < 3.7 ) { 501 fragments.add(pstart1,pstart2,0,1); 502 } else { 503 startReached = true; 504 } 505 } 506 507 if( ! endReached){ 508 509 double newRms2 = testAdd(ca1, ca2, fragments, pend1, pend2); 510 if ( newRms2 < 3.7) { 511 fragments.add(pend1,pend2,0,1); 512 } else { 513 endReached = true; 514 } 515 } 516 517 } 518 519 } 520 521 522 private double testAdd(Atom[] ca1, Atom[] ca2, JointFragments fragments, int pstart1, int pstart2) throws StructureException { 523 524 JointFragments frag = new JointFragments(); 525 frag.setIdxlist(fragments.getIdxlist()); 526 frag.add(pstart1, pstart2, 0,1); 527 528 return FragmentJoiner.getRMS(ca1, ca2, frag); 529 530 } 531 532} 533 534 535 536 537class JointFragmentsComparator 538 implements Comparator<JointFragments>, Serializable { 539 private static final long serialVersionUID = 1; 540 541 @Override 542 public int compare(JointFragments one, JointFragments two) { 543 544 545 int s1 = one.getIdxlist().size(); 546 int s2 = two.getIdxlist().size(); 547 548 double rms1 = one.getRms(); 549 double rms2 = two.getRms(); 550 if ( s1 > s2 ) { 551 return 1 ; 552 } 553 554 else if ( s1 < s2){ 555 return -1; 556 } 557 else{ 558 if ( rms1 < rms2) 559 return 1; 560 if ( rms1 > rms2) 561 return -1; 562 return 0; 563 } 564 } 565 566 567}