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.internal; 022 023import java.util.ArrayList; 024import java.util.List; 025 026import javax.vecmath.Matrix4d; 027 028import org.biojava.nbio.structure.Atom; 029import org.biojava.nbio.structure.Structure; 030import org.biojava.nbio.structure.StructureException; 031import org.biojava.nbio.structure.StructureIdentifier; 032import org.biojava.nbio.structure.StructureTools; 033import org.biojava.nbio.structure.align.ce.CECalculator; 034import org.biojava.nbio.structure.align.ce.CeCPMain; 035import org.biojava.nbio.structure.align.ce.MatrixListener; 036import org.biojava.nbio.structure.align.model.AFPChain; 037import org.biojava.nbio.structure.align.multiple.MultipleAlignment; 038import org.biojava.nbio.structure.align.util.AFPChainScorer; 039import org.biojava.nbio.structure.jama.Matrix; 040import org.biojava.nbio.structure.secstruc.SecStrucCalc; 041import org.biojava.nbio.structure.secstruc.SecStrucTools; 042import org.biojava.nbio.structure.symmetry.internal.CESymmParameters.SymmetryType; 043import org.biojava.nbio.structure.symmetry.utils.SymmetryTools; 044import org.slf4j.Logger; 045import org.slf4j.LoggerFactory; 046 047/** 048 * Identify the symmetries in a structure by running an alignment of the 049 * structure against itself disabling the diagonal of the identity alignment. 050 * <p> 051 * Iterating over previous results and disabling the diagonal of the previous 052 * alignments can also be done with this implementation, which will generate a 053 * set of self-alignments (disabled, because none improvements were shown, but 054 * can be turn on manually). 055 * <p> 056 * Multiple levels of symmetry can be analyzed by finding symmetries in repeats 057 * of previous results. This feature allows to find multiple symmetry axes. 058 * <p> 059 * The alignment is then refined to obtain a consistent alignment among all 060 * residues of the structure and organized into different parts, called 061 * symmetric repeats. 062 * <p> 063 * After refinement of the initial alignment, an optimization step can be used 064 * to improve the overall score of the repeat multiple alignment. 065 * 066 * @author Andreas Prlic 067 * @author Spencer Bliven 068 * @author Aleix Lafita 069 * @since 4.1.1 070 * 071 */ 072public class CeSymm { 073 074 /** 075 * Version History: 076 * <p> 077 * <ul> 078 * <li>1.0 - initial implementation of CE-Symm. 079 * <li>1.1 - enable multiple CE-Symm runs to calculate all self-alignments. 080 * <li>2.0 - refine the alignment for consistency of repeat definition. 081 * <li>2.1 - optimize the alignment to improve the score. 082 * <li>2.2 - run multiple symmetry levels recursively to find PG and 083 * hierarchical symmetries. 084 * </ul> 085 * </li> 086 */ 087 public static final String version = "2.2"; 088 public static final String algorithmName = "jCE-symm"; 089 private static final Logger logger = LoggerFactory.getLogger(CeSymm.class); 090 private final static boolean multiPass = false; // multiple self-alignments 091 092 /** 093 * Prevent instantiation. Static class. 094 */ 095 private CeSymm() { 096 } 097 098 private static Matrix align(AFPChain afpChain, Atom[] ca1, Atom[] ca2, 099 CESymmParameters params, Matrix origM, CECalculator calculator, 100 int counter) throws StructureException { 101 102 int fragmentLength = params.getWinSize(); 103 Atom[] ca2clone = StructureTools.cloneAtomArray(ca2); 104 105 int rows = ca1.length; 106 int cols = ca2.length; 107 108 // Matrix that tracks similarity of a fragment of length fragmentLength 109 // starting a position i,j. 110 111 int blankWindowSize = fragmentLength; 112 if (origM == null) { 113 114 // Build alignment ca1 to ca2-ca2 115 afpChain = calculator.extractFragments(afpChain, ca1, ca2clone); 116 117 origM = SymmetryTools.blankOutPreviousAlignment(afpChain, ca2, 118 rows, cols, calculator, null, blankWindowSize); 119 120 } else { 121 // we are doing an iteration on a previous alignment 122 // mask the previous alignment 123 origM = SymmetryTools.blankOutPreviousAlignment(afpChain, ca2, 124 rows, cols, calculator, origM, blankWindowSize); 125 } 126 127 Matrix clone = (Matrix) origM.clone(); 128 129 // that's the matrix to run the alignment on.. 130 calculator.setMatMatrix(clone.getArray()); 131 132 calculator.traceFragmentMatrix(afpChain, ca1, ca2clone); 133 134 final Matrix origMfinal = (Matrix) origM.clone(); 135 // Add a matrix listener to keep the blacked zones in max. 136 calculator.addMatrixListener(new MatrixListener() { 137 138 @Override 139 public double[][] matrixInOptimizer(double[][] max) { 140 141 // Check every entry of origM for blacked out regions 142 for (int i = 0; i < max.length; i++) { 143 for (int j = 0; j < max[i].length; j++) { 144 if (origMfinal.getArray()[i][j] > 1e9) { 145 max[i][j] = -origMfinal.getArray()[i][j]; 146 } 147 } 148 } 149 return max; 150 } 151 152 @Override 153 public boolean[][] initializeBreakFlag(boolean[][] brkFlag) { 154 155 return brkFlag; 156 } 157 }); 158 159 calculator.nextStep(afpChain, ca1, ca2clone); 160 161 afpChain.setAlgorithmName(algorithmName); 162 afpChain.setVersion(version); 163 164 afpChain.setDistanceMatrix(origM); 165 166 return origMfinal; 167 168 } 169 170 @SuppressWarnings("unused") 171 protected static CeSymmResult align(Atom[] atoms, CESymmParameters params) 172 throws StructureException { 173 174 CeSymmResult result = new CeSymmResult(); 175 result.setParams(params); 176 result.setAtoms(atoms); 177 178 // STEP 1: prepare all the information for the symmetry alignment 179 Atom[] ca2 = StructureTools.duplicateCA2(atoms); 180 int rows = atoms.length; 181 int cols = ca2.length; 182 183 if (rows == 0 || cols == 0) { 184 throw new StructureException("Aligning empty structure"); 185 } 186 187 Matrix origM = null; 188 AFPChain myAFP = new AFPChain(algorithmName); 189 CECalculator calculator = new CECalculator(params); 190 Matrix lastMatrix = null; 191 192 List<AFPChain> selfAlignments = new ArrayList<AFPChain>(); 193 AFPChain optimalAFP = null; 194 195 // STEP 2: perform the self-alignments of the structure 196 int i = 0; 197 do { 198 if (origM != null) 199 myAFP.setDistanceMatrix((Matrix) origM.clone()); 200 201 origM = align(myAFP, atoms, ca2, params, origM, calculator, i); 202 203 double tmScore2 = AFPChainScorer.getTMScore(myAFP, atoms, ca2); 204 myAFP.setTMScore(tmScore2); 205 206 AFPChain newAFP = (AFPChain) myAFP.clone(); 207 newAFP = CeCPMain.postProcessAlignment(newAFP, atoms, ca2, 208 calculator); 209 210 // Calculate and set the TM score for the newAFP alignment 211 double tmScore3 = AFPChainScorer.getTMScore(newAFP, atoms, ca2); 212 newAFP.setTMScore(tmScore3); 213 214 // Determine if the alignment is significant, stop if false 215 if (tmScore3 < params.getUnrefinedScoreThreshold()) { 216 // If it is the first alignment save it anyway 217 if (i == 0) 218 selfAlignments.add(newAFP); 219 // store final matrix 220 lastMatrix = newAFP.getDistanceMatrix().copy(); 221 break; 222 } 223 224 // If it is a symmetric alignment add it to the allAlignments list 225 selfAlignments.add(newAFP); 226 227 i++; 228 229 } while (i < params.getMaxSymmOrder() && multiPass); 230 231 // We reached the maximum order, so blank out the final alignment 232 if (lastMatrix == null && selfAlignments.size() > 1 && multiPass) { 233 AFPChain last = selfAlignments.get(selfAlignments.size() - 1); 234 lastMatrix = SymmetryTools.blankOutPreviousAlignment(last, ca2, 235 last.getCa1Length(), last.getCa2Length(), calculator, 236 origM, params.getWinSize()); 237 lastMatrix = lastMatrix.getMatrix(0, last.getCa1Length() - 1, 0, 238 last.getCa2Length() - 1); 239 } 240 241 // Extract the structure identifier 242 optimalAFP = selfAlignments.get(0); 243 StructureIdentifier id = atoms[0].getGroup().getChain().getStructure() 244 .getStructureIdentifier(); 245 optimalAFP.setName1(id.getIdentifier()); 246 optimalAFP.setName2(id.getIdentifier()); 247 248 // Store the optimal self-alignment 249 result.setSelfAlignment(optimalAFP); 250 result.setStructureId(id); 251 252 // Determine the symmetry Type or get the one in params 253 SymmetryType type = params.getSymmType(); 254 if (type == SymmetryType.AUTO) { 255 if (result.getSelfAlignment().getBlockNum() == 1) { 256 type = SymmetryType.OPEN; 257 logger.info("Open Symmetry detected"); 258 } else { 259 type = SymmetryType.CLOSED; 260 logger.info("Close Symmetry detected"); 261 } 262 } 263 264 // Do not try the refinement if the self-alignment is not significant 265 if (optimalAFP.getTMScore() < params.getUnrefinedScoreThreshold()){ 266 result.setNumRepeats(1); 267 return result; 268 } 269 270 // STEP 3: order detection & symmetry refinement, apply consistency 271 try { 272 // ORDER DETECTION 273 OrderDetector orderDetector = null; 274 int order = 1; 275 switch (params.getOrderDetectorMethod()) { 276 case USER_INPUT: 277 order = params.getUserOrder(); 278 break; 279 case SEQUENCE_FUNCTION: 280 // Does not work for OPEN alignments 281 if (type == SymmetryType.CLOSED) { 282 orderDetector = new SequenceFunctionOrderDetector( 283 params.getMaxSymmOrder(), 0.4f); 284 order = orderDetector.calculateOrder( 285 result.getSelfAlignment(), atoms); 286 break; 287 } 288 case ANGLE: 289 // Does not work for OPEN alignments 290 if (type == SymmetryType.CLOSED) { 291 orderDetector = new AngleOrderDetectorPlus( 292 params.getMaxSymmOrder()); 293 order = orderDetector.calculateOrder( 294 result.getSelfAlignment(), atoms); 295 break; 296 } 297 case GRAPH_COMPONENT: 298 orderDetector = new GraphComponentOrderDetector(); 299 order = orderDetector.calculateOrder(result.getSelfAlignment(), 300 atoms); 301 break; 302 } 303 result.setNumRepeats(order); 304 305 // REFINEMENT 306 SymmetryRefiner refiner = null; 307 switch (params.getRefineMethod()) { 308 case NOT_REFINED: 309 return result; 310 case SEQUENCE_FUNCTION: 311 // Does not work for OPEN alignments 312 if (type == SymmetryType.CLOSED) { 313 refiner = new SequenceFunctionRefiner(); 314 break; 315 } 316 case GRAPH_COMPONENT: 317 refiner = new GraphComponentRefiner(); 318 break; 319 } 320 321 MultipleAlignment msa = refiner.refine(result.getSelfAlignment(), 322 atoms, order); 323 324 // Refinement succeeded, store results 325 result.setMultipleAlignment(msa); 326 result.setNumRepeats(msa.size()); 327 result.setRefined(true); 328 329 } catch (RefinerFailedException e) { 330 logger.info("Refinement failed: " + e.getMessage()); 331 return result; 332 } 333 334 // STEP 4: symmetry axes 335 SymmetryAxes axes = new SymmetryAxes(); 336 int order = result.getMultipleAlignment().size(); 337 Matrix4d axis = result.getMultipleAlignment().getBlockSet(0) 338 .getTransformations().get(1); 339 axes.addAxis(axis, order, type); 340 341 result.setAxes(axes); 342 return result; 343 } 344 345 /** 346 * Analyze the symmetries of the input Atom array using the DEFAULT 347 * parameters. 348 * 349 * @param atoms 350 * representative Atom array of the Structure 351 * @return CeSymmResult 352 * @throws StructureException 353 */ 354 public static CeSymmResult analyze(Atom[] atoms) throws StructureException { 355 CESymmParameters params = new CESymmParameters(); 356 return analyze(atoms, params); 357 } 358 359 /** 360 * Analyze the symmetries of the input Atom array using the provided 361 * parameters. 362 * 363 * @param atoms 364 * representative Atom array of the Structure 365 * @param param 366 * CeSymmParameters bean 367 * @return CeSymmResult 368 * @throws StructureException 369 */ 370 public static CeSymmResult analyze(Atom[] atoms, CESymmParameters params) 371 throws StructureException { 372 373 if (atoms.length < 1) 374 throw new IllegalArgumentException("Empty Atom array given."); 375 376 // If the SSE information is needed, we calculate it if the user did not 377 if (params.getSSEThreshold() > 0) { 378 Structure s = atoms[0].getGroup().getChain().getStructure(); 379 if (SecStrucTools.getSecStrucInfo(s).isEmpty()) { 380 logger.info("Calculating Secondary Structure..."); 381 SecStrucCalc ssp = new SecStrucCalc(); 382 ssp.calculate(s, true); 383 } 384 } 385 386 CeSymmIterative iter = new CeSymmIterative(params); 387 CeSymmResult result = iter.execute(atoms); 388 389 if (result.isRefined()) { 390 // Optimize the global alignment freely once more (final step) 391 if (params.getOptimization() && result.getSymmLevels() > 1) { 392 try { 393 SymmOptimizer optimizer = new SymmOptimizer(result); 394 MultipleAlignment optimized = optimizer.optimize(); 395 // Set the optimized MultipleAlignment and the axes 396 result.setMultipleAlignment(optimized); 397 } catch (RefinerFailedException e) { 398 logger.info("Final optimization failed:" + e.getMessage()); 399 } 400 } 401 result.getMultipleAlignment().getEnsemble() 402 .setStructureIdentifiers(result.getRepeatsID()); 403 } 404 return result; 405 } 406 407 /** 408 * Analyze a single level of symmetry. 409 * 410 * @param atoms 411 * Atom array of the current level 412 * @return CeSymmResult 413 * @throws StructureException 414 */ 415 public static CeSymmResult analyzeLevel(Atom[] atoms, 416 CESymmParameters params) throws StructureException { 417 418 if (atoms.length < 1) 419 throw new IllegalArgumentException("Empty Atom array given."); 420 421 CeSymmResult result = align(atoms, params); 422 423 if (result.isRefined()) { 424 // STEP 5: symmetry alignment optimization 425 if (result.getParams().getOptimization()) { 426 try { 427 MultipleAlignment msa = result.getMultipleAlignment(); 428 SymmOptimizer optimizer = new SymmOptimizer(result); 429 msa = optimizer.optimize(); 430 result.setMultipleAlignment(msa); 431 } catch (RefinerFailedException e) { 432 logger.debug("Optimization failed:" + e.getMessage()); 433 } 434 } 435 } 436 return result; 437 } 438 439}