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.addModel(new ArrayList<Chain>(1)); 522 s.setStructureIdentifier(repeatsId.get(i)); 523 524 Block align = symmetry.getMultipleAlignment().getBlock(0); 525 526 // Get the start and end of the repeat 527 // Repeats are always sequential blocks 528 int res1 = align.getStartResidue(i); 529 int res2 = align.getFinalResidue(i); 530 531 // All atoms from the repeat, used for ligand search 532 // AA have an average of 8.45 atoms, so guess capacity with that 533 List<Atom> repeat = new ArrayList<>(Math.max(9*(res2-res1+1),9)); 534 // speedy chain lookup 535 Chain prevChain = null; 536 for(int k=res1;k<=res2; k++) { 537 Group g = atoms[k].getGroup(); 538 prevChain = StructureTools.addGroupToStructure(s, g, 0, prevChain,true); 539 repeat.addAll(g.getAtoms()); 540 } 541 542 543 List<Group> ligands = StructureTools.getLigandsByProximity( 544 allGroups, 545 repeat.toArray(new Atom[repeat.size()]), 546 StructureTools.DEFAULT_LIGAND_PROXIMITY_CUTOFF); 547 548 logger.warn("Adding {} ligands to {}",ligands.size(), symmetry.getMultipleAlignment().getStructureIdentifier(i)); 549 for( Group ligand : ligands) { 550 prevChain = StructureTools.addGroupToStructure(s, ligand, 0, prevChain,true); 551 } 552 553 repeats.add(s); 554 } 555 return repeats; 556 } 557 558 /** 559 * Method that converts a repeats symmetric alignment into an alignment of 560 * whole structures. 561 * <p> 562 * Example: if the structure has repeats A,B and C, the original alignment 563 * is A-B-C, and the returned alignment is ABC-BCA-CAB. 564 * 565 * @param symm 566 * CeSymmResult 567 * @return MultipleAlignment of the full structure superpositions 568 */ 569 public static MultipleAlignment toFullAlignment(CeSymmResult symm) { 570 571 if (!symm.isRefined()) 572 throw new IllegalArgumentException("The symmetry result " 573 + "is not refined, repeats cannot be defined"); 574 575 MultipleAlignment full = symm.getMultipleAlignment().clone(); 576 577 for (int str = 1; str < full.size(); str++) { 578 // Create a new Block with swapped AlignRes (move first to last) 579 Block b = full.getBlock(full.getBlocks().size() - 1).clone(); 580 b.getAlignRes().add(b.getAlignRes().get(0)); 581 b.getAlignRes().remove(0); 582 full.getBlockSet(0).getBlocks().add(b); 583 } 584 return full; 585 } 586 587 /** 588 * Method that converts a symmetry alignment into an alignment of the 589 * repeats only, as new independent structures. 590 * <p> 591 * This method changes the structure identifiers, the Atom arrays and 592 * re-scles the aligned residues in the Blocks corresponding to those 593 * changes. 594 * <p> 595 * Application: display superimposed repeats in Jmol. 596 * 597 * @param result 598 * CeSymmResult of symmetry 599 * @return MultipleAlignment of the repeats 600 * @throws StructureException 601 */ 602 public static MultipleAlignment toRepeatsAlignment(CeSymmResult result) 603 throws StructureException { 604 605 if (!result.isRefined()) 606 throw new IllegalArgumentException("The symmetry result " 607 + "is not refined, repeats cannot be defined"); 608 609 MultipleAlignment msa = result.getMultipleAlignment(); 610 MultipleAlignmentEnsemble newEnsemble = msa.getEnsemble().clone(); 611 612 List<Structure> repSt = SymmetryTools.divideStructure(result); 613 614 MultipleAlignment repeats = newEnsemble.getMultipleAlignment(0); 615 Block block = repeats.getBlock(0); 616 List<Atom[]> atomArrays = new ArrayList<Atom[]>(); 617 618 for (Structure s : repSt) 619 atomArrays.add(StructureTools.getRepresentativeAtomArray(s)); 620 621 newEnsemble.setAtomArrays(atomArrays); 622 623 for (int su = 0; su < block.size(); su++) { 624 Integer start = block.getStartResidue(su); 625 626 // Normalize aligned residues 627 for (int res = 0; res < block.length(); res++) { 628 Integer residue = block.getAlignRes().get(su).get(res); 629 if (residue != null) 630 residue -= start; 631 block.getAlignRes().get(su).set(res, residue); 632 } 633 } 634 635 return repeats; 636 } 637 638 /** 639 * Converts a refined symmetry AFPChain alignment into the standard 640 * representation of symmetry in a MultipleAlignment, that contains the 641 * entire Atom array of the strcuture and the symmetric repeats are orgaized 642 * in different rows in a single Block. 643 * 644 * @param symm 645 * AFPChain created with a symmetry algorithm and refined 646 * @param atoms 647 * Atom array of the entire structure 648 * @return MultipleAlignment format of the symmetry 649 * @throws StructureException 650 */ 651 public static MultipleAlignment fromAFP(AFPChain symm, Atom[] atoms) 652 throws StructureException { 653 654 if (!symm.getAlgorithmName().contains("symm")) { 655 throw new IllegalArgumentException( 656 "The input alignment is not a symmetry alignment."); 657 } 658 659 MultipleAlignmentEnsemble e = new MultipleAlignmentEnsembleImpl(symm, 660 atoms, atoms, false); 661 e.setAtomArrays(new ArrayList<Atom[]>()); 662 StructureIdentifier name = null; 663 if (e.getStructureIdentifiers() != null) { 664 if (!e.getStructureIdentifiers().isEmpty()) 665 name = e.getStructureIdentifiers().get(0); 666 } else 667 name = atoms[0].getGroup().getChain().getStructure() 668 .getStructureIdentifier(); 669 670 e.setStructureIdentifiers(new ArrayList<StructureIdentifier>()); 671 672 MultipleAlignment result = new MultipleAlignmentImpl(); 673 BlockSet bs = new BlockSetImpl(result); 674 Block b = new BlockImpl(bs); 675 b.setAlignRes(new ArrayList<List<Integer>>()); 676 677 int order = symm.getBlockNum(); 678 for (int su = 0; su < order; su++) { 679 List<Integer> residues = e.getMultipleAlignment(0).getBlock(su) 680 .getAlignRes().get(0); 681 b.getAlignRes().add(residues); 682 e.getStructureIdentifiers().add(name); 683 e.getAtomArrays().add(atoms); 684 } 685 e.getMultipleAlignments().set(0, result); 686 result.setEnsemble(e); 687 688 CoreSuperimposer imposer = new CoreSuperimposer(); 689 imposer.superimpose(result); 690 updateSymmetryScores(result); 691 692 return result; 693 } 694 695 /** 696 * Given a symmetry result, it calculates the overall global symmetry, 697 * factoring out the alignment and detection steps of 698 * {@link QuatSymmetryDetector} algorithm. 699 * 700 * @param result 701 * symmetry result 702 * @return global symmetry results 703 * @throws StructureException 704 */ 705 public static QuatSymmetryResults getQuaternarySymmetry(CeSymmResult result) 706 throws StructureException { 707 708 // Obtain the subunits of the repeats 709 List<Atom[]> atoms = toRepeatsAlignment(result).getAtomArrays(); 710 List<Subunit> subunits = atoms.stream() 711 .map(a -> new Subunit(a, null, null, null)) 712 .collect(Collectors.toList()); 713 714 // The clustering thresholds are set to 0 so that all always merged 715 SubunitClustererParameters cp = new SubunitClustererParameters(); 716 cp.setClustererMethod(SubunitClustererMethod.STRUCTURE); 717 cp.setRMSDThreshold(10.0); 718 cp.setStructureCoverageThreshold(0.0); 719 720 QuatSymmetryParameters sp = new QuatSymmetryParameters(); 721 722 QuatSymmetryResults gSymmetry = QuatSymmetryDetector 723 .calcGlobalSymmetry(subunits, sp, cp); 724 725 return gSymmetry; 726 } 727 728 /** 729 * Returns the List of Groups of the corresponding representative Atom 730 * array. The representative Atom array needs to fulfill: no two Atoms are 731 * from the same Group and Groups are sequential (connected in the original 732 * Structure), except if they are from different Chains. 733 * 734 * @param rAtoms 735 * array of representative Atoms (CA, P, etc). 736 * @return List of Groups 737 */ 738 public static List<Group> getGroups(Atom[] rAtoms) { 739 740 List<Group> groups = new ArrayList<Group>(rAtoms.length); 741 742 for (Atom a : rAtoms) { 743 Group g = a.getGroup(); 744 if (g != null) 745 groups.add(g); 746 else 747 logger.info("Group not found for representative Atom {}", a); 748 } 749 return groups; 750 } 751 752 /** 753 * Calculates the set of symmetry operation Matrices (transformations) of 754 * the new alignment, based on the symmetry relations in the SymmetryAxes 755 * object. It sets the transformations to the input MultipleAlignment and 756 * SymmetryAxes objects. If the SymmetryAxes object is null, the 757 * superposition of the repeats is done without symmetry constraints. 758 * <p> 759 * This method also sets the scores (RMSD and TM-score) after the new 760 * superposition has been updated. 761 * 762 * @param axes 763 * SymmetryAxes object. It will be modified. 764 * @param msa 765 * MultipleAlignment. It will be modified. 766 */ 767 public static void updateSymmetryTransformation(SymmetryAxes axes, 768 MultipleAlignment msa) throws StructureException { 769 770 List<List<Integer>> block = msa.getBlocks().get(0).getAlignRes(); 771 int length = block.get(0).size(); 772 773 if (axes != null) { 774 for (int level = 0; level < axes.getNumLevels(); level++) { 775 776 // Calculate the aligned atom arrays to superimpose 777 List<Atom> list1 = new ArrayList<Atom>(); 778 List<Atom> list2 = new ArrayList<Atom>(); 779 780 for (int firstRepeat : axes.getFirstRepeats(level)) { 781 782 Matrix4d transform = axes.getRepeatTransform(firstRepeat); 783 784 List<List<Integer>> relation = axes.getRepeatRelation( 785 level, firstRepeat); 786 787 for (int index = 0; index < relation.get(0).size(); index++) { 788 int p1 = relation.get(0).get(index); 789 int p2 = relation.get(1).get(index); 790 791 for (int k = 0; k < length; k++) { 792 Integer pos1 = block.get(p1).get(k); 793 Integer pos2 = block.get(p2).get(k); 794 if (pos1 != null && pos2 != null) { 795 Atom a = (Atom) msa.getAtomArrays().get(p1)[pos1] 796 .clone(); 797 Atom b = (Atom) msa.getAtomArrays().get(p2)[pos2] 798 .clone(); 799 Calc.transform(a, transform); 800 Calc.transform(b, transform); 801 list1.add(a); 802 list2.add(b); 803 } 804 } 805 } 806 } 807 808 Atom[] arr1 = list1.toArray(new Atom[list1.size()]); 809 Atom[] arr2 = list2.toArray(new Atom[list2.size()]); 810 811 // Calculate the new transformation information 812 if (arr1.length > 0 && arr2.length > 0) { 813 Matrix4d axis = SuperPositions.superpose( 814 Calc.atomsToPoints(arr1), 815 Calc.atomsToPoints(arr2)); 816 axes.updateAxis(level, axis); 817 } 818 819 // Get the transformations from the SymmetryAxes 820 List<Matrix4d> transformations = new ArrayList<Matrix4d>(); 821 for (int su = 0; su < msa.size(); su++) { 822 transformations.add(axes.getRepeatTransform(su)); 823 } 824 msa.getBlockSet(0).setTransformations(transformations); 825 } 826 } else { 827 MultipleSuperimposer imposer = new CoreSuperimposer(); 828 imposer.superimpose(msa); 829 } 830 updateSymmetryScores(msa); 831 } 832 833 /** 834 * Update the scores (TM-score and RMSD) of a symmetry multiple alignment. 835 * This method does not redo the superposition of the alignment. 836 * 837 * @param symm 838 * Symmetry Multiple Alignment of Repeats 839 * @throws StructureException 840 */ 841 public static void updateSymmetryScores(MultipleAlignment symm) 842 throws StructureException { 843 844 // Multiply by the order of symmetry to normalize score 845 double tmScore = MultipleAlignmentScorer.getAvgTMScore(symm) 846 * symm.size(); 847 double rmsd = MultipleAlignmentScorer.getRMSD(symm); 848 849 symm.putScore(MultipleAlignmentScorer.AVGTM_SCORE, tmScore); 850 symm.putScore(MultipleAlignmentScorer.RMSD, rmsd); 851 } 852 853 /** 854 * Returns the representative Atom Array of the first model, if the 855 * structure is NMR, or the Array for each model, if it is a biological 856 * assembly with multiple models. 857 * 858 * @param structure 859 * @return representative Atom[] 860 */ 861 public static Atom[] getRepresentativeAtoms(Structure structure) { 862 863 if (structure.isNmr()) 864 return StructureTools.getRepresentativeAtomArray(structure); 865 866 else { 867 868 // Get Atoms of all models 869 List<Atom> atomList = new ArrayList<Atom>(); 870 for (int m = 0; m < structure.nrModels(); m++) { 871 for (Chain c : structure.getModel(m)) 872 atomList.addAll(Arrays.asList(StructureTools 873 .getRepresentativeAtomArray(c))); 874 } 875 return atomList.toArray(new Atom[0]); 876 } 877 878 } 879 880 /** 881 * Find valid symmetry orders for a given stoichiometry. For instance, an 882 * A6B4 protein would give [1,2] because (A6B4)1 and (A3B2)2 are valid 883 * decompositions. 884 * 885 * @param stoichiometry 886 * List giving the number of copies in each Subunit cluster 887 * @return The common factors of the stoichiometry 888 */ 889 public static List<Integer> getValidFolds(List<Integer> stoichiometry) { 890 891 List<Integer> denominators = new ArrayList<Integer>(); 892 893 if (stoichiometry.isEmpty()) 894 return denominators; 895 896 int nChains = Collections.max(stoichiometry); 897 898 // Remove duplicate stoichiometries 899 Set<Integer> nominators = new TreeSet<Integer>(stoichiometry); 900 901 // find common denominators 902 for (int d = 1; d <= nChains; d++) { 903 boolean isDivisable = true; 904 for (Integer n : nominators) { 905 if (n % d != 0) { 906 isDivisable = false; 907 break; 908 } 909 } 910 if (isDivisable) { 911 denominators.add(d); 912 } 913 } 914 return denominators; 915 } 916 917}