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