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