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