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