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.MultipleAlignmentTools; 057import org.biojava.nbio.structure.align.multiple.util.MultipleSuperimposer; 058import org.biojava.nbio.structure.cluster.Subunit; 059import org.biojava.nbio.structure.cluster.SubunitCluster; 060import org.biojava.nbio.structure.cluster.SubunitClustererMethod; 061import org.biojava.nbio.structure.cluster.SubunitClustererParameters; 062import org.biojava.nbio.structure.geometry.SuperPositions; 063import org.biojava.nbio.structure.jama.Matrix; 064import org.biojava.nbio.structure.symmetry.core.QuatSymmetryDetector; 065import org.biojava.nbio.structure.symmetry.core.QuatSymmetryParameters; 066import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults; 067import org.biojava.nbio.structure.symmetry.core.Stoichiometry; 068import org.biojava.nbio.structure.symmetry.internal.CeSymmResult; 069import org.biojava.nbio.structure.symmetry.internal.SymmetryAxes; 070import org.jgrapht.Graph; 071import org.jgrapht.graph.DefaultEdge; 072import org.jgrapht.graph.SimpleGraph; 073import org.slf4j.Logger; 074import org.slf4j.LoggerFactory; 075 076/** 077 * Utility methods for symmetry (quaternary and internal) detection and result 078 * manipulation. 079 * 080 * @author Spencer Bliven 081 * @author Aleix Lafita 082 * @author Peter Rose 083 * 084 */ 085public class SymmetryTools { 086 087 private static final Logger logger = LoggerFactory 088 .getLogger(SymmetryTools.class); 089 090 /** Prevent instantiation. */ 091 private SymmetryTools() { 092 } 093 094 /** 095 * Returns the "reset value" for graying out the main diagonal. If we're 096 * blanking out the main diagonal, this value is always Integer.MIN_VALUE. 097 * <p> 098 * This is possible if {@code gradientPolyCoeff = Integer.MIN_VALUE} and 099 * {@code gradientExpCoeff = 0}. 100 * 101 * @param unpenalizedScore 102 * @param nResFromMainDiag 103 * @param gradientPolyCoeff 104 * @param gradientExpCoeff 105 * @return 106 */ 107 private static double getResetVal(double unpenalizedScore, 108 double nResFromMainDiag, double[] gradientPolyCoeff, 109 double gradientExpCoeff) { 110 111 if (Double.isNaN(unpenalizedScore)) 112 return 0; // what else? 113 114 // We can actually return a positive value if this is high enough 115 double updateVal = unpenalizedScore; 116 updateVal -= gradientExpCoeff * Math.pow(Math.E, -nResFromMainDiag); 117 for (int p = 0; p < gradientPolyCoeff.length; p++) { 118 updateVal -= gradientPolyCoeff[gradientPolyCoeff.length - 1 - p] 119 * Math.pow(nResFromMainDiag, -p); 120 } 121 return updateVal; 122 } 123 124 /** 125 * Grays out the main diagonal of a duplicated distance matrix. 126 * 127 * @param ca2 128 * @param rows 129 * Number of rows 130 * @param cols 131 * Number of original columns 132 * @param calculator 133 * Used to get the matrix if origM is null 134 * @param origM 135 * starting matrix. If null, uses 136 * {@link CECalculator#getMatMatrix()} 137 * @param blankWindowSize 138 * Width of section to gray out 139 * @param gradientPolyCoeff 140 * @param gradientExpCoeff 141 * @return 142 */ 143 public static Matrix grayOutCEOrig(Atom[] ca2, int rows, int cols, 144 CECalculator calculator, Matrix origM, int blankWindowSize, 145 double[] gradientPolyCoeff, double gradientExpCoeff) { 146 147 if (origM == null) { 148 origM = new Matrix(calculator.getMatMatrix()); 149 } 150 151 // symmetry hack, disable main diagonal 152 153 for (int i = 0; i < rows; i++) { 154 for (int j = 0; j < cols; j++) { 155 int diff = Math.abs(i - j); 156 157 double resetVal = getResetVal(origM.get(i, j), diff, 158 gradientPolyCoeff, gradientExpCoeff); 159 160 if (diff < blankWindowSize) { 161 origM.set(i, j, origM.get(i, j) + resetVal); 162 163 } 164 int diff2 = Math.abs(i - (j - ca2.length / 2)); // other side 165 166 double resetVal2 = getResetVal(origM.get(i, j), diff2, 167 gradientPolyCoeff, gradientExpCoeff); 168 169 if (diff2 < blankWindowSize) { 170 origM.set(i, j, origM.get(i, j) + resetVal2); 171 172 } 173 } 174 } 175 return origM; 176 } 177 178 public static Matrix grayOutPreviousAlignment(AFPChain afpChain, 179 Atom[] ca2, int rows, int cols, CECalculator calculator, 180 Matrix max, int blankWindowSize, double[] gradientPolyCoeff, 181 double gradientExpCoeff) { 182 183 max = grayOutCEOrig(ca2, rows, cols, calculator, max, blankWindowSize, 184 gradientPolyCoeff, gradientExpCoeff); 185 186 double[][] dist1 = calculator.getDist1(); 187 double[][] dist2 = calculator.getDist2(); 188 189 int[][][] optAln = afpChain.getOptAln(); 190 int blockNum = afpChain.getBlockNum(); 191 192 int[] optLen = afpChain.getOptLen(); 193 194 // ca2 is circularly permutated 195 int breakPoint = ca2.length / 2; 196 for (int bk = 0; bk < blockNum; bk++) { 197 198 for (int i = 0; i < optLen[bk]; i++) { 199 int pos1 = optAln[bk][0][i]; 200 int pos2 = optAln[bk][1][i]; 201 202 int dist = blankWindowSize / 2; 203 int start1 = Math.max(pos1 - dist, 0); 204 int start2 = Math.max(pos2 - dist, 0); 205 int end1 = Math.min(pos1 + dist, rows - 1); 206 int end2 = Math.min(pos2 + dist, cols - 1); 207 208 for (int i1 = start1; i1 < end1; i1++) { 209 210 // blank diagonal of dist1 211 for (int k = 0; k < blankWindowSize / 2; k++) { 212 if (i1 - k >= 0) { 213 double resetVal = getResetVal( 214 max.get(i1 - k, i1 - k), 0, 215 gradientPolyCoeff, gradientExpCoeff); 216 dist1[i1 - k][i1 - k] = resetVal; 217 } else if (i1 + k < rows) { 218 double resetVal = getResetVal( 219 max.get(i1 + k, i1 + k), 0, 220 gradientPolyCoeff, gradientExpCoeff); 221 dist1[i1 + k][i1 + k] = resetVal; 222 } 223 224 } 225 226 for (int j2 = start2; j2 < end2; j2++) { 227 double resetVal = getResetVal(max.get(i1, j2), 228 Math.abs(i1 - j2), gradientPolyCoeff, 229 gradientExpCoeff); 230 max.set(i1, j2, resetVal); 231 if (j2 < breakPoint) { 232 double resetVal2 = getResetVal( 233 max.get(i1, j2 + breakPoint), 234 Math.abs(i1 - (j2 + breakPoint)), 235 gradientPolyCoeff, gradientExpCoeff); 236 max.set(i1, j2 + breakPoint, resetVal2); 237 } else { 238 double resetVal2 = getResetVal( 239 max.get(i1, j2 - breakPoint), 240 Math.abs(i1 - (j2 - breakPoint)), 241 gradientPolyCoeff, gradientExpCoeff); 242 max.set(i1, j2 - breakPoint, resetVal2); 243 } 244 for (int k = 0; k < blankWindowSize / 2; k++) { 245 if (j2 - k >= 0) { 246 if (j2 - k < breakPoint) { 247 double resetVal2 = getResetVal( 248 max.get(j2 - k, j2 - k), 0, 249 gradientPolyCoeff, gradientExpCoeff); 250 dist2[j2 - k][j2 - k] = resetVal2; 251 } else { 252 double resetVal2 = getResetVal(max.get(j2 253 - k - breakPoint, j2 - k), 0, 254 gradientPolyCoeff, gradientExpCoeff); 255 dist2[j2 - k - breakPoint][j2 - k 256 - breakPoint] = resetVal2; 257 } 258 } else if (j2 + k < cols) { 259 if (j2 + k < breakPoint) { 260 double resetVal2 = getResetVal( 261 max.get(j2 + k, j2 + k), 0, 262 gradientPolyCoeff, gradientExpCoeff); 263 dist2[j2 + k][j2 + k] = resetVal2; 264 } else { 265 double resetVal2 = getResetVal(max.get(j2 266 + k - breakPoint, j2 + k), 0, 267 gradientPolyCoeff, gradientExpCoeff); 268 dist2[j2 + k - breakPoint][j2 + k 269 - breakPoint] = resetVal2; 270 } 271 } 272 } 273 } 274 } 275 276 } 277 } 278 calculator.setDist1(dist1); 279 calculator.setDist2(dist2); 280 return max; 281 282 } 283 284 public Matrix getDkMatrix(Atom[] ca1, Atom[] ca2, int fragmentLength, 285 double[] dist1, double[] dist2, int rows, int cols) { 286 287 Matrix diffDistMax = Matrix.identity(ca1.length, ca2.length); 288 289 for (int i = 0; i < rows; i++) { 290 double score1 = 0; 291 for (int x = 0; x < fragmentLength; x++) { 292 score1 += dist1[i + x]; 293 } 294 for (int j = 0; j < cols; j++) { 295 double score2 = 0; 296 for (int y = 0; y < fragmentLength; y++) { 297 score2 += dist2[j + y]; 298 } 299 300 // if the intramolecular distances are very similar 301 // the two scores should be similar, 302 // i.e. the difference is close to 0 303 diffDistMax.set(i, j, Math.abs(score1 - score2)); 304 } 305 } 306 307 // symmetry hack, disable main diagonal 308 309 for (int i = 0; i < rows; i++) { 310 for (int j = 0; j < cols; j++) { 311 int diff = Math.abs(i - j); 312 313 if (diff < 15) { 314 diffDistMax.set(i, j, 99); 315 } 316 int diff2 = Math.abs(i - (j - ca2.length / 2)); 317 if (diff2 < 15) { 318 diffDistMax.set(i, j, 99); 319 } 320 } 321 } 322 return diffDistMax; 323 324 } 325 326 public static Matrix blankOutPreviousAlignment(AFPChain afpChain, 327 Atom[] ca2, int rows, int cols, CECalculator calculator, 328 Matrix max, int blankWindowSize) { 329 return grayOutPreviousAlignment(afpChain, ca2, rows, cols, calculator, 330 max, blankWindowSize, new double[] { Integer.MIN_VALUE }, 0.0); 331 } 332 333 public static Matrix blankOutCEOrig(Atom[] ca2, int rows, int cols, 334 CECalculator calculator, Matrix origM, int blankWindowSize) { 335 return grayOutCEOrig(ca2, rows, cols, calculator, origM, 336 blankWindowSize, new double[] { Integer.MIN_VALUE }, 0.0); 337 } 338 339 public static Matrix getDkMatrix(Atom[] ca1, Atom[] ca2, int k, 340 int fragmentLength) { 341 342 double[] dist1 = AlignUtils.getDiagonalAtK(ca1, k); 343 double[] dist2 = AlignUtils.getDiagonalAtK(ca2, k); 344 345 int rows = ca1.length - fragmentLength - k + 1; 346 int cols = ca2.length - fragmentLength - k + 1; 347 348 // Matrix that tracks similarity of a fragment of length fragmentLength 349 // starting a position i,j. 350 351 Matrix m2 = new Matrix(rows, cols); 352 353 for (int i = 0; i < rows; i++) { 354 double score1 = 0; 355 for (int x = 0; x < fragmentLength; x++) { 356 score1 += dist1[i + x]; 357 } 358 for (int j = 0; j < cols; j++) { 359 double score2 = 0; 360 for (int y = 0; y < fragmentLength; y++) { 361 score2 += dist2[j + y]; 362 } 363 364 // if the intramolecular distances are very similar 365 // the two scores should be similar, 366 // i.e. the difference is close to 0 367 m2.set(i, j, Math.abs(score1 - score2)); 368 } 369 } 370 return m2; 371 } 372 373 public static boolean[][] blankOutBreakFlag(AFPChain afpChain, Atom[] ca2, 374 int rows, int cols, CECalculator calculator, boolean[][] breakFlag, 375 int blankWindowSize) { 376 377 int[][][] optAln = afpChain.getOptAln(); 378 int blockNum = afpChain.getBlockNum(); 379 380 int[] optLen = afpChain.getOptLen(); 381 382 // ca2 is circularly permutated at this point. 383 int breakPoint = ca2.length / 2; 384 385 for (int bk = 0; bk < blockNum; bk++) { 386 387 // Matrix m= afpChain.getBlockRotationMatrix()[bk]; 388 // Atom shift = afpChain.getBlockShiftVector()[bk]; 389 for (int i = 0; i < optLen[bk]; i++) { 390 int pos1 = optAln[bk][0][i]; 391 int pos2 = optAln[bk][1][i]; 392 // blank out area around these positions... 393 394 int dist = blankWindowSize; 395 int start1 = Math.max(pos1 - dist, 0); 396 int start2 = Math.max(pos2 - dist, 0); 397 int end1 = Math.min(pos1 + dist, rows - 1); 398 int end2 = Math.min(pos2 + dist, cols - 1); 399 400 for (int i1 = start1; i1 < end1; i1++) { 401 402 for (int j2 = start2; j2 < end2; j2++) { 403 // System.out.println(i1 + " " + j2 + " (***)"); 404 breakFlag[i1][j2] = true; 405 if (j2 < breakPoint) { 406 breakFlag[i1][j2 + breakPoint] = true; 407 } 408 } 409 } 410 411 } 412 } 413 414 return breakFlag; 415 } 416 417 /** 418 * Returns the <em>magnitude</em> of the angle between the first and second 419 * blocks of {@code afpChain}, measured in degrees. This is always a 420 * positive value (unsigned). 421 * 422 * @param afpChain 423 * @param ca1 424 * @param ca2 425 * @return 426 */ 427 public static double getAngle(AFPChain afpChain, Atom[] ca1, Atom[] ca2) { 428 Matrix rotation = afpChain.getBlockRotationMatrix()[0]; 429 return Math.acos(rotation.trace() - 1) * 180 / Math.PI; 430 } 431 432 /** 433 * Converts a set of AFP alignments into a Graph of aligned residues, where 434 * each vertex is a residue and each edge means the connection between the 435 * two residues in one of the alignments. 436 * 437 * @param afps 438 * List of AFPChains 439 * @param atoms 440 * Atom array of the symmetric structure 441 * @param undirected 442 * if true, the graph is undirected 443 * 444 * @return adjacency List of aligned residues 445 */ 446 public static List<List<Integer>> buildSymmetryGraph(List<AFPChain> afps, 447 Atom[] atoms, boolean undirected) { 448 449 List<List<Integer>> graph = new ArrayList<>(); 450 451 for (int n = 0; n < atoms.length; n++) { 452 graph.add(new ArrayList<Integer>()); 453 } 454 455 for (int k = 0; k < afps.size(); k++) { 456 for (int i = 0; i < afps.get(k).getOptAln().length; i++) { 457 for (int j = 0; j < afps.get(k).getOptAln()[i][0].length; j++) { 458 Integer res1 = afps.get(k).getOptAln()[i][0][j]; 459 Integer res2 = afps.get(k).getOptAln()[i][1][j]; 460 graph.get(res1).add(res2); 461 if (undirected) 462 graph.get(res2).add(res1); 463 } 464 } 465 } 466 return graph; 467 } 468 469 /** 470 * Converts a self alignment into a directed jGraphT of aligned residues, 471 * where each vertex is a residue and each edge means the equivalence 472 * between the two residues in the self-alignment. 473 * 474 * @param selfAlignment 475 * AFPChain 476 * 477 * @return alignment Graph 478 */ 479 public static Graph<Integer, DefaultEdge> buildSymmetryGraph( 480 AFPChain selfAlignment) { 481 482 Graph<Integer, DefaultEdge> graph = new SimpleGraph<Integer, DefaultEdge>( 483 DefaultEdge.class); 484 485 for (int i = 0; i < selfAlignment.getOptAln().length; i++) { 486 for (int j = 0; j < selfAlignment.getOptAln()[i][0].length; j++) { 487 Integer res1 = selfAlignment.getOptAln()[i][0][j]; 488 Integer res2 = selfAlignment.getOptAln()[i][1][j]; 489 graph.addVertex(res1); 490 graph.addVertex(res2); 491 graph.addEdge(res1, res2); 492 } 493 } 494 return graph; 495 } 496 497 /** 498 * Method that converts the symmetric units of a structure into different 499 * structures, so that they can be individually visualized. 500 * 501 * @param symmetry 502 * CeSymmResult 503 * @throws StructureException 504 * @return List of structures, by repeat index sequentially 505 * 506 */ 507 public static List<Structure> divideStructure(CeSymmResult symmetry) 508 throws StructureException { 509 510 if (!symmetry.isRefined()) 511 throw new IllegalArgumentException("The symmetry result " 512 + "is not refined, repeats cannot be defined"); 513 514 int order = symmetry.getMultipleAlignment().size(); 515 Atom[] atoms = symmetry.getAtoms(); 516 Set<Group> allGroups = StructureTools.getAllGroupsFromSubset(atoms, GroupType.HETATM); 517 List<StructureIdentifier> repeatsId = symmetry.getRepeatsID(); 518 List<Structure> repeats = new ArrayList<>(order); 519 520 // Create new structure containing the repeat atoms 521 for (int i = 0; i < order; i++) { 522 523 Structure s = new StructureImpl(); 524 s.addModel(new ArrayList<Chain>(1)); 525 s.setStructureIdentifier(repeatsId.get(i)); 526 527 Block align = symmetry.getMultipleAlignment().getBlock(0); 528 529 // Get the start and end of the repeat 530 // Repeats are always sequential blocks 531 int res1 = align.getStartResidue(i); 532 int res2 = align.getFinalResidue(i); 533 534 // All atoms from the repeat, used for ligand search 535 // AA have an average of 8.45 atoms, so guess capacity with that 536 List<Atom> repeat = new ArrayList<>(Math.max(9*(res2-res1+1),9)); 537 // speedy chain lookup 538 Chain prevChain = null; 539 for(int k=res1;k<=res2; k++) { 540 Group g = atoms[k].getGroup(); 541 prevChain = StructureTools.addGroupToStructure(s, g, 0, prevChain,true); 542 repeat.addAll(g.getAtoms()); 543 } 544 545 546 List<Group> ligands = StructureTools.getLigandsByProximity( 547 allGroups, 548 repeat.toArray(new Atom[repeat.size()]), 549 StructureTools.DEFAULT_LIGAND_PROXIMITY_CUTOFF); 550 551 logger.warn("Adding {} ligands to {}",ligands.size(), symmetry.getMultipleAlignment().getStructureIdentifier(i)); 552 for( Group ligand : ligands) { 553 prevChain = StructureTools.addGroupToStructure(s, ligand, 0, prevChain,true); 554 } 555 556 repeats.add(s); 557 } 558 return repeats; 559 } 560 561 /** 562 * Method that converts a repeats symmetric alignment into an alignment of 563 * whole structures. 564 * <p> 565 * Example: if the structure has repeats A,B and C, the original alignment 566 * is A-B-C, and the returned alignment is ABC-BCA-CAB. 567 * 568 * @param symm 569 * CeSymmResult 570 * @return MultipleAlignment of the full structure superpositions 571 */ 572 public static MultipleAlignment toFullAlignment(CeSymmResult symm) { 573 574 if (!symm.isRefined()) 575 throw new IllegalArgumentException("The symmetry result " 576 + "is not refined, repeats cannot be defined"); 577 578 MultipleAlignment full = symm.getMultipleAlignment().clone(); 579 580 for (int str = 1; str < full.size(); str++) { 581 // Create a new Block with swapped AlignRes (move first to last) 582 Block b = full.getBlock(full.getBlocks().size() - 1).clone(); 583 b.getAlignRes().add(b.getAlignRes().get(0)); 584 b.getAlignRes().remove(0); 585 full.getBlockSet(0).getBlocks().add(b); 586 } 587 return full; 588 } 589 590 /** 591 * Method that converts a symmetry alignment into an alignment of the 592 * repeats only, as new independent structures. 593 * <p> 594 * This method changes the structure identifiers, the Atom arrays and 595 * re-scles the aligned residues in the Blocks corresponding to those 596 * changes. 597 * <p> 598 * Application: display superimposed repeats in Jmol. 599 * 600 * @param result 601 * CeSymmResult of symmetry 602 * @return MultipleAlignment of the repeats 603 * @throws StructureException 604 */ 605 public static MultipleAlignment toRepeatsAlignment(CeSymmResult result) 606 throws StructureException { 607 608 if (!result.isRefined()) 609 throw new IllegalArgumentException("The symmetry result " 610 + "is not refined, repeats cannot be defined"); 611 612 MultipleAlignment msa = result.getMultipleAlignment(); 613 MultipleAlignmentEnsemble newEnsemble = msa.getEnsemble().clone(); 614 615 List<Structure> repSt = SymmetryTools.divideStructure(result); 616 617 MultipleAlignment repeats = newEnsemble.getMultipleAlignment(0); 618 Block block = repeats.getBlock(0); 619 List<Atom[]> atomArrays = new ArrayList<>(); 620 621 for (Structure s : repSt) 622 atomArrays.add(StructureTools.getRepresentativeAtomArray(s)); 623 624 newEnsemble.setAtomArrays(atomArrays); 625 626 for (int su = 0; su < block.size(); su++) { 627 Integer start = block.getStartResidue(su); 628 629 // Normalize aligned residues 630 for (int res = 0; res < block.length(); res++) { 631 Integer residue = block.getAlignRes().get(su).get(res); 632 if (residue != null) 633 residue -= start; 634 block.getAlignRes().get(su).set(res, residue); 635 } 636 } 637 638 return repeats; 639 } 640 641 /** 642 * Converts a refined symmetry AFPChain alignment into the standard 643 * representation of symmetry in a MultipleAlignment, that contains the 644 * entire Atom array of the strcuture and the symmetric repeats are orgaized 645 * in different rows in a single Block. 646 * 647 * @param symm 648 * AFPChain created with a symmetry algorithm and refined 649 * @param atoms 650 * Atom array of the entire structure 651 * @return MultipleAlignment format of the symmetry 652 * @throws StructureException 653 */ 654 public static MultipleAlignment fromAFP(AFPChain symm, Atom[] atoms) 655 throws StructureException { 656 657 if (!symm.getAlgorithmName().contains("symm")) { 658 throw new IllegalArgumentException( 659 "The input alignment is not a symmetry alignment."); 660 } 661 662 MultipleAlignmentEnsemble e = new MultipleAlignmentEnsembleImpl(symm, 663 atoms, atoms, false); 664 e.setAtomArrays(new ArrayList<Atom[]>()); 665 StructureIdentifier name = null; 666 if (e.getStructureIdentifiers() != null) { 667 if (!e.getStructureIdentifiers().isEmpty()) 668 name = e.getStructureIdentifiers().get(0); 669 } else 670 name = atoms[0].getGroup().getChain().getStructure() 671 .getStructureIdentifier(); 672 673 e.setStructureIdentifiers(new ArrayList<StructureIdentifier>()); 674 675 MultipleAlignment result = new MultipleAlignmentImpl(); 676 BlockSet bs = new BlockSetImpl(result); 677 Block b = new BlockImpl(bs); 678 b.setAlignRes(new ArrayList<List<Integer>>()); 679 680 int order = symm.getBlockNum(); 681 for (int su = 0; su < order; su++) { 682 List<Integer> residues = e.getMultipleAlignment(0).getBlock(su) 683 .getAlignRes().get(0); 684 b.getAlignRes().add(residues); 685 e.getStructureIdentifiers().add(name); 686 e.getAtomArrays().add(atoms); 687 } 688 e.getMultipleAlignments().set(0, result); 689 result.setEnsemble(e); 690 691 CoreSuperimposer imposer = new CoreSuperimposer(); 692 imposer.superimpose(result); 693 updateSymmetryScores(result); 694 695 return result; 696 } 697 698 /** 699 * Given a symmetry result, it calculates the overall global symmetry, 700 * factoring out the alignment and detection steps of 701 * {@link QuatSymmetryDetector} algorithm. 702 * 703 * @param result 704 * symmetry result 705 * @return global symmetry results 706 * @throws StructureException 707 */ 708 public static QuatSymmetryResults getQuaternarySymmetry(CeSymmResult result) 709 throws StructureException { 710 711 // Obtain the subunits of the repeats 712 MultipleAlignment msa = toRepeatsAlignment(result); 713 List<Atom[]> atoms = msa.getAtomArrays(); 714 List<Subunit> subunits = atoms.stream() 715 .map(a -> new Subunit(a, null, null, null)) 716 .collect(Collectors.toList()); 717 List<List<Integer>> eqr = MultipleAlignmentTools.getEquivalentResidues(msa, true); 718 SubunitCluster cluster = new SubunitCluster(subunits, eqr); 719 Stoichiometry composition = new Stoichiometry(Arrays.asList(cluster)); 720 721 QuatSymmetryParameters sp = new QuatSymmetryParameters(); 722 QuatSymmetryResults gSymmetry = QuatSymmetryDetector 723 .calcGlobalSymmetry(composition, sp); 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<>(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<>(); 778 List<Atom> list2 = new ArrayList<>(); 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<>(); 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<>(); 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<>(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}