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