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