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 May 21, 2006 021 * 022 */ 023package org.biojava.nbio.structure.align; 024 025import org.biojava.nbio.structure.*; 026import org.biojava.nbio.structure.align.ce.GuiWrapper; 027import org.biojava.nbio.structure.align.helper.AlignTools; 028import org.biojava.nbio.structure.align.helper.JointFragments; 029import org.biojava.nbio.structure.align.pairwise.*; 030import org.biojava.nbio.structure.io.PDBFileParser; 031import org.biojava.nbio.structure.io.PDBFileReader; 032import org.biojava.nbio.structure.jama.Matrix; 033import org.slf4j.Logger; 034import org.slf4j.LoggerFactory; 035 036import java.io.FileOutputStream; 037import java.io.InputStream; 038import java.io.PrintStream; 039import java.util.ArrayList; 040import java.util.Collections; 041import java.util.Comparator; 042import java.util.List; 043 044 045/** 046 * Perform a pairwise protein structure superimposition. 047 * 048 * <p> 049 * The algorithm is a distance matrix based, rigid body protein structure superimposition. 050 * It is based on a variation of the PSC++ algorithm provided by Peter Lackner 051 * (Peter.Lackner@sbg.ac.at, personal communication) . 052 * </p> 053 * 054 * 055 * 056 * <h2>Example</h2> 057 * <pre> 058 * public void run(){ 059 060 // first load two example structures 061 {@link InputStream} inStream1 = this.getClass().getResourceAsStream("/files/5pti.pdb"); 062 {@link InputStream} inStream2 = this.getClass().getResourceAsStream("/files/1tap.pdb"); 063 064 {@link Structure} structure1 = null; 065 {@link Structure} structure2 = null; 066 067 {@link PDBFileParser} pdbpars = new {@link PDBFileParser}(); 068 structure1 = pdbpars.parsePDBFile(inStream1) ; 069 structure2 = pdbpars.parsePDBFile(inStream2); 070 071 072 // calculate structure superimposition for two complete structures 073 {@link StructurePairAligner} aligner = new {@link StructurePairAligner}(); 074 075 076 // align the full 2 structures with default parameters. 077 // see StructurePairAligner for more options and how to align 078 // any set of Atoms 079 aligner.align(structure1,structure2); 080 081 {@link AlternativeAlignment}[] aligs = aligner.getAlignments(); 082 {@link AlternativeAlignment} a = aligs[0]; 083 System.out.println(a); 084 085 //display the alignment in Jmol 086 087 // first get an artificial structure for the alignment 088 {@link Structure} artificial = a.getAlignedStructure(structure1, structure2); 089 090 091 // and then send it to Jmol (only will work if Jmol is in the Classpath) 092 093 BiojavaJmol jmol = new BiojavaJmol(); 094 jmol.setTitle(artificial.getName()); 095 jmol.setStructure(artificial); 096 097 // color the two structures 098 099 100 jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " + 101 "select not protein and not solvent; spacefill on;"); 102 jmol.evalString("select *"+"/1 ; color red; model 1; "); 103 104 105 // now color the equivalent residues ... 106 107 String[] pdb1 = a.getPDBresnum1(); 108 for (String res : pdb1 ){ 109 jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;"); 110 } 111 112 jmol.evalString("select *"+"/2; color blue; model 2;"); 113 String[] pdb2 = a.getPDBresnum2(); 114 for (String res :pdb2 ){ 115 jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;"); 116 } 117 118 119 // now show both models again. 120 jmol.evalString("model 0;"); 121 122 } 123 * </pre> 124 * 125 * 126 * 127 * @author Andreas Prlic 128 * @author Peter Lackner 129 * @since 1.4 130 * @version %I% %G% 131 */ 132public class StructurePairAligner { 133 134 private final static Logger logger = LoggerFactory.getLogger(StructurePairAligner.class); 135 136 AlternativeAlignment[] alts; 137 Matrix distanceMatrix; 138 StrucAligParameters params; 139 FragmentPair[] fragPairs; 140 141 List<AlignmentProgressListener> listeners = new ArrayList<AlignmentProgressListener>(); 142 143 public StructurePairAligner() { 144 super(); 145 params = StrucAligParameters.getDefaultParameters(); 146 reset(); 147 alts = new AlternativeAlignment[0]; 148 distanceMatrix = new Matrix(0,0); 149 } 150 151 public void addProgressListener(AlignmentProgressListener li){ 152 listeners.add(li); 153 } 154 155 public void clearListeners(){ 156 listeners.clear(); 157 } 158 159 160 /** example usage of this class 161 * 162 * @param args 163 */ 164 public static void main(String[] args) throws Exception { 165 // UPDATE THE FOLLOWING LINES TO MATCH YOUR SETUP 166 167 PDBFileReader pdbr = new PDBFileReader(); 168 pdbr.setPath("/Users/andreas/WORK/PDB/"); 169 170 171 //String pdb1 = "1crl"; 172 //String pdb2 = "1ede"; 173 174 String pdb1 = "1buz"; 175 String pdb2 = "1ali"; 176 String outputfile = "/tmp/alig_"+pdb1+"_"+pdb2+".pdb"; 177 178 // NO NEED TO DO CHANGE ANYTHING BELOW HERE... 179 180 StructurePairAligner sc = new StructurePairAligner(); 181 182 183 // step1 : read molecules 184 185 logger.info("aligning {} vs. {}", pdb1, pdb2); 186 187 Structure s1 = pdbr.getStructureById(pdb1); 188 Structure s2 = pdbr.getStructureById(pdb2); 189 190 // step 2 : do the calculations 191 sc.align(s1,s2); 192 193 194 AlternativeAlignment[] aligs = sc.getAlignments(); 195 196 //cluster similar results together 197 ClusterAltAligs.cluster(aligs); 198 199 200 // print the result: 201 // the AlternativeAlignment object gives also access to rotation matrices / shift vectors. 202 for (AlternativeAlignment aa : aligs) { 203 logger.info("Alternative Alignment: ", aa); 204 } 205 206 207 208 // convert AlternativeAlignemnt 1 to PDB file, so it can be opened with a viewer (e.g. Jmol, Rasmol) 209 210 if ( aligs.length > 0) { 211 AlternativeAlignment aa1 =aligs[0]; 212 String pdbstr = aa1.toPDB(s1,s2); 213 214 logger.info("writing alignment to {}", outputfile); 215 FileOutputStream out= new FileOutputStream(outputfile); 216 PrintStream p = new PrintStream( out ); 217 218 p.println (pdbstr); 219 220 p.close(); 221 out.close(); 222 } 223 224 225 // display the alignment in Jmol 226 // only will work if Jmol is in the Classpath 227 228 if ( aligs.length > 0) { 229 230 if (! GuiWrapper.isGuiModuleInstalled()){ 231 logger.error("Could not find structure-gui modules in classpath, please install first!"); 232 } 233 234 } 235 236 237 } 238 239 240 241 private void reset(){ 242 alts = new AlternativeAlignment[0]; 243 distanceMatrix = new Matrix(0,0); 244 fragPairs = new FragmentPair[0]; 245 246 } 247 248 249 /** get the results of step 1 - the FragmentPairs used for seeding the alignment 250 * @return a FragmentPair[] array 251 */ 252 253 public FragmentPair[] getFragmentPairs() { 254 return fragPairs; 255 } 256 257 258 259 public void setFragmentPairs(FragmentPair[] fragPairs) { 260 this.fragPairs = fragPairs; 261 } 262 263 264 /** return the alternative alignments that can be found for the two structures 265 * 266 * @return AlternativeAlignment[] array 267 */ 268 public AlternativeAlignment[] getAlignments() { 269 return alts; 270 } 271 272 /** return the difference of distance matrix between the two structures 273 * 274 * @return a Matrix 275 */ 276 public Matrix getDistMat(){ 277 return distanceMatrix; 278 } 279 280 /** get the parameters. 281 * 282 * @return the Parameters. 283 */ 284 public StrucAligParameters getParams() { 285 return params; 286 } 287 288 /** set the parameters to be used for the algorithm 289 * 290 * @param params the Parameter object 291 */ 292 public void setParams(StrucAligParameters params) { 293 this.params = params; 294 } 295 296 297 /** check if debug mode is set on 298 * 299 * @return debug flag 300 */ 301 @Deprecated 302 public boolean isDebug() { 303 return false; 304 } 305 306 /** set the debug flag 307 * 308 * @param debug flag 309 */ 310 @Deprecated 311 public void setDebug(boolean debug) { 312 } 313 314 /** Calculate the alignment between the two full structures with default parameters 315 * 316 * @param s1 317 * @param s2 318 * @throws StructureException 319 */ 320 public void align(Structure s1, Structure s2) 321 throws StructureException { 322 323 align(s1,s2,params); 324 } 325 326 /** Calculate the alignment between the two full structures with user provided parameters 327 * 328 * @param s1 329 * @param s2 330 * @param params 331 * @throws StructureException 332 */ 333 public void align(Structure s1, Structure s2, StrucAligParameters params) 334 throws StructureException { 335 // step 1 convert the structures to Atom Arrays 336 337 338 Atom[] ca1 = getAlignmentAtoms(s1); 339 Atom[] ca2 = getAlignmentAtoms(s2); 340 341 notifyStartingAlignment(s1.getName(),ca1,s2.getName(),ca2); 342 align(ca1,ca2,params); 343 } 344 345 346 /** Align two chains from the structures. Uses default parameters. 347 * 348 * @param s1 349 * @param chainId1 350 * @param s2 351 * @param chainId2 352 */ 353 public void align(Structure s1, String chainId1, Structure s2, String chainId2) 354 throws StructureException{ 355 align(s1,chainId1,s2,chainId2, params); 356 } 357 358 /** Aligns two chains from the structures using user provided parameters. 359 * 360 * @param s1 361 * @param chainId1 362 * @param s2 363 * @param chainId2 364 * @param params 365 * @throws StructureException 366 */ 367 public void align(Structure s1, String chainId1, Structure s2, String chainId2, StrucAligParameters params) 368 throws StructureException{ 369 reset(); 370 this.params = params; 371 372 Chain c1 = s1.getChainByPDB(chainId1); 373 Chain c2 = s2.getChainByPDB(chainId2); 374 375 Structure s3 = new StructureImpl(); 376 s3.addChain(c1); 377 378 Structure s4 = new StructureImpl(); 379 s4.addChain(c2); 380 381 Atom[] ca1 = getAlignmentAtoms(s3); 382 Atom[] ca2 = getAlignmentAtoms(s4); 383 384 notifyStartingAlignment(s1.getName(),ca1,s2.getName(),ca2); 385 align(ca1,ca2,params); 386 } 387 388 /** Returns the atoms that are being used for the alignment. (E.g. Calpha only, etc.) 389 * 390 * @param s 391 * @return an array of Atoms objects 392 */ 393 public Atom[] getAlignmentAtoms(Structure s){ 394 String[] atomNames = params.getUsedAtomNames(); 395 return StructureTools.getAtomArray(s,atomNames); 396 } 397 398 /** calculate the protein structure superimposition, between two sets of atoms. 399 * 400 * 401 * 402 * @param ca1 set of Atoms of structure 1 403 * @param ca2 set of Atoms of structure 2 404 * @param params the parameters to use for the alignment 405 * @throws StructureException 406 */ 407 public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params) 408 throws StructureException { 409 410 411 reset(); 412 this.params = params; 413 414 long timeStart = System.currentTimeMillis(); 415 416// step 1 get all Diagonals of length X that are similar between both structures 417 logger.debug(" length atoms1:" + ca1.length); 418 logger.debug(" length atoms2:" + ca2.length); 419 420 logger.debug("step 1 - get fragments with similar intramolecular distances "); 421 422 int k = params.getDiagonalDistance(); 423 int k2 = params.getDiagonalDistance2(); 424 int fragmentLength = params.getFragmentLength(); 425 426 if ( ca1.length < (fragmentLength + 1) ) { 427 throw new StructureException("structure 1 too short ("+ca1.length+"), can not align"); 428 } 429 if ( ca2.length < (fragmentLength + 1) ){ 430 throw new StructureException("structure 2 too short ("+ca2.length+"), can not align"); 431 } 432 int rows = ca1.length - fragmentLength + 1; 433 int cols = ca2.length - fragmentLength + 1; 434 distanceMatrix = new Matrix(rows,cols,0.0); 435 436 double[] dist1 = AlignTools.getDiagonalAtK(ca1, k); 437 438 double[] dist2 = AlignTools.getDiagonalAtK(ca2, k); 439 double[] dist3 = new double[0]; 440 double[] dist4 = new double[0]; 441 if ( k2 > 0) { 442 dist3 = AlignTools.getDiagonalAtK(ca1, k2); 443 dist4 = AlignTools.getDiagonalAtK(ca2, k2); 444 } 445 446 double[][] utmp = new double[][] {{0,0,1}}; 447 Atom unitvector = new AtomImpl(); 448 unitvector.setCoords(utmp[0]); 449 450 List<FragmentPair> fragments = new ArrayList<FragmentPair>(); 451 452 for ( int i = 0 ; i< rows; i++){ 453 454 Atom[] catmp1 = AlignTools.getFragment( ca1, i, fragmentLength); 455 Atom center1 = AlignTools.getCenter( ca1, i, fragmentLength); 456 457 for ( int j = 0 ; j < cols ; j++){ 458 459 double rdd1 = AlignTools.rms_dk_diag(dist1,dist2,i,j,fragmentLength,k); 460 double rdd2 = 0; 461 if ( k2 > 0) 462 rdd2 = AlignTools.rms_dk_diag(dist3,dist4,i,j,fragmentLength,k2); 463 double rdd = rdd1 + rdd2; 464 distanceMatrix.set(i,j,rdd); 465 466 467 if ( rdd < params.getFragmentMiniDistance()) { 468 FragmentPair f = new FragmentPair(fragmentLength,i,j); 469 try { 470 471 Atom[] catmp2 = AlignTools.getFragment(ca2, j, fragmentLength); 472 Atom center2 = AlignTools.getCenter(ca2,j,fragmentLength); 473 474 f.setCenter1(center1); 475 f.setCenter2(center2); 476 477 SVDSuperimposer svd = new SVDSuperimposer(catmp1,catmp2); 478 Matrix rotmat = svd.getRotation(); 479 f.setRot(rotmat); 480 481 Atom aunitv = (Atom)unitvector.clone(); 482 Calc.rotate(aunitv,rotmat); 483 f.setUnitv(aunitv); 484 485 boolean doNotAdd = false; 486 if ( params.reduceInitialFragments()) { 487 doNotAdd = FragmentJoiner.reduceFragments(fragments,f, distanceMatrix); 488 489 } 490 if ( doNotAdd) 491 continue; 492 493 fragments.add(f); 494 495 } catch (StructureException e){ 496 logger.error("Exception: ", e); 497 break; 498 } 499 } 500 } 501 } 502 503 notifyFragmentListeners(fragments); 504 505 FragmentPair[] fp = fragments.toArray(new FragmentPair[fragments.size()]); 506 setFragmentPairs(fp); 507 508 logger.debug(" got # fragment pairs: {}", fp.length); 509 510 logger.debug("step 2 - join fragments"); 511 512 // step 2 combine them to possible models 513 FragmentJoiner joiner = new FragmentJoiner(); 514 515 516 JointFragments[] frags; 517 518 if ( params.isJoinFast()) { 519 // apply the quick alignment procedure. 520 // less quality in alignments, better for DB searches... 521 frags = joiner.approach_ap3(ca1,ca2,fp,params); 522 523 joiner.extendFragments(ca1,ca2,frags,params); 524 525 } else if ( params.isJoinPlo()){ 526 // this approach by StrComPy (peter lackner): 527 frags = joiner.frag_pairwise_compat(fp, 528 params.getAngleDiff(), 529 params.getFragCompat(), 530 params.getMaxrefine()); 531 532 } else { 533 534 // my first implementation 535 frags = joiner.approach_ap3( 536 ca1,ca2, fp, params); 537 } 538 539 notifyJointFragments(frags); 540 541 logger.debug(" number joint fragments: ", frags.length); 542 543 logger.debug("step 3 - refine alignments"); 544 545 List<AlternativeAlignment> aas = new ArrayList<AlternativeAlignment>(); 546 for ( int i = 0 ; i < frags.length;i++){ 547 JointFragments f = frags[i]; 548 AlternativeAlignment a = new AlternativeAlignment(); 549 a.apairs_from_idxlst(f); 550 a.setAltAligNumber(i+1); 551 a.setDistanceMatrix(distanceMatrix); 552 553 try { 554 if ( params.getMaxIter() > 0 ){ 555 556 a.refine(params,ca1,ca2); 557 } 558 else { 559 560 a.finish(params,ca1,ca2); 561 562 } 563 } catch (StructureException e){ 564 logger.error("Refinement of fragment {} failed", i, e); 565 } 566 a.calcScores(ca1,ca2); 567 aas.add(a); 568 } 569 570 571 // sort the alternative alignments 572 Comparator<AlternativeAlignment> comp = new AltAligComparator(); 573 Collections.sort(aas,comp); 574 Collections.reverse(aas); 575 576 alts = aas.toArray(new AlternativeAlignment[aas.size()]); 577 // do final numbering of alternative solutions 578 int aanbr = 0; 579 for (AlternativeAlignment a : alts) { 580 aanbr++; 581 a.setAltAligNumber(aanbr); 582 } 583 584 logger.debug("total calculation time: {} ms.", (System.currentTimeMillis()-timeStart)); 585 } 586 587 private void notifyStartingAlignment(String name1, Atom[] ca1, String name2, Atom[] ca2){ 588 for (AlignmentProgressListener li : listeners){ 589 li.startingAlignment(name1, ca1, name2, ca2); 590 } 591 } 592 593 private void notifyFragmentListeners(List<FragmentPair> fragments){ 594 595 for (AlignmentProgressListener li : listeners){ 596 li.calculatedFragmentPairs(fragments); 597 } 598 599 } 600 601 private void notifyJointFragments(JointFragments[] fragments){ 602 for (AlignmentProgressListener li : listeners){ 603 li.jointFragments(fragments); 604 } 605 } 606 607}