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 025 026import org.biojava.nbio.structure.*; 027import org.biojava.nbio.structure.align.StrucAligParameters; 028import org.biojava.nbio.structure.align.helper.AligMatEl; 029import org.biojava.nbio.structure.align.helper.IndexPair; 030import org.biojava.nbio.structure.align.helper.JointFragments; 031import org.biojava.nbio.structure.geometry.Matrices; 032import org.biojava.nbio.structure.geometry.SuperPositions; 033import org.biojava.nbio.structure.jama.Matrix; 034import org.slf4j.Logger; 035import org.slf4j.LoggerFactory; 036 037import java.io.Serializable; 038import java.text.DecimalFormat; 039import java.util.List; 040 041 042import javax.vecmath.Matrix4d; 043 044/** 045 * Implements a class which handles one possible (alternative) solution. 046 * 047 * Alternative alignments arise from different seed 048 * alignments or seed FPairs. The AltAlg class contains methods 049 * for refinement (Dynamic Programming based) and filtering 050 * (i.e. removing probably wrongly matched APairs). In the refinement 051 * phase, different seed alignments can converge to the same solution. 052 * 053 * @author Andreas Prlic, 054 * @author Peter Lackner (original Python and C code) 055 * @since 3:04:26 PM 056 * @version %I% %G% 057 */ 058public class AlternativeAlignment implements Serializable{ 059 060 061 private static final long serialVersionUID = -6226717654562221241L; 062 063 private int[] idx1; 064 private int[] idx2; 065 private String[] pdbresnum1; 066 private String[] pdbresnum2; 067 //short[] alig1; 068 //short[] alig2; 069 070 private int nfrags; 071 private Atom center; 072 private Matrix rot; 073 private Atom tr; 074 075 076 // the scores... 077 private int gaps0; 078 private int eqr0; 079 private int rms0; 080 private int joined; 081 private int percId; 082 private int cluster; 083 private float score; 084 private IndexPair[] aligpath; 085 private int fromia; 086 private Matrix currentRotMatrix; 087 private Atom currentTranMatrix; 088 089 private double rms; 090 091 private Matrix distanceMatrix; 092 093 public static final Logger logger = LoggerFactory.getLogger(AlternativeAlignment.class); 094 095 096 public AlternativeAlignment() { 097 super(); 098 099 nfrags = 0; 100 aligpath = new IndexPair[0]; 101 //alig1 = new short[0]; 102 //alig2 = new short[0]; 103 idx1 = new int[0]; 104 idx2 = new int[0]; 105 106 center = new AtomImpl(); 107 rot = null; 108 tr = new AtomImpl(); 109 eqr0 = -99; 110 rms0 = 99; 111 joined = 0; 112 gaps0 = -99; 113 fromia = 0; 114 115 currentRotMatrix = new Matrix(0,0); 116 currentTranMatrix = new AtomImpl(); 117 118 distanceMatrix = new Matrix(0,0); 119 } 120 121 122 123 124 /** print the idx positions of this alignment 125 * 126 * @return a String representation 127 */ 128 @Override 129 public String toString(){ 130 DecimalFormat d2 = new DecimalFormat(); 131 // the result can be localized. To change this and enforce UK local do... 132 //(DecimalFormat)NumberFormat.getInstance(java.util.Locale.UK); 133 d2.setMaximumIntegerDigits(3); 134 d2.setMinimumFractionDigits(2); 135 d2.setMaximumFractionDigits(2); 136 StringBuffer s = new StringBuffer(); 137 s.append("#" + getAltAligNumber() + 138 " cluster:" + cluster + 139 " eqr:" + getEqr() + 140 " rmsd:" + d2.format(getRmsd()) + 141 " %id:" + getPercId() + 142 " gaps:" + getGaps() + 143 " score:" + d2.format(score) ); 144 145 return s.toString(); 146 } 147 148 149 /** get the number of the cluster this alignment belongs to 150 * 151 * @return an int giving the number of the cluster 152 */ 153 public int getCluster() { 154 return cluster; 155 } 156 157 158 159 /** set the number of the cluster this alignment belongs to. 160 * All alignments in a cluster are quite similar. 161 * 162 * @param cluster the number of the cluster 163 */ 164 public void setCluster(int cluster) { 165 this.cluster = cluster; 166 } 167 168 169 170 171 public double getRmsd() { 172 return rms; 173 } 174 175 /** the rms in the structurally equivalent regions 176 * 177 * @param rms 178 */ 179 public void setRms(double rms) { 180 this.rms = rms; 181 } 182 183 /** the alignment score 184 * 185 * @return the score of this alignment 186 */ 187 public float getScore() { 188 return score; 189 } 190 191 public void setScore(float score) { 192 this.score = score; 193 } 194 195 196 /** return the number of gaps in this alignment 197 * 198 * @return the number of Gaps 199 */ 200 public int getGaps(){ 201 return gaps0; 202 } 203 204 /** returns the number of euqivalent residues in this alignment 205 * 206 * @return the number of equivalent residues 207 */ 208 public int getEqr(){ 209 return eqr0 ; 210 } 211 212 /** the positions of the structure equivalent positions in atom set 1 213 * 214 * @return the array of the positions 215 */ 216 public int[] getIdx1(){ 217 return idx1; 218 } 219 220 /** the positions of the structure equivalent atoms in atom set 2 221 * 222 * @return the array of the positions 223 */ 224 public int[] getIdx2(){ 225 return idx2; 226 } 227 228 public int getPercId() { 229 return percId; 230 } 231 232 public void setPercId(int percId) { 233 this.percId = percId; 234 } 235 236 /** Set apairs according to a seed position. 237 * 238 * @param l 239 * @param i 240 * @param j 241 */ 242 public void apairs_from_seed(int l,int i, int j){ 243 aligpath = new IndexPair[l]; 244 idx1 = new int[l]; 245 idx2 = new int[l]; 246 for (int x = 0 ; x < l ; x++) { 247 idx1[x]=i+x; 248 idx2[x]=j+x; 249 aligpath[x] = new IndexPair((short)(i+x),(short)(j+x)); 250 } 251 } 252 253 /** Set apairs according to a list of (i,j) tuples. 254 * 255 * @param jf a JoingFragment 256 */ 257 public void apairs_from_idxlst(JointFragments jf) { 258 List<int[]> il = jf.getIdxlist(); 259 //System.out.println("Alt Alig apairs_from_idxlst"); 260 261 aligpath = new IndexPair[il.size()]; 262 idx1 = new int[il.size()]; 263 idx2 = new int[il.size()]; 264 for (int i =0 ; i < il.size();i++) { 265 int[] p = il.get(i); 266 //System.out.print(" idx1 " + p[0] + " idx2 " + p[1]); 267 idx1[i] = p[0]; 268 idx2[i] = p[1]; 269 aligpath[i] = new IndexPair((short)p[0],(short)p[1]); 270 271 } 272 eqr0 = idx1.length; 273 //System.out.println("eqr " + eqr0); 274 gaps0 = count_gaps(idx1,idx2); 275 276 } 277 278 279 /** returns the sequential number of this alternative alignment 280 * 281 * @return the sequential number of this alternative alignment 282 */ 283 public int getAltAligNumber() { 284 return fromia; 285 } 286 287 public void setAltAligNumber(int fromia) { 288 this.fromia = fromia; 289 } 290 291 /** rotate and shift atoms with currentRotMatrix and current Tranmatrix 292 * 293 * @param ca 294 */ 295 private void rotateShiftAtoms(Atom[] ca){ 296 297 for (int i = 0 ; i < ca.length; i++){ 298 Atom c = ca[i]; 299 300 Calc.rotate(c,currentRotMatrix); 301 Calc.shift(c,currentTranMatrix); 302 //System.out.println("after " + c); 303 ca[i] = c; 304 } 305 //System.out.println("after " + ca[0]); 306 } 307 308 public void finish(StrucAligParameters params,Atom[]ca1,Atom[]ca2) throws StructureException{ 309 310 Atom[] ca3 = new Atom[ca2.length]; 311 for ( int i = 0 ; i < ca2.length;i++){ 312 ca3[i] = (Atom) ca2[i].clone(); 313 314 } 315 // do the inital superpos... 316 317 super_pos_alig(ca1,ca3,idx1,idx2,true); 318 rotateShiftAtoms(ca3); 319 320 calcScores(ca1,ca2); 321 logger.debug("eqr " + eqr0 + " " + gaps0 + " " +idx1[0] + " " +idx1[1]); 322 323 getPdbRegions(ca1,ca2); 324 325 } 326 327 public static Matrix getDistanceMatrix(Atom[] ca1, Atom[]ca2){ 328 329 int r = ca1.length; 330 int c = ca2.length; 331 332 Matrix out = new Matrix(r,c); 333 334 for (int i=0; i<r; i++) { 335 Atom a1 = ca1[i]; 336 for (int j=0;j<c;j++){ 337 Atom b1 = ca2[j]; 338 339 340 double d = Calc.getDistance(a1,b1); 341 out.set(i,j,d); 342 343 } 344 } 345 return out; 346 } 347 348 private Alignable getInitalStrCompAlignment( 349 Atom[] ca1, 350 Atom[]ca2, 351 StrucAligParameters params 352 ){ 353 354 int rows = ca1.length; 355 int cols = ca2.length; 356 357 float gapOpen = params.getGapOpen(); 358 float gapExtension = params.getGapExtension(); 359 float co = params.getCreate_co(); 360 361 Alignable al = new StrCompAlignment(rows,cols); 362 al.setGapExtCol(gapExtension); 363 al.setGapExtRow(gapExtension); 364 al.setGapOpenCol(gapOpen); 365 al.setGapOpenRow(gapOpen); 366 367 AligMatEl[][] aligmat = al.getAligMat(); 368 369 //System.out.println(rows + " " + cols ); 370 371 aligmat[0][cols] = new AligMatEl(); 372 aligmat[rows][0] = new AligMatEl(); 373 aligmat[rows][cols] = new AligMatEl(); 374 375 for (int i = 0; i < rows; i++) { 376 Atom a1 = ca1[i]; 377 378 aligmat[i][0] = new AligMatEl(); 379 //aligmat[i][cols] = new AligMatEl(); 380 381 for (int j = 0; j < cols; j++) { 382 //System.out.println("index " + i + " " + j); 383 aligmat[0][j] = new AligMatEl(); 384 385 Atom b1 = ca2[j]; 386 double d = 999; 387 388 d = Calc.getDistance(a1,b1); 389 390 391 AligMatEl e = new AligMatEl(); 392 if (d > co) { 393 e.setValue(0); 394 } else { 395 double s = 2 * co / ( 1 + ( d/co) * (d/co)) - co; 396 e.setValue((int) Math.round(Gotoh.ALIGFACTOR * s)); 397 } 398 aligmat[i+1][j+1] = e; 399 } 400 } 401 402 return al; 403 } 404 405 /*private Alignable getAlignableFromSimMatrix (Matrix sim,StrucAligParameters params){ 406 //System.out.println("align_NPE"); 407 408 float gapOpen = params.getGapOpen(); 409 float gapExtension = params.getGapExtension(); 410 411 int rows = sim.getRowDimension(); 412 int cols = sim.getColumnDimension(); 413 414 Alignable al = new StrCompAlignment(rows,cols); 415 al.setGapExtCol(gapExtension); 416 al.setGapExtRow(gapExtension); 417 al.setGapOpenCol(gapOpen); 418 al.setGapOpenRow(gapOpen); 419 //System.out.println("size of aligmat: " + rows+1 + " " + cols+1); 420 //AligMatEl[][] aligmat = new AligMatEl[rows+1][cols+1]; 421 AligMatEl[][] aligmat = al.getAligMat(); 422 423 for (int i = 0; i < rows; i++) { 424 for (int j = 0; j < cols; j++) { 425 426 int e=0; 427 //if ( ( i < rows) && 428 // ( j < cols)) { 429 //TODO: the ALIGFACTOR calc should be hidden in Gotoh!! 430 431 e = (int)Math.round(Gotoh.ALIGFACTOR * sim.get(i,j)); 432 //} 433 //System.out.println(e); 434 //AligMatEl am = new AligMatEl(); 435 aligmat[i+1][j+1].setValue(e); 436 //.setValue(e); 437 //am.setTrack(new IndexPair((short)-99,(short)-99)); 438 //igmat[i+1][j+1] = am; 439 440 } 441 } 442 //al.setAligMat(aligmat); 443 444 445 return al; 446 }*/ 447 /** Refinement procedure based on superposition and dynamic programming. 448 * Performs an iterative refinement. Several methods apply such a procedure, 449 * e.g. CE or ProSup. Here we additionally test for circular permutation, 450 * which are in the same frame of superposition as the optimal alignment. 451 * This feature may be switched off by setting permsize to -1. 452 * 453 * @param params the parameters 454 * @param ca1 atoms of structure 1 455 * @param ca2 atoms of structure 2 456 * @throws StructureException 457 */ 458 459 public void refine(StrucAligParameters params,Atom[]ca1,Atom[]ca2) throws StructureException{ 460 // System.out.println("refine Alternative Alignment #"+ getAltAligNumber()+" l1:" + ca1.length + " l2:" + ca2.length); 461 // for ( int i= 0 ; i < idx1.length;i++){ 462 // System.out.println(idx1[i] + " " + idx2[i]); 463 // } 464 465 // avoid any operations on the original Atoms ... 466 Atom[] ca3 = new Atom[ca2.length]; 467 for ( int i = 0 ; i < ca2.length;i++){ 468 ca3[i] = (Atom) ca2[i].clone(); 469 470 } 471 //Atom[] ca3 = (Atom[]) ca2.clone(); 472 473 // do the inital superpos... 474 475 super_pos_alig(ca1,ca3,idx1,idx2,false); 476 477 //JmolDisplay jmol1 = new JmolDisplay(ca1,ca2,"after first alig " ); 478 //jmol1.selectEquiv(idx1,idx2); 479 //int[] jmoltmp1 = idx1; 480 //int[] jmoltmp2 = idx2; 481 482 int lenalt =idx1.length; 483 int lenneu = aligpath.length; 484 //Matrix rmsalt = currentRotMatrix; 485 486 487 488 int ml = Math.max(ca1.length, ca3.length); 489 idx1 = new int[ml]; 490 idx2 = new int[ml]; 491 492 //int m1l = ca1.length; 493 //int m2l = ca3.length; 494 495 //Matrix sim = new Matrix(ca1.length,ca3.length,0); 496 497 int maxiter = params.getMaxIter(); 498 for (int iter = 0 ; iter< maxiter; iter++){ 499 500 float subscore = 0.0f; 501 502 rotateShiftAtoms(ca3); 503 504 // do the dynamic alignment... 505 Alignable ali = getInitalStrCompAlignment(ca1,ca3, params); 506 507 new Gotoh(ali); 508 509 this.score = ali.getScore(); 510 //System.out.println("score " + score); 511 IndexPair[] path = ali.getPath(); 512 513 int pathsize = ali.getPathSize(); 514 515 516 // now for a superimposable permutation. First we need to check the size and 517 // position of the quadrant left out by the optimal path 518 int firsta = path[0].getRow(); 519 int firstb = path[0].getCol(); 520 int lasta = path[pathsize-1].getRow(); 521 int lastb = path[pathsize-1].getCol(); 522 523 int quada_beg, quada_end, quadb_beg,quadb_end; 524 525 if ( firsta == 0){ 526 quada_beg = lasta +1; 527 quada_end = ca1.length -1; 528 quadb_beg = 0; 529 quadb_end = firstb-1; 530 } else { 531 quada_beg = 0; 532 quada_end = firsta-1; 533 quadb_beg = lastb+1; 534 quadb_end = ca3.length -1; 535 } 536 537 // if the unaligned quadrant is larger than permsize 538 int permsize = params.getPermutationSize(); 539 540 if ( permsize > 0 && 541 (quada_end - quada_beg >= permsize) && 542 ( quadb_end - quadb_beg >= permsize)) { 543 AligMatEl[][] aligmat = ali.getAligMat(); 544 Matrix submat = new Matrix(quada_end-quada_beg, quadb_end - quadb_beg) ; 545 //then we copy the score values of the quadrant into a new matrix 546 547 //System.out.println(quada_beg + " " + quada_end + " " +quadb_beg+" " + quadb_end); 548 //System.out.println(sim.getRowDimension() + " " + sim.getColumnDimension()); 549 Atom[] tmp1 = new Atom[quada_end - quada_beg ]; 550 Atom[] tmp2 = new Atom[quadb_end - quadb_beg ]; 551 552 for ( int s = quada_beg; s < quada_end; s++){ 553 //System.out.println(s + " " +( quada_end-s)); 554 tmp1[quada_end - s -1] = ca1[s]; 555 for ( int t = quadb_beg; t < quadb_end; t++){ 556 if ( s == quada_beg ) 557 tmp2[quadb_end - t -1 ] = ca3[t]; 558 559 //System.out.println(s+" " +t); 560 //double val = sim.get(s,t); 561 double val = aligmat[s][t].getValue(); 562 submat.set(s-quada_beg,t-quadb_beg,val); 563 } 564 } 565 // and perform a dp alignment again. (Note, that we fix the superposition). 566 567 //Alignable subali = getAlignableFromSimMatrix(submat,params); 568 Alignable subali = getInitalStrCompAlignment(tmp1, tmp2, params); 569 subscore = subali.getScore(); 570 //System.out.println("join score" + score + " " + subscore); 571 this.score = score+subscore; 572 573 IndexPair[] subpath = subali.getPath(); 574 int subpathsize = subali.getPathSize(); 575 for ( int p=0;p<subpath.length;p++){ 576 IndexPair sp = subpath[p]; 577 sp.setRow((short)(sp.getRow()+quada_beg)); 578 sp.setCol((short)(sp.getCol()+quadb_beg)); 579 } 580 581 // now we join the two solutions 582 if ( subpathsize > permsize){ 583 IndexPair[] wholepath = new IndexPair[pathsize+subpathsize]; 584 for ( int t=0; t < pathsize; t++){ 585 wholepath[t] = path[t]; 586 } 587 for ( int t=0 ; t < subpathsize; t++){ 588 wholepath[t+pathsize] = subpath[t]; 589 } 590 pathsize += subpathsize; 591 path = wholepath; 592 ali.setPath(path); 593 ali.setPathSize(pathsize); 594 } 595 596 } 597 598 599 int j =0 ; 600 //System.out.println("pathsize,path " + pathsize + " " + path.length); 601 602 for (int i=0; i < pathsize; i++){ 603 int x = path[i].getRow(); 604 int y = path[i].getCol(); 605 //System.out.println(x+ " " + y); 606 double d = Calc.getDistance(ca1[x], ca3[y]); 607 //double d = dismat.get(x,y); 608 // now we apply the evaluation distance cutoff 609 if ( d < params.getEvalCutoff()){ 610 //System.out.println("j:"+j+ " " + x+" " + y + " " + d ); 611 612 idx1[j] = x; 613 idx2[j] = y; 614 j+=1; 615 } 616 } 617 618 lenneu = j; 619 //System.out.println("lenalt,neu:" + lenalt + " " + lenneu); 620 621 622 int[] tmpidx1 = new int[j]; 623 int[] tmpidx2 = new int[j]; 624 //String idx1str ="idx1: "; 625 //String idx2str ="idx2: "; 626 for (int i = 0 ; i<j;i++){ 627 //idx1str += idx1[i]+ " "; 628 //idx2str += idx2[i]+ " "; 629 tmpidx1[i] = idx1[i]; 630 tmpidx2[i] = idx2[i]; 631 } 632 //System.out.println(idx1str); 633 //System.out.println(idx2str); 634 //jmol.selectEquiv(idx1,idx2); 635 636 //if ( j > 0) 637 // System.out.println("do new superimpos " + idx1.length + " " + idx2.length + " " + idx1[0] + " " + idx2[0]); 638 super_pos_alig(ca1,ca3,tmpidx1,tmpidx2,false); 639 640 this.aligpath = path; 641 642 if (lenneu == lenalt) 643 break; 644 645 } 646 647 648 649 idx1 = (int[]) FragmentJoiner.resizeArray(idx1,lenneu); 650 idx2 = (int[]) FragmentJoiner.resizeArray(idx2,lenneu); 651 // new ... get rms... 652 // now use the original atoms to get the rotmat relative to the original structure... 653 super_pos_alig(ca1,ca2,idx1,idx2,true); 654 eqr0 = idx1.length; 655 gaps0 = count_gaps(idx1,idx2); 656 657 getPdbRegions(ca1,ca2); 658 //System.out.println("eqr " + eqr0 + " aligpath len:"+aligpath.length+ " gaps:" + gaps0 + " score " + score); 659 } 660 661 private void getPdbRegions(Atom[] ca1, Atom[] ca2){ 662 pdbresnum1 = new String[idx1.length]; 663 pdbresnum2 = new String[idx2.length]; 664 665 for (int i =0 ; i < idx1.length;i++){ 666 Atom a1 = ca1[idx1[i]]; 667 Atom a2 = ca2[idx2[i]]; 668 Group p1 = a1.getGroup(); 669 Group p2 = a2.getGroup(); 670 Chain c1 = p1.getChain(); 671 Chain c2 = p2.getChain(); 672 673 String cid1 = c1.getId(); 674 String cid2 = c2.getId(); 675 676 String pdb1 = p1.getResidueNumber().toString(); 677 String pdb2 = p2.getResidueNumber().toString(); 678 679 680 if ( ! cid1.equals(" ")) 681 pdb1 += ":" + cid1; 682 683 684 if ( ! cid2.equals(" ")) 685 pdb2 += ":" + cid2; 686 687 688 pdbresnum1[i] = pdb1; 689 pdbresnum2[i] = pdb2; 690 } 691 } 692 693 694 public String[] getPDBresnum1() { 695 return pdbresnum1; 696 } 697 698 699 700 701 public void setPDBresnum1(String[] pdbresnum1) { 702 this.pdbresnum1 = pdbresnum1; 703 } 704 705 706 707 708 public String[] getPDBresnum2() { 709 return pdbresnum2; 710 } 711 712 713 714 715 public void setPDBresnum2(String[] pdbresnum2) { 716 this.pdbresnum2 = pdbresnum2; 717 } 718 719 720 721 722 /** 723 * Count the number of gaps in an alignment represented by idx1,idx2. 724 * 725 * @param i1 726 * @param i2 727 * @return the number of gaps in this alignment 728 */ 729 private int count_gaps(int[] i1, int[] i2){ 730 731 int i0 = i1[0]; 732 int j0 = i2[0]; 733 int gaps = 0; 734 for (int i =1 ; i<i1.length;i++ ){ 735 if ( Math.abs(i1[i]-i0) != 1 || 736 ( Math.abs(i2[i]-j0) != 1)){ 737 gaps +=1; 738 } 739 i0 = i1[i]; 740 j0 = i2[i]; 741 } 742 743 return gaps; 744 } 745 746 747 748 public void calculateSuperpositionByIdx(Atom[] ca1, Atom[] ca2)throws StructureException { 749 750 super_pos_alig(ca1,ca2,idx1,idx2,false); 751 752 } 753 754 /** Superimposes two molecules according to residue index list idx1 and idx2. 755 * Does not change the original coordinates. 756 * as an internal result the rotation matrix and shift vectors for are set 757 * 758 * @param ca1 Atom set 1 759 * @param ca2 Atom set 2 760 * @param idx1 idx positions in set1 761 * @param idx2 idx positions in set2 762 * @param getRMS a flag if the RMS should be calculated 763 * @throws StructureException 764 */ 765 766 767 768 private void super_pos_alig(Atom[]ca1,Atom[]ca2,int[] idx1, int[] idx2, boolean getRMS) throws StructureException{ 769 770 //System.out.println("superpos alig "); 771 Atom[] ca1subset = new Atom[idx1.length]; 772 Atom[] ca2subset = new Atom[idx2.length]; 773 774 for (int i = 0 ; i < idx1.length;i++){ 775 //System.out.println("idx1 "+ idx1[i] + " " + idx2[i]); 776 int pos1 = idx1[i]; 777 int pos2 = idx2[i]; 778 779 ca1subset[i] = ca1[pos1]; 780 ca2subset[i] = (Atom) ca2[pos2].clone(); 781 } 782 783 Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1subset), 784 Calc.atomsToPoints(ca2subset)); 785 this.currentRotMatrix = Matrices.getRotationJAMA(trans); 786 this.currentTranMatrix = Calc.getTranslationVector(trans); 787 //currentRotMatrix.print(3,3); 788 if ( getRMS) { 789 rotateShiftAtoms(ca2subset); 790 this.rms = Calc.rmsd(ca1subset,ca2subset); 791 } 792 793 794 } 795 796 797 /** returns the rotation matrix that needs to be applied to structure 2 to rotate on structure 1 798 * 799 * @return the rotation Matrix 800 */ 801 public Matrix getRotationMatrix(){ 802 return currentRotMatrix; 803 } 804 805 /** returns the shift vector that has to be applied on structure to to shift on structure one 806 * 807 * @return the shift vector 808 */ 809 public Atom getShift(){ 810 return currentTranMatrix; 811 } 812 813 814 815 /** calculates scores for this alignment ( %id ) 816 * @param ca1 set of Atoms for molecule 1 817 * @param ca2 set of Atoms for molecule 2 818 819 */ 820 public void calcScores(Atom[] ca1, Atom[] ca2){ 821 eqr0 = idx1.length; 822 gaps0 = count_gaps(idx1,idx2); 823 824 percId = 0; 825 // calc the % id 826 for (int i=0 ; i< idx1.length; i++){ 827 Atom a1 = ca1[idx1[i]]; 828 Atom a2 = ca2[idx2[i]]; 829 830 Group g1 = a1.getGroup(); 831 Group g2 = a2.getGroup(); 832 if ( g1.getPDBName().equals(g2.getPDBName())){ 833 percId++; 834 } 835 836 } 837 } 838 839 840 841 /** create an artifical Structure object that contains the two 842 * structures superimposed onto each other. Each structure is in a separate model. 843 * Model 1 is structure 1 and Model 2 is structure 2. 844 * 845 * @param s1 the first structure. its coordinates will not be changed 846 * @param s2 the second structure, it will be cloned and the cloned coordinates will be rotated according to the alignment results. 847 * @return composite structure containing the 2 aligned structures as a models 1 and 2 848 */ 849 public Structure getAlignedStructure(Structure s1, Structure s2){ 850 // do not change the original coords .. 851 Structure s3 = s2.clone(); 852 853 currentRotMatrix.print(3,3); 854 855 Calc.rotate(s3, currentRotMatrix); 856 Calc.shift( s3, currentTranMatrix); 857 858 Structure newpdb = new StructureImpl(); 859 newpdb.setPDBCode("Java"); 860 newpdb.setName("Aligned with BioJava"); 861 862 863 newpdb.addModel(s1.getChains(0)); 864 newpdb.addModel(s3.getChains(0)); 865 866 return newpdb; 867 868 } 869 870 /** converts the alignment to a PDB file 871 * each of the structures will be represented as a model. 872 * 873 * 874 * @param s1 875 * @param s2 876 * @return a PDB file as a String 877 */ 878 public String toPDB(Structure s1, Structure s2){ 879 880 881 Structure newpdb = getAlignedStructure(s1, s2); 882 883 return newpdb.toPDB(); 884 } 885 886 887 888 /** The distance matrix this alignment is based on 889 * 890 * @return a Matrix object. 891 */ 892 public Matrix getDistanceMatrix() { 893 return distanceMatrix; 894 } 895 896 897 898 /** The distance matrix this alignment is based on 899 * 900 * @param distanceMatrix 901 */ 902 public void setDistanceMatrix(Matrix distanceMatrix) { 903 this.distanceMatrix = distanceMatrix; 904 } 905 906 907 908 909 public IndexPair[] getPath() { 910 return aligpath; 911 } 912 913 914 915} 916