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.align.multiple.util;
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.Calc;
030import org.biojava.nbio.structure.StructureException;
031import org.biojava.nbio.structure.StructureTools;
032import org.biojava.nbio.structure.align.multiple.Block;
033import org.biojava.nbio.structure.align.multiple.BlockSet;
034import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
035import org.biojava.nbio.structure.geometry.SuperPositionSVD;
036import org.biojava.nbio.structure.geometry.SuperPositions;
037
038/**
039 * Superimposes the core aligned residues of every structure in a
040 * {@link MultipleAlignment} onto a reference structure. This method
041 * can eliminate the pairwise similarities of some structures to the
042 * reference, when doing the superposition, taking into account only
043 * those shared parts between the structures.
044 * <p>
045 * Performs a global superposition of the MultipleAlignment in case
046 * there is only one {@link BlockSet}, and a superposition for every
047 * BlockSet in case there is more than one (flexible alignment).
048 * <p>
049 * This class uses the {@link SuperPositionSVD} algorithm.
050 *
051 * @author Aleix Lafita
052 * @since 4.2.0
053 *
054 */
055public class CoreSuperimposer implements MultipleSuperimposer {
056
057        private int reference;
058
059        /**
060         * Default Constructor.
061         * Uses the first structure as the reference.
062         */
063        public CoreSuperimposer() {
064                this(0);
065        }
066
067        /**
068         * Constructor using a specified structure as reference.
069         *
070         * @param reference Index of the structure to use as a reference
071         *                      (it has to be > 0)
072         */
073        public CoreSuperimposer(int reference) {
074                if (reference<0) {
075                        throw new IllegalArgumentException(
076                                        "reference index has to be positive, but was "+reference);
077                }
078                this.reference = reference;
079        }
080
081        @Override
082        public void superimpose(MultipleAlignment alignment)
083                        throws StructureException {
084
085                //Check for inconsistencies in the alignment
086                if(alignment.getEnsemble() == null) {
087                        throw new NullPointerException("No ensemble set for this alignment."
088                                        + " Structure information cannot be obtained.");
089                }
090                if (alignment.size() < 1) {
091                        throw new IndexOutOfBoundsException(
092                                        "No aligned structures, alignment size == 0.");
093                }
094                if (alignment.getCoreLength() < 1){
095                        throw new IndexOutOfBoundsException(
096                                        "Alignment too short, core alignment length < 1.");
097                }
098
099                List<Atom[]> atomArrays = alignment.getAtomArrays();
100                if (atomArrays.size() <= reference) {
101                        throw new IndexOutOfBoundsException(String.format(
102                                        "Invalid reference structure: requested %d but "
103                                                        + "only %d structures.",
104                                                        reference,atomArrays.size()));
105                }
106
107                alignment.clear();
108
109                //Calculate BlockSet transformations
110                for (BlockSet bs:alignment.getBlockSets()){
111
112                        //Block transformations
113                        List<Matrix4d> transforms =
114                                        new ArrayList<Matrix4d>(atomArrays.size());
115
116                        //Loop through structures
117                        for (int i=0; i<atomArrays.size(); i++){
118
119                                if( i == reference) {
120                                        //Identity operation
121                                        Matrix4d ident = new Matrix4d();
122                                        ident.setIdentity();
123                                        transforms.add(ident);
124                                        continue;
125                                }
126
127                                Atom[] ref = atomArrays.get(reference);
128                                Atom[] curr = atomArrays.get(i);
129
130                                List<Atom> atomSet1 = new ArrayList<Atom>();
131                                List<Atom> atomSet2 = new ArrayList<Atom>();
132
133                                for( Block blk : bs.getBlocks() ) {
134
135                                        if( blk.size() != atomArrays.size()) {
136                                                throw new IllegalStateException(String.format(
137                                                                "Mismatched block length. Expected %d "
138                                                                                + "structures, found %d.",
139                                                                                atomArrays.size(),blk.size() ));
140                                        }
141
142                                        List<Integer> corePositions =
143                                                        MultipleAlignmentTools.getCorePositions(blk);
144
145                                        //Loop through all aligned residues
146                                        for (int j=0; j<blk.length(); j++){
147                                                //Check that the position is in the core
148                                                if (!corePositions.contains(j)) continue;
149
150                                                Integer pos1 = blk.getAlignRes().get(reference).get(j);
151                                                Integer pos2 = blk.getAlignRes().get(i).get(j);
152
153                                                if (pos1==null || pos2==null) continue;
154
155                                                atomSet1.add(ref[pos1]);
156                                                atomSet2.add(curr[pos2]);
157                                        }
158                                }
159                                Atom[] array1 = atomSet1.toArray(new Atom[atomSet1.size()]);
160                                Atom[] array2 = atomSet2.toArray(new Atom[atomSet2.size()]);
161
162                                array2 = StructureTools.cloneAtomArray(array2);
163
164                                //From the superimposer we obtain the rotation and translation
165                                Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(array1),
166                                                Calc.atomsToPoints(array2));
167                                transforms.add(trans);
168                        }
169                        //Set transformation of the BlockSet
170                        bs.setTransformations(transforms);
171                }
172        }
173}