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 */ 021package org.biojava.nbio.structure.symmetry.utils; 022 023import java.util.ArrayList; 024import java.util.Arrays; 025import java.util.Collections; 026import java.util.List; 027 028import javax.vecmath.AxisAngle4d; 029import javax.vecmath.Matrix4d; 030import javax.vecmath.Point3d; 031 032import org.biojava.nbio.structure.Atom; 033import org.biojava.nbio.structure.Calc; 034import org.biojava.nbio.structure.Chain; 035import org.biojava.nbio.structure.ChainImpl; 036import org.biojava.nbio.structure.Group; 037import org.biojava.nbio.structure.SVDSuperimposer; 038import org.biojava.nbio.structure.Structure; 039import org.biojava.nbio.structure.StructureException; 040import org.biojava.nbio.structure.StructureIdentifier; 041import org.biojava.nbio.structure.StructureImpl; 042import org.biojava.nbio.structure.StructureTools; 043import org.biojava.nbio.structure.align.ce.CECalculator; 044import org.biojava.nbio.structure.align.helper.AlignTools; 045import org.biojava.nbio.structure.align.model.AFPChain; 046import org.biojava.nbio.structure.align.multiple.Block; 047import org.biojava.nbio.structure.align.multiple.BlockImpl; 048import org.biojava.nbio.structure.align.multiple.BlockSet; 049import org.biojava.nbio.structure.align.multiple.BlockSetImpl; 050import org.biojava.nbio.structure.align.multiple.MultipleAlignment; 051import org.biojava.nbio.structure.align.multiple.MultipleAlignmentEnsemble; 052import org.biojava.nbio.structure.align.multiple.MultipleAlignmentEnsembleImpl; 053import org.biojava.nbio.structure.align.multiple.MultipleAlignmentImpl; 054import org.biojava.nbio.structure.align.multiple.util.CoreSuperimposer; 055import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentScorer; 056import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentTools; 057import org.biojava.nbio.structure.align.multiple.util.MultipleSuperimposer; 058import org.biojava.nbio.structure.jama.Matrix; 059import org.biojava.nbio.structure.symmetry.core.QuatSymmetryDetector; 060import org.biojava.nbio.structure.symmetry.core.QuatSymmetryParameters; 061import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults; 062import org.biojava.nbio.structure.symmetry.core.Subunits; 063import org.biojava.nbio.structure.symmetry.internal.CeSymmResult; 064import org.biojava.nbio.structure.symmetry.internal.SymmetryAxes; 065import org.jgrapht.UndirectedGraph; 066import org.jgrapht.graph.DefaultEdge; 067import org.jgrapht.graph.SimpleGraph; 068import org.slf4j.Logger; 069import org.slf4j.LoggerFactory; 070 071/** 072 * Utility methods for the internal symmetry identification and manipulation. 073 * <p> 074 * Methods include: blank out regions of DP Matrix, build symmetry graphs, get 075 * rotation symmetry angles, split repeats in quaternary structure chains, 076 * convert between symmetry formats (full, repeats, rotations), determine if two 077 * symmetry axes are equivalent, get groups from representative Atoms. 078 * 079 * @author Spencer Bliven 080 * @author Aleix Lafita 081 * 082 */ 083public class SymmetryTools { 084 085 private static final Logger logger = LoggerFactory 086 .getLogger(SymmetryTools.class); 087 088 /** Prevent instantiation. */ 089 private SymmetryTools() { 090 } 091 092 /** 093 * Returns the "reset value" for graying out the main diagonal. If we're 094 * blanking out the main diagonal, this value is always Integer.MIN_VALUE. 095 * <p> 096 * This is possible if {@code gradientPolyCoeff = Integer.MIN_VALUE} and 097 * {@code gradientExpCoeff = 0}. 098 * 099 * @param unpenalizedScore 100 * @param nResFromMainDiag 101 * @param gradientPolyCoeff 102 * @param gradientExpCoeff 103 * @return 104 */ 105 private static double getResetVal(double unpenalizedScore, 106 double nResFromMainDiag, double[] gradientPolyCoeff, 107 double gradientExpCoeff) { 108 109 if (unpenalizedScore == Double.NaN) 110 return 0; // what else? 111 112 // We can actually return a positive value if this is high enough 113 double updateVal = unpenalizedScore; 114 updateVal -= gradientExpCoeff * Math.pow(Math.E, -nResFromMainDiag); 115 for (int p = 0; p < gradientPolyCoeff.length; p++) { 116 updateVal -= gradientPolyCoeff[gradientPolyCoeff.length - 1 - p] 117 * Math.pow(nResFromMainDiag, -p); 118 } 119 return updateVal; 120 } 121 122 /** 123 * Grays out the main diagonal of a duplicated distance matrix. 124 * 125 * @param ca2 126 * @param rows 127 * Number of rows 128 * @param cols 129 * Number of original columns 130 * @param calculator 131 * Used to get the matrix if origM is null 132 * @param origM 133 * starting matrix. If null, uses 134 * {@link CECalculator#getMatMatrix()} 135 * @param blankWindowSize 136 * Width of section to gray out 137 * @param gradientPolyCoeff 138 * @param gradientExpCoeff 139 * @return 140 */ 141 public static Matrix grayOutCEOrig(Atom[] ca2, int rows, int cols, 142 CECalculator calculator, Matrix origM, int blankWindowSize, 143 double[] gradientPolyCoeff, double gradientExpCoeff) { 144 145 if (origM == null) { 146 origM = new Matrix(calculator.getMatMatrix()); 147 } 148 149 // symmetry hack, disable main diagonal 150 151 for (int i = 0; i < rows; i++) { 152 for (int j = 0; j < cols; j++) { 153 int diff = Math.abs(i - j); 154 155 double resetVal = getResetVal(origM.get(i, j), diff, 156 gradientPolyCoeff, gradientExpCoeff); 157 158 if (diff < blankWindowSize) { 159 origM.set(i, j, origM.get(i, j) + resetVal); 160 161 } 162 int diff2 = Math.abs(i - (j - ca2.length / 2)); // other side 163 164 double resetVal2 = getResetVal(origM.get(i, j), diff2, 165 gradientPolyCoeff, gradientExpCoeff); 166 167 if (diff2 < blankWindowSize) { 168 origM.set(i, j, origM.get(i, j) + resetVal2); 169 170 } 171 } 172 } 173 return origM; 174 } 175 176 public static Matrix grayOutPreviousAlignment(AFPChain afpChain, 177 Atom[] ca2, int rows, int cols, CECalculator calculator, 178 Matrix max, int blankWindowSize, double[] gradientPolyCoeff, 179 double gradientExpCoeff) { 180 181 max = grayOutCEOrig(ca2, rows, cols, calculator, max, blankWindowSize, 182 gradientPolyCoeff, gradientExpCoeff); 183 184 double[][] dist1 = calculator.getDist1(); 185 double[][] dist2 = calculator.getDist2(); 186 187 int[][][] optAln = afpChain.getOptAln(); 188 int blockNum = afpChain.getBlockNum(); 189 190 int[] optLen = afpChain.getOptLen(); 191 192 // ca2 is circularly permutated 193 int breakPoint = ca2.length / 2; 194 for (int bk = 0; bk < blockNum; bk++) { 195 196 for (int i = 0; i < optLen[bk]; i++) { 197 int pos1 = optAln[bk][0][i]; 198 int pos2 = optAln[bk][1][i]; 199 200 int dist = blankWindowSize / 2; 201 int start1 = Math.max(pos1 - dist, 0); 202 int start2 = Math.max(pos2 - dist, 0); 203 int end1 = Math.min(pos1 + dist, rows - 1); 204 int end2 = Math.min(pos2 + dist, cols - 1); 205 206 for (int i1 = start1; i1 < end1; i1++) { 207 208 // blank diagonal of dist1 209 for (int k = 0; k < blankWindowSize / 2; k++) { 210 if (i1 - k >= 0) { 211 double resetVal = getResetVal( 212 max.get(i1 - k, i1 - k), 0, 213 gradientPolyCoeff, gradientExpCoeff); 214 dist1[i1 - k][i1 - k] = resetVal; 215 } else if (i1 + k < rows) { 216 double resetVal = getResetVal( 217 max.get(i1 + k, i1 + k), 0, 218 gradientPolyCoeff, gradientExpCoeff); 219 dist1[i1 + k][i1 + k] = resetVal; 220 } 221 222 } 223 224 for (int j2 = start2; j2 < end2; j2++) { 225 double resetVal = getResetVal(max.get(i1, j2), 226 Math.abs(i1 - j2), gradientPolyCoeff, 227 gradientExpCoeff); 228 max.set(i1, j2, resetVal); 229 if (j2 < breakPoint) { 230 double resetVal2 = getResetVal( 231 max.get(i1, j2 + breakPoint), 232 Math.abs(i1 - (j2 + breakPoint)), 233 gradientPolyCoeff, gradientExpCoeff); 234 max.set(i1, j2 + breakPoint, resetVal2); 235 } else { 236 double resetVal2 = getResetVal( 237 max.get(i1, j2 - breakPoint), 238 Math.abs(i1 - (j2 - breakPoint)), 239 gradientPolyCoeff, gradientExpCoeff); 240 max.set(i1, j2 - breakPoint, resetVal2); 241 } 242 for (int k = 0; k < blankWindowSize / 2; k++) { 243 if (j2 - k >= 0) { 244 if (j2 - k < breakPoint) { 245 double resetVal2 = getResetVal( 246 max.get(j2 - k, j2 - k), 0, 247 gradientPolyCoeff, gradientExpCoeff); 248 dist2[j2 - k][j2 - k] = resetVal2; 249 } else { 250 double resetVal2 = getResetVal(max.get(j2 251 - k - breakPoint, j2 - k), 0, 252 gradientPolyCoeff, gradientExpCoeff); 253 dist2[j2 - k - breakPoint][j2 - k 254 - breakPoint] = resetVal2; 255 } 256 } else if (j2 + k < cols) { 257 if (j2 + k < breakPoint) { 258 double resetVal2 = getResetVal( 259 max.get(j2 + k, j2 + k), 0, 260 gradientPolyCoeff, gradientExpCoeff); 261 dist2[j2 + k][j2 + k] = resetVal2; 262 } else { 263 double resetVal2 = getResetVal(max.get(j2 264 + k - breakPoint, j2 + k), 0, 265 gradientPolyCoeff, gradientExpCoeff); 266 dist2[j2 + k - breakPoint][j2 + k 267 - breakPoint] = resetVal2; 268 } 269 } 270 } 271 } 272 } 273 274 } 275 } 276 calculator.setDist1(dist1); 277 calculator.setDist2(dist2); 278 return max; 279 280 } 281 282 public Matrix getDkMatrix(Atom[] ca1, Atom[] ca2, int fragmentLength, 283 double[] dist1, double[] dist2, int rows, int cols) { 284 285 Matrix diffDistMax = Matrix.identity(ca1.length, ca2.length); 286 287 for (int i = 0; i < rows; i++) { 288 double score1 = 0; 289 for (int x = 0; x < fragmentLength; x++) { 290 score1 += dist1[i + x]; 291 } 292 for (int j = 0; j < cols; j++) { 293 double score2 = 0; 294 for (int y = 0; y < fragmentLength; y++) { 295 score2 += dist2[j + y]; 296 } 297 298 // if the intramolecular distances are very similar 299 // the two scores should be similar, 300 // i.e. the difference is close to 0 301 diffDistMax.set(i, j, Math.abs(score1 - score2)); 302 } 303 } 304 305 // symmetry hack, disable main diagonal 306 307 for (int i = 0; i < rows; i++) { 308 for (int j = 0; j < cols; j++) { 309 int diff = Math.abs(i - j); 310 311 if (diff < 15) { 312 diffDistMax.set(i, j, 99); 313 } 314 int diff2 = Math.abs(i - (j - ca2.length / 2)); 315 if (diff2 < 15) { 316 diffDistMax.set(i, j, 99); 317 } 318 } 319 } 320 return diffDistMax; 321 322 } 323 324 public static Matrix blankOutPreviousAlignment(AFPChain afpChain, 325 Atom[] ca2, int rows, int cols, CECalculator calculator, 326 Matrix max, int blankWindowSize) { 327 return grayOutPreviousAlignment(afpChain, ca2, rows, cols, calculator, 328 max, blankWindowSize, new double[] { Integer.MIN_VALUE }, 0.0); 329 } 330 331 public static Matrix blankOutCEOrig(Atom[] ca2, int rows, int cols, 332 CECalculator calculator, Matrix origM, int blankWindowSize) { 333 return grayOutCEOrig(ca2, rows, cols, calculator, origM, 334 blankWindowSize, new double[] { Integer.MIN_VALUE }, 0.0); 335 } 336 337 public static Matrix getDkMatrix(Atom[] ca1, Atom[] ca2, int k, 338 int fragmentLength) { 339 340 double[] dist1 = AlignTools.getDiagonalAtK(ca1, k); 341 double[] dist2 = AlignTools.getDiagonalAtK(ca2, k); 342 343 int rows = ca1.length - fragmentLength - k + 1; 344 int cols = ca2.length - fragmentLength - k + 1; 345 346 // Matrix that tracks similarity of a fragment of length fragmentLength 347 // starting a position i,j. 348 349 Matrix m2 = new Matrix(rows, cols); 350 351 for (int i = 0; i < rows; i++) { 352 double score1 = 0; 353 for (int x = 0; x < fragmentLength; x++) { 354 score1 += dist1[i + x]; 355 } 356 for (int j = 0; j < cols; j++) { 357 double score2 = 0; 358 for (int y = 0; y < fragmentLength; y++) { 359 score2 += dist2[j + y]; 360 } 361 362 // if the intramolecular distances are very similar 363 // the two scores should be similar, 364 // i.e. the difference is close to 0 365 m2.set(i, j, Math.abs(score1 - score2)); 366 } 367 } 368 return m2; 369 } 370 371 public static boolean[][] blankOutBreakFlag(AFPChain afpChain, Atom[] ca2, 372 int rows, int cols, CECalculator calculator, boolean[][] breakFlag, 373 int blankWindowSize) { 374 375 int[][][] optAln = afpChain.getOptAln(); 376 int blockNum = afpChain.getBlockNum(); 377 378 int[] optLen = afpChain.getOptLen(); 379 380 // ca2 is circularly permutated at this point. 381 int breakPoint = ca2.length / 2; 382 383 for (int bk = 0; bk < blockNum; bk++) { 384 385 // Matrix m= afpChain.getBlockRotationMatrix()[bk]; 386 // Atom shift = afpChain.getBlockShiftVector()[bk]; 387 for (int i = 0; i < optLen[bk]; i++) { 388 int pos1 = optAln[bk][0][i]; 389 int pos2 = optAln[bk][1][i]; 390 // blank out area around these positions... 391 392 int dist = blankWindowSize; 393 int start1 = Math.max(pos1 - dist, 0); 394 int start2 = Math.max(pos2 - dist, 0); 395 int end1 = Math.min(pos1 + dist, rows - 1); 396 int end2 = Math.min(pos2 + dist, cols - 1); 397 398 for (int i1 = start1; i1 < end1; i1++) { 399 400 for (int j2 = start2; j2 < end2; j2++) { 401 // System.out.println(i1 + " " + j2 + " (***)"); 402 breakFlag[i1][j2] = true; 403 if (j2 < breakPoint) { 404 breakFlag[i1][j2 + breakPoint] = true; 405 } 406 } 407 } 408 409 } 410 } 411 412 return breakFlag; 413 } 414 415 /** 416 * Returns the <em>magnitude</em> of the angle between the first and second 417 * blocks of {@code afpChain}, measured in degrees. This is always a 418 * positive value (unsigned). 419 * 420 * @param afpChain 421 * @param ca1 422 * @param ca2 423 * @return 424 */ 425 public static double getAngle(AFPChain afpChain, Atom[] ca1, Atom[] ca2) { 426 Matrix rotation = afpChain.getBlockRotationMatrix()[0]; 427 return Math.acos(rotation.trace() - 1) * 180 / Math.PI; 428 } 429 430 /** 431 * Converts a set of AFP alignments into a Graph of aligned residues, where 432 * each vertex is a residue and each edge means the connection between the 433 * two residues in one of the alignments. 434 * 435 * @param afps 436 * List of AFPChains 437 * @param atoms 438 * Atom array of the symmetric structure 439 * @param undirected 440 * if true, the graph is undirected 441 * 442 * @return adjacency List of aligned residues 443 */ 444 public static List<List<Integer>> buildSymmetryGraph(List<AFPChain> afps, 445 Atom[] atoms, boolean undirected) { 446 447 List<List<Integer>> graph = new ArrayList<List<Integer>>(); 448 449 for (int n = 0; n < atoms.length; n++) { 450 graph.add(new ArrayList<Integer>()); 451 } 452 453 for (int k = 0; k < afps.size(); k++) { 454 for (int i = 0; i < afps.get(k).getOptAln().length; i++) { 455 for (int j = 0; j < afps.get(k).getOptAln()[i][0].length; j++) { 456 Integer res1 = afps.get(k).getOptAln()[i][0][j]; 457 Integer res2 = afps.get(k).getOptAln()[i][1][j]; 458 graph.get(res1).add(res2); 459 if (undirected) 460 graph.get(res2).add(res1); 461 } 462 } 463 } 464 return graph; 465 } 466 467 /** 468 * Converts a self alignment into a directed jGraphT of aligned residues, 469 * where each vertex is a residue and each edge means the equivalence 470 * between the two residues in the self-alignment. 471 * 472 * @param selfAlignment 473 * AFPChain 474 * 475 * @return alignment Graph 476 */ 477 public static UndirectedGraph<Integer, DefaultEdge> buildSymmetryGraph( 478 AFPChain selfAlignment) { 479 480 UndirectedGraph<Integer, DefaultEdge> graph = new SimpleGraph<Integer, DefaultEdge>( 481 DefaultEdge.class); 482 483 for (int i = 0; i < selfAlignment.getOptAln().length; i++) { 484 for (int j = 0; j < selfAlignment.getOptAln()[i][0].length; j++) { 485 Integer res1 = selfAlignment.getOptAln()[i][0][j]; 486 Integer res2 = selfAlignment.getOptAln()[i][1][j]; 487 graph.addVertex(res1); 488 graph.addVertex(res2); 489 graph.addEdge(res1, res2); 490 } 491 } 492 return graph; 493 } 494 495 /** 496 * Method that converts the symmetric units of a structure into different 497 * structures, so that they can be individually visualized. 498 * 499 * @param symmetry 500 * CeSymmResult 501 * @throws StructureException 502 * @result List of structures, by repeat index sequentially 503 * 504 */ 505 public static List<Structure> divideStructure(CeSymmResult symmetry) 506 throws StructureException { 507 508 if (!symmetry.isRefined()) 509 throw new IllegalArgumentException("The symmetry result " 510 + "is not refined, repeats cannot be defined"); 511 512 int order = symmetry.getMultipleAlignment().size(); 513 Atom[] atoms = StructureTools.cloneAtomArray(symmetry.getAtoms()); 514 List<StructureIdentifier> repeatsId = symmetry.getRepeatsID(); 515 List<Structure> repeats = new ArrayList<Structure>(order); 516 517 // Create new structure containing the repeat atoms 518 for (int i = 0; i < order; i++) { 519 520 Structure s = new StructureImpl(); 521 s.setStructureIdentifier(repeatsId.get(i)); 522 523 Block align = symmetry.getMultipleAlignment().getBlock(0); 524 525 // Get the start and end of the repeat 526 int res1 = align.getStartResidue(i); 527 int res2 = align.getFinalResidue(i); 528 529 Atom[] repeat = Arrays.copyOfRange(atoms, res1, res2 + 1); 530 531 Chain newCh = new ChainImpl(); 532 newCh.setChainID(repeat[0].getGroup().getChainId()); 533 534 for (int k = 0; k < repeat.length; k++) { 535 Group g = (Group) repeat[k].getGroup().clone(); 536 newCh.addGroup(g); 537 } 538 s.addChain(newCh); 539 repeats.add(s); 540 } 541 return repeats; 542 } 543 544 /** 545 * Method that converts a repeats symmetric alignment into an alignment of 546 * whole structures. 547 * <p> 548 * Example: if the structure has repeats A,B and C, the original alignment 549 * is A-B-C, and the returned alignment is ABC-BCA-CAB. 550 * 551 * @param symm 552 * CeSymmResult 553 * @return MultipleAlignment of the full structure superpositions 554 */ 555 public static MultipleAlignment toFullAlignment(CeSymmResult symm) { 556 557 if (!symm.isRefined()) 558 throw new IllegalArgumentException("The symmetry result " 559 + "is not refined, repeats cannot be defined"); 560 561 MultipleAlignment full = symm.getMultipleAlignment().clone(); 562 563 for (int str = 1; str < full.size(); str++) { 564 // Create a new Block with swapped AlignRes (move first to last) 565 Block b = full.getBlock(full.getBlocks().size() - 1).clone(); 566 b.getAlignRes().add(b.getAlignRes().get(0)); 567 b.getAlignRes().remove(0); 568 full.getBlockSet(0).getBlocks().add(b); 569 } 570 return full; 571 } 572 573 /** 574 * Method that converts a symmetry alignment into an alignment of the 575 * repeats only, as new independent structures. 576 * <p> 577 * This method changes the structure identifiers, the Atom arrays and 578 * re-scles the aligned residues in the Blocks corresponding to those 579 * changes. 580 * <p> 581 * Application: display superimposed repeats in Jmol. 582 * 583 * @param result 584 * CeSymmResult of symmetry 585 * @return MultipleAlignment of the repeats 586 * @throws StructureException 587 */ 588 public static MultipleAlignment toRepeatsAlignment(CeSymmResult result) 589 throws StructureException { 590 591 if (!result.isRefined()) 592 throw new IllegalArgumentException("The symmetry result " 593 + "is not refined, repeats cannot be defined"); 594 595 MultipleAlignment msa = result.getMultipleAlignment(); 596 MultipleAlignmentEnsemble newEnsemble = msa.getEnsemble().clone(); 597 598 List<Structure> repSt = SymmetryTools.divideStructure(result); 599 600 MultipleAlignment repeats = newEnsemble.getMultipleAlignment(0); 601 Block block = repeats.getBlock(0); 602 List<Atom[]> atomArrays = new ArrayList<Atom[]>(); 603 604 for (Structure s : repSt) 605 atomArrays.add(StructureTools.getRepresentativeAtomArray(s)); 606 607 newEnsemble.setAtomArrays(atomArrays); 608 609 for (int su = 0; su < block.size(); su++) { 610 Integer start = block.getStartResidue(su); 611 612 // Normalize aligned residues 613 for (int res = 0; res < block.length(); res++) { 614 Integer residue = block.getAlignRes().get(su).get(res); 615 if (residue != null) 616 residue -= start; 617 block.getAlignRes().get(su).set(res, residue); 618 } 619 } 620 621 return repeats; 622 } 623 624 /** 625 * Converts a refined symmetry AFPChain alignment into the standard 626 * representation of symmetry in a MultipleAlignment, that contains the 627 * entire Atom array of the strcuture and the symmetric repeats are orgaized 628 * in different rows in a single Block. 629 * 630 * @param symm 631 * AFPChain created with a symmetry algorithm and refined 632 * @param atoms 633 * Atom array of the entire structure 634 * @return MultipleAlignment format of the symmetry 635 * @throws StructureException 636 */ 637 public static MultipleAlignment fromAFP(AFPChain symm, Atom[] atoms) 638 throws StructureException { 639 640 if (!symm.getAlgorithmName().contains("symm")) { 641 throw new IllegalArgumentException( 642 "The input alignment is not a symmetry alignment."); 643 } 644 645 MultipleAlignmentEnsemble e = new MultipleAlignmentEnsembleImpl(symm, 646 atoms, atoms, false); 647 e.setAtomArrays(new ArrayList<Atom[]>()); 648 StructureIdentifier name = null; 649 if (e.getStructureIdentifiers() != null) { 650 if (!e.getStructureIdentifiers().isEmpty()) 651 name = e.getStructureIdentifiers().get(0); 652 } else 653 name = atoms[0].getGroup().getChain().getStructure() 654 .getStructureIdentifier(); 655 656 e.setStructureIdentifiers(new ArrayList<StructureIdentifier>()); 657 658 MultipleAlignment result = new MultipleAlignmentImpl(); 659 BlockSet bs = new BlockSetImpl(result); 660 Block b = new BlockImpl(bs); 661 b.setAlignRes(new ArrayList<List<Integer>>()); 662 663 int order = symm.getBlockNum(); 664 for (int su = 0; su < order; su++) { 665 List<Integer> residues = e.getMultipleAlignment(0).getBlock(su) 666 .getAlignRes().get(0); 667 b.getAlignRes().add(residues); 668 e.getStructureIdentifiers().add(name); 669 e.getAtomArrays().add(atoms); 670 } 671 e.getMultipleAlignments().set(0, result); 672 result.setEnsemble(e); 673 674 CoreSuperimposer imposer = new CoreSuperimposer(); 675 imposer.superimpose(result); 676 updateSymmetryScores(result); 677 678 return result; 679 } 680 681 /** 682 * Determines if two symmetry axis are equivalent inside the error 683 * threshold. It only takes into account the direction of the vector where 684 * the rotation is made: the angle and translation are not taken into 685 * account. 686 * 687 * @param axis1 688 * @param axis2 689 * @param epsilon 690 * error allowed in the axis comparison 691 * @return true if equivalent, false otherwise 692 */ 693 @Deprecated 694 public static boolean equivalentAxes(Matrix4d axis1, Matrix4d axis2, 695 double epsilon) { 696 697 AxisAngle4d rot1 = new AxisAngle4d(); 698 rot1.set(axis1); 699 AxisAngle4d rot2 = new AxisAngle4d(); 700 rot2.set(axis2); 701 702 // rot1.epsilonEquals(rot2, error); //that also compares angle 703 // L-infinite distance without comparing the angle (epsilonEquals) 704 List<Double> sameDir = new ArrayList<Double>(); 705 sameDir.add(Math.abs(rot1.x - rot2.x)); 706 sameDir.add(Math.abs(rot1.y - rot2.y)); 707 sameDir.add(Math.abs(rot1.z - rot2.z)); 708 709 List<Double> otherDir = new ArrayList<Double>(); 710 otherDir.add(Math.abs(rot1.x + rot2.x)); 711 otherDir.add(Math.abs(rot1.y + rot2.y)); 712 otherDir.add(Math.abs(rot1.z + rot2.z)); 713 714 Double error = Math.min(Collections.max(sameDir), 715 Collections.max(otherDir)); 716 717 return error < epsilon; 718 } 719 720 /** 721 * Given a symmetry result, it calculates the overall global symmetry, 722 * factoring out the alignment and detection steps of 723 * {@link QuatSymmetryDetector} algorithm. 724 * 725 * @param result 726 * symmetry result 727 * @return global symmetry results 728 * @throws StructureException 729 */ 730 public static QuatSymmetryResults getQuaternarySymmetry(CeSymmResult result) 731 throws StructureException { 732 733 // Obtain the clusters of aligned Atoms and repeat variables 734 List<Atom[]> alignedCA = result.getMultipleAlignment().getAtomArrays(); 735 MultipleAlignment msa = result.getMultipleAlignment(); 736 List<Integer> corePos = MultipleAlignmentTools.getCorePositions(result 737 .getMultipleAlignment().getBlock(0)); 738 739 List<Point3d[]> caCoords = new ArrayList<Point3d[]>(); 740 List<Integer> folds = new ArrayList<Integer>(); 741 List<Boolean> pseudo = new ArrayList<Boolean>(); 742 List<String> chainIds = new ArrayList<String>(); 743 List<Integer> models = new ArrayList<Integer>(); 744 List<Double> seqIDmin = new ArrayList<Double>(); 745 List<Double> seqIDmax = new ArrayList<Double>(); 746 List<Integer> clusterIDs = new ArrayList<Integer>(); 747 int fold = 1; 748 749 for (int str = 0; str < alignedCA.size(); str++) { 750 Atom[] array = alignedCA.get(str); 751 List<Point3d> points = new ArrayList<Point3d>(); 752 List<Integer> alignedRes = msa.getBlock(0).getAlignRes().get(str); 753 for (int pos = 0; pos < alignedRes.size(); pos++) { 754 Integer residue = alignedRes.get(pos); 755 if (residue == null) 756 continue; 757 else if (!corePos.contains(pos)) 758 continue; 759 Atom a = array[residue]; 760 points.add(new Point3d(a.getCoords())); 761 } 762 caCoords.add(points.toArray(new Point3d[points.size()])); 763 if (alignedCA.size() % fold == 0) { 764 folds.add(fold); // the folds are the common denominators 765 } 766 fold++; 767 pseudo.add(false); 768 chainIds.add(alignedCA.get(str)[0].getGroup().getChainId()); 769 models.add(0); 770 seqIDmax.add(1.0); 771 seqIDmin.add(1.0); 772 clusterIDs.add(0); 773 } 774 775 // Create directly the repeats, because we know the aligned CA 776 Subunits globalSubunits = new Subunits(caCoords, clusterIDs, pseudo, 777 seqIDmin, seqIDmax, folds, chainIds, models); 778 779 // Quaternary Symmetry Detection 780 QuatSymmetryParameters param = new QuatSymmetryParameters(); 781 782 QuatSymmetryResults gSymmetry = QuatSymmetryDetector.calcQuatSymmetry( 783 globalSubunits, param); 784 785 return gSymmetry; 786 } 787 788 /** 789 * Returns true a symmetry multiple alignment has been refined, false 790 * otherwise. 791 * <p> 792 * For a refined alignment only one Block with no repeated residues is 793 * necessary. Sufficient condition is not tested (only known from the 794 * algorithm or CeSymmResult). 795 * 796 * @param symm 797 * the symmetry alignment 798 * @return true if the alignment is refined 799 */ 800 @Deprecated 801 public static boolean isRefined(MultipleAlignment symm) { 802 803 if (symm.getBlocks().size() > 1) { 804 return false; 805 } else if (symm.size() < 2) 806 return false; 807 else { 808 List<Integer> alreadySeen = new ArrayList<Integer>(); 809 List<List<Integer>> align = symm.getBlock(0).getAlignRes(); 810 for (int str = 0; str < symm.size(); str++) { 811 for (int res = 0; res < align.get(str).size(); res++) { 812 Integer residue = align.get(str).get(res); 813 if (residue == null) 814 continue; 815 if (alreadySeen.contains(residue)) { 816 return false; 817 } else { 818 alreadySeen.add(residue); 819 } 820 } 821 } // end of all repeats 822 return true; 823 } 824 } 825 826 /** 827 * Returns true if the symmetry alignment is significant, false otherwise. 828 * <p> 829 * For a symmetry alignment to be significant, the alignment has to be 830 * refined and the TM-score has to be higher than the threshold. 831 * <p> 832 * It is recommended to use the {@link CeSymmResult#isSignificant()} method 833 * instead. 834 * 835 * @param msa 836 * @param symmetryThreshold 837 * @return 838 * @throws StructureException 839 */ 840 @Deprecated 841 public static boolean isSignificant(MultipleAlignment msa, 842 double symmetryThreshold) throws StructureException { 843 844 // Order/refinement check 845 if (!SymmetryTools.isRefined(msa)) 846 return false; 847 848 // TM-score cutoff 849 double tm = 0.0; 850 if (msa.getScore(MultipleAlignmentScorer.AVGTM_SCORE) == null) 851 tm = MultipleAlignmentScorer.getAvgTMScore(msa); 852 else 853 tm = msa.getScore(MultipleAlignmentScorer.AVGTM_SCORE); 854 855 if (tm < symmetryThreshold) 856 return false; 857 858 return true; 859 } 860 861 /** 862 * Returns the List of Groups of the corresponding representative Atom 863 * array. The representative Atom array needs to fulfill: no two Atoms are 864 * from the same Group and Groups are sequential (connected in the original 865 * Structure), except if they are from different Chains. 866 * 867 * @param rAtoms 868 * array of representative Atoms (CA, P, etc). 869 * @return List of Groups 870 */ 871 public static List<Group> getGroups(Atom[] rAtoms) { 872 873 List<Group> groups = new ArrayList<Group>(rAtoms.length); 874 875 for (Atom a : rAtoms) { 876 Group g = a.getGroup(); 877 if (g != null) 878 groups.add(g); 879 else 880 logger.info("Group not found for representative Atom {}", a); 881 } 882 return groups; 883 } 884 885 /** 886 * Calculates the set of symmetry operation Matrices (transformations) of 887 * the new alignment, based on the symmetry relations in the SymmetryAxes 888 * object. It sets the transformations to the input MultipleAlignment and 889 * SymmetryAxes objects. If the SymmetryAxes object is null, the 890 * superposition of the repeats is done without symmetry constraints. 891 * <p> 892 * This method also sets the scores (RMSD and TM-score) after the new 893 * superposition has been updated. 894 * 895 * @param axes 896 * SymmetryAxes object. It will be modified. 897 * @param msa 898 * MultipleAlignment. It will be modified. 899 */ 900 public static void updateSymmetryTransformation(SymmetryAxes axes, 901 MultipleAlignment msa) throws StructureException { 902 903 List<List<Integer>> block = msa.getBlocks().get(0).getAlignRes(); 904 int length = block.get(0).size(); 905 906 if (axes != null) { 907 for (int level = 0; level < axes.getNumLevels(); level++) { 908 909 // Calculate the aligned atom arrays to superimpose 910 List<Atom> list1 = new ArrayList<Atom>(); 911 List<Atom> list2 = new ArrayList<Atom>(); 912 913 for (int firstRepeat : axes.getFirstRepeats(level)) { 914 915 Matrix4d transform = axes.getRepeatTransform(firstRepeat); 916 917 List<List<Integer>> relation = axes.getRepeatRelation( 918 level, firstRepeat); 919 920 for (int index = 0; index < relation.get(0).size(); index++) { 921 int p1 = relation.get(0).get(index); 922 int p2 = relation.get(1).get(index); 923 924 for (int k = 0; k < length; k++) { 925 Integer pos1 = block.get(p1).get(k); 926 Integer pos2 = block.get(p2).get(k); 927 if (pos1 != null && pos2 != null) { 928 Atom a = (Atom) msa.getAtomArrays().get(p1)[pos1].clone(); 929 Atom b = (Atom) msa.getAtomArrays().get(p2)[pos2].clone(); 930 Calc.transform(a, transform); 931 Calc.transform(b, transform); 932 list1.add(a); 933 list2.add(b); 934 } 935 } 936 } 937 } 938 939 Atom[] arr1 = list1.toArray(new Atom[list1.size()]); 940 Atom[] arr2 = list2.toArray(new Atom[list2.size()]); 941 942 // Calculate the new transformation information 943 if (arr1.length > 0 && arr2.length > 0) { 944 SVDSuperimposer svd = new SVDSuperimposer(arr1, arr2); 945 Matrix4d axis = svd.getTransformation(); 946 axes.updateAxis(level, axis); 947 } 948 949 // Get the transformations from the SymmetryAxes 950 List<Matrix4d> transformations = new ArrayList<Matrix4d>(); 951 for (int su = 0; su < msa.size(); su++) { 952 transformations.add(axes.getRepeatTransform(su)); 953 } 954 msa.getBlockSet(0).setTransformations(transformations); 955 } 956 } else { 957 MultipleSuperimposer imposer = new CoreSuperimposer(); 958 imposer.superimpose(msa); 959 } 960 updateSymmetryScores(msa); 961 } 962 963 /** 964 * Update the scores (TM-score and RMSD) of a symmetry multiple alignment. 965 * This method does not redo the superposition of the alignment. 966 * 967 * @param symm 968 * Symmetry Multiple Alignment of Repeats 969 * @throws StructureException 970 */ 971 public static void updateSymmetryScores(MultipleAlignment symm) 972 throws StructureException { 973 974 // Multiply by the order of symmetry to normalize score 975 double tmScore = MultipleAlignmentScorer.getAvgTMScore(symm) 976 * symm.size(); 977 double rmsd = MultipleAlignmentScorer.getRMSD(symm); 978 979 symm.putScore(MultipleAlignmentScorer.AVGTM_SCORE, tmScore); 980 symm.putScore(MultipleAlignmentScorer.RMSD, rmsd); 981 } 982 983 /** 984 * Returns the representative Atom Array of the first model, if the 985 * structure is NMR, or the Array for each model, if it is a biological 986 * assembly with multiple models. 987 * 988 * @param structure 989 * @return representative Atom[] 990 */ 991 public static Atom[] getRepresentativeAtoms(Structure structure) { 992 993 if (structure.isNmr()) 994 return StructureTools.getRepresentativeAtomArray(structure); 995 996 else { 997 998 // Get Atoms of all models 999 List<Atom> atomList = new ArrayList<Atom>(); 1000 for (int m = 0; m < structure.nrModels(); m++) { 1001 for (Chain c : structure.getModel(m)) 1002 atomList.addAll(Arrays.asList(StructureTools 1003 .getRepresentativeAtomArray(c))); 1004 } 1005 return atomList.toArray(new Atom[0]); 1006 } 1007 1008 } 1009 1010}