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}