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.Arrays; 025import java.util.Collections; 026import java.util.List; 027import java.util.Set; 028 029import javax.vecmath.Matrix4d; 030 031import org.biojava.nbio.structure.Atom; 032import org.biojava.nbio.structure.StructureException; 033import org.biojava.nbio.structure.StructureIdentifier; 034import org.biojava.nbio.structure.align.multiple.Block; 035import org.biojava.nbio.structure.align.multiple.BlockImpl; 036import org.biojava.nbio.structure.align.multiple.BlockSet; 037import org.biojava.nbio.structure.align.multiple.BlockSetImpl; 038import org.biojava.nbio.structure.align.multiple.MultipleAlignment; 039import org.biojava.nbio.structure.align.multiple.MultipleAlignmentImpl; 040import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentScorer; 041import org.biojava.nbio.structure.secstruc.SecStrucElement; 042import org.biojava.nbio.structure.secstruc.SecStrucTools; 043import org.biojava.nbio.structure.secstruc.SecStrucType; 044import org.biojava.nbio.structure.symmetry.internal.CESymmParameters.RefineMethod; 045import org.biojava.nbio.structure.symmetry.internal.CESymmParameters.SymmetryType; 046import org.biojava.nbio.structure.symmetry.utils.SymmetryTools; 047import org.jgrapht.Graph; 048import org.jgrapht.alg.connectivity.ConnectivityInspector; 049import org.jgrapht.graph.DefaultEdge; 050import org.jgrapht.graph.SimpleGraph; 051import org.slf4j.Logger; 052import org.slf4j.LoggerFactory; 053 054/** 055 * Iterative version of CeSymm that aims at identifying all symmetry axis of a 056 * structure. 057 * <p> 058 * Works in the following way: 059 * <ul> 060 * <li>Run CeSymm on the original structure. 061 * <li>Calculate the symmetric unit boundaries. 062 * <li>Run CeSymm on one of the symmetric units to find further symmetries. 063 * <li>Repeat the last two steps until no more significant results are found. 064 * <li>Map back all residues in a multiple alignment of the repeats. 065 * </ul> 066 * </li> 067 * 068 * @author Aleix Lafita 069 * @since 4.1.1 070 * 071 */ 072public class CeSymmIterative { 073 074 private static final Logger logger = LoggerFactory 075 .getLogger(CeSymmIterative.class); 076 077 private CESymmParameters params; 078 private Graph<Integer, DefaultEdge> alignGraph; // cumulative 079 private List<CeSymmResult> levels; // symmetry at each level 080 081 /** 082 * For the iterative algorithm to work properly the refinement and 083 * optimization options should be turned on, because the alignment has to be 084 * consistent at every recursive step. 085 * 086 * @param param 087 * CeSymm parameters, make sure they are cloned 088 */ 089 public CeSymmIterative(CESymmParameters param) { 090 params = param; 091 alignGraph = new SimpleGraph<Integer, DefaultEdge>(DefaultEdge.class); 092 levels = new ArrayList<CeSymmResult>(); 093 } 094 095 /** 096 * This method uses iteratively CeSymm to calculate all symmetries in the 097 * input array of atoms and organize them in a multiple alignment of the 098 * repeats. 099 * 100 * @param atoms 101 * atoms 102 * @return CeSymmResult 103 * 104 * @throws StructureException 105 */ 106 public CeSymmResult execute(Atom[] atoms) throws StructureException { 107 108 // First iterate through all levels and then reconstruct all repeats 109 iterate(atoms); 110 return reconstructSymmResult(atoms); 111 112 } 113 114 /** 115 * This method runs iteratively the analysis of one level of symmetry with 116 * CeSymm on the input Atom array until no more symmetries exist. 117 * 118 * @param atoms 119 * representative Atom array of the Structure 120 * @return true if any symmetry was found, false if asymmetric 121 * @throws StructureException 122 */ 123 private void iterate(Atom[] atoms) throws StructureException { 124 125 logger.debug("Starting new iteration..."); 126 127 // Return if the Atom array is too short 128 if ((atoms.length <= params.getWinSize() 129 || atoms.length <= params.getMinCoreLength()) 130 && !levels.isEmpty()) { 131 logger.debug("Aborting iteration due to insufficient Atom " 132 + "array length: %d", atoms.length); 133 return; 134 } 135 136 // Return if the maximum levels of symmetry have been reached 137 if (params.getSymmLevels() > 0) { 138 if (levels.size() == params.getSymmLevels()) 139 return; 140 } 141 142 // Perform one level CeSymm alignment 143 CeSymmResult result = CeSymm.analyzeLevel(atoms, params); 144 145 if (params.getRefineMethod() == RefineMethod.NOT_REFINED 146 || !result.isSignificant()) { 147 if (levels.isEmpty()) 148 levels.add(result); 149 return; 150 } 151 152 // Generate the Atoms of one of the symmetric repeat 153 Integer start = null; 154 int it = 0; 155 while (start == null) { 156 start = result.getMultipleAlignment().getBlocks().get(0) 157 .getAlignRes().get(0).get(it); 158 it++; 159 } 160 Integer end = null; 161 it = result.getMultipleAlignment().getBlocks().get(0).getAlignRes() 162 .get(0).size() - 1; 163 while (end == null) { 164 end = result.getMultipleAlignment().getBlocks().get(0) 165 .getAlignRes().get(0).get(it); 166 it--; 167 } 168 Atom[] atomsR = Arrays.copyOfRange(atoms, start, end + 1); 169 170 // Check the SSE requirement 171 if (countHelixStrandSSE(atomsR) < params.getSSEThreshold()) { 172 if (levels.isEmpty()) 173 levels.add(result); 174 return; 175 } 176 177 // If symmetric store the residue dependencies in alignment graph 178 Block b = result.getMultipleAlignment().getBlock(0); 179 for (int pos = 0; pos < b.length(); pos++) { 180 for (int su = 0; su < b.size() - 1; su++) { 181 Integer pos1 = b.getAlignRes().get(su).get(pos); 182 Integer pos2 = b.getAlignRes().get(su + 1).get(pos); 183 // Add edge from lower to higher positions 184 if (pos1 != null && pos2 != null) { 185 alignGraph.addVertex(pos1); 186 alignGraph.addVertex(pos2); 187 alignGraph.addEdge(pos1, pos2); 188 } 189 } 190 } 191 192 // Iterate further on those Atoms (of the first repeat only) 193 levels.add(result); 194 iterate(atomsR); 195 } 196 197 /** 198 * After all the analysis iterations have finished, the final Result object 199 * is reconstructed using the cumulative alignment graph. 200 * 201 * @param atoms 202 * the original structure atoms 203 * @return CeSymmResult reconstructed symmetry result 204 * @throws StructureException 205 */ 206 private CeSymmResult reconstructSymmResult(Atom[] atoms) 207 throws StructureException { 208 209 // If one level, nothing to build or calculate 210 if (levels.size() == 1) 211 return levels.get(0); 212 213 CeSymmResult result = new CeSymmResult(); 214 result.setSelfAlignment(levels.get(0).getSelfAlignment()); 215 result.setStructureId(levels.get(0).getStructureId()); 216 result.setAtoms(levels.get(0).getAtoms()); 217 result.setParams(levels.get(0).getParams()); 218 219 // Initialize a new multiple alignment 220 MultipleAlignment msa = new MultipleAlignmentImpl(); 221 msa.getEnsemble().setAtomArrays(new ArrayList<Atom[]>()); 222 msa.getEnsemble().setStructureIdentifiers( 223 new ArrayList<StructureIdentifier>()); 224 msa.getEnsemble().setAlgorithmName(CeSymm.algorithmName); 225 msa.getEnsemble().setVersion(CeSymm.version); 226 227 BlockSet bs = new BlockSetImpl(msa); 228 Block b = new BlockImpl(bs); 229 b.setAlignRes(new ArrayList<List<Integer>>()); 230 231 // Calculate the connected groups of the alignment graph 232 ConnectivityInspector<Integer, DefaultEdge> inspector = new ConnectivityInspector<Integer, DefaultEdge>( 233 alignGraph); 234 List<Set<Integer>> comps = inspector.connectedSets(); 235 List<ResidueGroup> groups = new ArrayList<ResidueGroup>(comps.size()); 236 for (Set<Integer> comp : comps) 237 groups.add(new ResidueGroup(comp)); 238 239 // Calculate the total number of repeats 240 int order = 1; 241 for (CeSymmResult sr : levels) 242 order *= sr.getMultipleAlignment().size(); 243 for (int su = 0; su < order; su++) 244 b.getAlignRes().add(new ArrayList<Integer>()); 245 246 // Construct the resulting MultipleAlignment from ResidueGroups 247 for (ResidueGroup group : groups) { 248 if (group.order() != order) 249 continue; 250 group.combineWith(b.getAlignRes()); 251 } 252 253 // The reconstruction failed, so the top level is returned 254 if (b.length() == 0) 255 return levels.get(0); 256 257 for (int su = 0; su < order; su++) { 258 Collections.sort(b.getAlignRes().get(su)); 259 msa.getEnsemble().getAtomArrays().add(atoms); 260 msa.getEnsemble().getStructureIdentifiers() 261 .add(result.getStructureId()); 262 } 263 264 result.setMultipleAlignment(msa); 265 result.setRefined(true); 266 result.setNumRepeats(order); 267 268 SymmetryAxes axes = recoverAxes(result); 269 result.setAxes(axes); 270 271 // Set the transformations and scores of the final alignment 272 SymmetryTools 273 .updateSymmetryTransformation(result.getAxes(), msa); 274 double tmScore = MultipleAlignmentScorer.getAvgTMScore(msa) 275 * msa.size(); 276 double rmsd = MultipleAlignmentScorer.getRMSD(msa); 277 msa.putScore(MultipleAlignmentScorer.AVGTM_SCORE, tmScore); 278 msa.putScore(MultipleAlignmentScorer.RMSD, rmsd); 279 280 return result; 281 } 282 283 /** 284 * The symmetry axes of each level are recovered after the symmetry analysis 285 * iterations have finished, using the stored MultipleAlignment at each 286 * symmetry level. 287 * @return SymmetryAxes 288 */ 289 private SymmetryAxes recoverAxes(CeSymmResult result) { 290 291 SymmetryAxes axes = new SymmetryAxes(); 292 293 for (int m = 0; m < levels.size(); m++) { 294 295 MultipleAlignment align = levels.get(m).getMultipleAlignment(); 296 Matrix4d axis = align.getBlockSet(0).getTransformations().get(1); 297 SymmetryType type = levels.get(m).getAxes().getElementaryAxis(0).getSymmType(); 298 int order = align.size(); 299 300 axes.addAxis(axis, order, type); 301 } 302 return axes; 303 } 304 305 /** 306 * Calculate the number of helix and strand SSE of a repeat. 307 * 308 * @param atoms 309 * Atom array of the repeat found 310 * @return int number of helix or strand SSE 311 */ 312 private static int countHelixStrandSSE(Atom[] atoms) { 313 314 List<SecStrucElement> sses = SecStrucTools 315 .getSecStrucElements(SymmetryTools.getGroups(atoms)); 316 int count = 0; 317 318 //keep track of different helix types 319 boolean helix = false; 320 int hEnd = 0; 321 322 for (SecStrucElement sse : sses) { 323 SecStrucType t = sse.getType(); 324 if (t.isBetaStrand()) { 325 helix = false; 326 count++; 327 } else if (t.isHelixType()){ 328 if (helix){ 329 // If this helix is contiguous to the previous 330 if (sse.getRange().getStart().getSeqNum() + 1 == hEnd) 331 hEnd = sse.getRange().getEnd().getSeqNum(); 332 else 333 count++; 334 } else 335 count++; 336 } else 337 helix = false; 338 } 339 return count; 340 } 341 342}