Source code

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.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}