001/* This class is based on the original FATCAT implementation by
002 * <pre>
003 * Yuzhen Ye & Adam Godzik (2003)
004 * Flexible structure alignment by chaining aligned fragment pairs allowing twists.
005 * Bioinformatics vol.19 suppl. 2. ii246-ii255.
006 * https://www.ncbi.nlm.nih.gov/pubmed/14534198
007 * </pre>
008 *
009 * Thanks to Yuzhen Ye and A. Godzik for granting permission to freely use and redistribute this code.
010 *
011 * This code may be freely distributed and modified under the
012 * terms of the GNU Lesser General Public Licence.  This should
013 * be distributed with the code.  If you do not have a copy,
014 * see:
015 *
016 *      http://www.gnu.org/copyleft/lesser.html
017 *
018 * Copyright for this code is held jointly by the individual
019 * authors.  These should be listed in @author doc comments.
020 *
021 *
022 * Created on Jun 17, 2009
023 * Created by Andreas Prlic - RCSB PDB
024 *
025 */
026
027package org.biojava.nbio.structure.align;
028
029import org.biojava.nbio.structure.*;
030import org.biojava.nbio.structure.align.model.AFP;
031import org.biojava.nbio.structure.align.model.AFPChain;
032import org.biojava.nbio.structure.geometry.Matrices;
033import org.biojava.nbio.structure.geometry.SuperPositions;
034import org.biojava.nbio.structure.jama.Matrix;
035import org.slf4j.Logger;
036import org.slf4j.LoggerFactory;
037
038import java.util.ArrayList;
039import java.util.List;
040
041import javax.vecmath.Matrix4d;
042
043//import org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol;
044
045public class AFPTwister {
046        private final static Logger logger = LoggerFactory
047                        .getLogger(AFPTwister.class);
048
049        // private static final boolean showAlignmentSteps = false;
050
051        /**
052         * calculate the total rmsd of the blocks output a merged pdb file for both
053         * proteins protein 1, in chain A protein 2 is twisted according to the
054         * twists detected, in chain B
055         *
056         * @return twisted Groups
057         */
058        public static Group[] twistPDB(AFPChain afpChain, Atom[] ca1, Atom[] ca2)
059                        throws StructureException {
060                // --------------------------------------------------------
061
062                if (afpChain.isShortAlign())
063                        return new Group[0];
064
065                List<AFP> afpSet = afpChain.getAfpSet();
066
067                int blockNum = afpChain.getBlockNum();
068
069                int i, b2, e2;
070
071                // superimposing according to the initial AFP-chaining
072                Atom[] origCA = StructureTools.cloneAtomArray(ca2);
073                Atom[] iniTwistPdb = StructureTools.cloneAtomArray(ca2);
074
075                int[] blockResSize = afpChain.getBlockResSize();
076
077                int[][][] blockResList = afpChain.getBlockResList();
078
079                int[] afpChainList = afpChain.getAfpChainList();
080                int[] block2Afp = afpChain.getBlock2Afp();
081                int[] blockSize = afpChain.getBlockSize();
082                int[] focusAfpList = afpChain.getFocusAfpList();
083
084                if (focusAfpList == null) {
085                        focusAfpList = new int[afpChain.getMinLen()];
086                        afpChain.setFocusAfpList(focusAfpList);
087                }
088
089                int focusAfpn = 0;
090                e2 = 0;
091                b2 = 0;
092
093                logger.debug("blockNUm at twister: ", blockNum);
094
095                for (int bk = 0; bk < blockNum; bk++) {
096
097                        // THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
098                        // copies the atoms over to iniTwistPdb later on in modifyCod
099
100                        transformOrigPDB(blockResSize[bk], blockResList[bk][0],
101                                        blockResList[bk][1], ca1, ca2, null, -1);
102
103                        // transform pro2 according to comparison of pro1 and pro2 at give
104                        // residues
105                        if (bk > 0) {
106                                b2 = e2;
107                        }
108
109                        logger.debug("b2 is {} before  modifyCon", b2);
110
111                        if (bk < (blockNum - 1)) {
112
113                                // bend at the middle of two consecutive AFPs
114                                int afpPos = afpChainList[block2Afp[bk] + blockSize[bk] - 1];
115                                AFP a1 = afpSet.get(afpPos);
116                                e2 = a1.getP2();
117
118                                int afpPos2 = afpChainList[block2Afp[bk + 1]];
119                                AFP a2 = afpSet.get(afpPos2);
120                                e2 = (a2.getP2() - e2) / 2 + e2;
121                                logger.debug("e2 : {}", e2);
122                        } else {
123                                // last one is until the end...
124                                e2 = ca2.length;
125                        }
126
127                        // this copies the coordinates over into iniTwistPdb
128                        cloneAtomRange(iniTwistPdb, ca2, b2, e2);
129
130                        // bound[bk] = e2;
131                        for (i = 0; i < blockSize[bk]; i++) {
132                                focusAfpList[focusAfpn] = afpChainList[block2Afp[bk] + i];
133                                focusAfpn++;
134                        }
135                }
136
137                int focusResn = afp2Res(afpChain, focusAfpn, focusAfpList, 0);
138
139                afpChain.setTotalLenIni(focusResn);
140
141                logger.debug(String.format("calrmsdini for %d residues", focusResn));
142
143                double totalRmsdIni = calCaRmsd(ca1, iniTwistPdb, focusResn,
144                                afpChain.getFocusRes1(), afpChain.getFocusRes2());
145
146                afpChain.setTotalRmsdIni(totalRmsdIni);
147                logger.debug("got iniRMSD: {}", totalRmsdIni);
148                if (totalRmsdIni == 5.76611141613097) {
149                        logger.debug("{}", afpChain.getAlnseq1());
150                        logger.debug("{}", afpChain.getAlnsymb());
151                        logger.debug("{}", afpChain.getAlnseq2());
152                }
153
154                afpChain.setFocusAfpList(focusAfpList);
155                afpChain.setBlock2Afp(block2Afp);
156                afpChain.setAfpChainList(afpChainList);
157
158                return twistOptimized(afpChain, ca1, origCA);
159        }
160
161        /**
162         * superimposing according to the optimized alignment
163         *
164         * @param afpChain
165         * @param ca1
166         * @param ca2
167         * @return Group array twisted.
168         * @throws StructureException
169         */
170        public static Group[] twistOptimized(AFPChain afpChain, Atom[] ca1,
171                        Atom[] ca2) throws StructureException {
172
173                Atom[] optTwistPdb = new Atom[ca2.length];
174
175                int gPos = -1;
176                for (Atom a : ca2) {
177                        gPos++;
178                        optTwistPdb[gPos] = a;
179                }
180
181                int blockNum = afpChain.getBlockNum();
182
183                int b2 = 0;
184                int e2 = 0;
185                int focusResn = 0;
186                int[] focusRes1 = afpChain.getFocusRes1();
187                int[] focusRes2 = afpChain.getFocusRes2();
188
189                if (focusRes1 == null) {
190                        focusRes1 = new int[afpChain.getCa1Length()];
191                        afpChain.setFocusRes1(focusRes1);
192                }
193                if (focusRes2 == null) {
194                        focusRes2 = new int[afpChain.getCa2Length()];
195                        afpChain.setFocusRes2(focusRes2);
196                }
197
198                int[] optLen = afpChain.getOptLen();
199                int[][][] optAln = afpChain.getOptAln();
200
201                for (int bk = 0; bk < blockNum; bk++) {
202                        // THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
203                        // copies the atoms over to iniTwistPdb later on in modifyCod
204                        transformOrigPDB(optLen[bk], optAln[bk][0], optAln[bk][1], ca1,
205                                        ca2, afpChain, bk);
206
207                        // transform pro2 according to comparison of pro1 and pro2 at give
208                        // residues
209                        if (bk > 0) {
210                                b2 = e2;
211                        }
212                        if (bk < blockNum - 1) { // bend at the middle of two consecutive
213                                                                                // blocks
214                                e2 = optAln[bk][1][optLen[bk] - 1];
215                                e2 = (optAln[bk + 1][1][0] - e2) / 2 + e2;
216                        } else {
217                                e2 = ca2.length;
218                        }
219                        cloneAtomRange(optTwistPdb, ca2, b2, e2);
220                        for (int i = 0; i < optLen[bk]; i++) {
221                                focusRes1[focusResn] = optAln[bk][0][i];
222                                focusRes2[focusResn] = optAln[bk][1][i];
223                                focusResn++;
224                        }
225                }
226                int totalLenOpt = focusResn;
227                logger.debug("calrmsdopt for {} residues", focusResn);
228                double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1,
229                                focusRes2);
230                logger.debug("got opt RMSD: {}", totalRmsdOpt);
231                int optLength = afpChain.getOptLength();
232
233                if (totalLenOpt != optLength) {
234                        logger.warn("Final alignment length is different {} {}",
235                                        totalLenOpt, optLength);
236                }
237                logger.debug("final alignment length {}, rmsd {}", focusResn,
238                                totalRmsdOpt);
239
240                afpChain.setTotalLenOpt(totalLenOpt);
241                afpChain.setTotalRmsdOpt(totalRmsdOpt);
242
243                return StructureTools.cloneGroups(optTwistPdb);
244
245        }
246
247        /**
248         * transform the coordinates in the ca2 according to the superimposing of
249         * the given position pairs. No Cloning, transforms input atoms.
250         */
251        // orig name: transPdb
252        private static void transformOrigPDB(int n, int[] res1, int[] res2,
253                        Atom[] ca1, Atom[] ca2, AFPChain afpChain, int blockNr)
254                        throws StructureException {
255                logger.debug(
256                                "transforming original coordinates {} len1: {} res1: {} len2: {} res2: {}",
257                                n, ca1.length, res1.length, ca2.length, res2.length);
258
259                Atom[] cod1 = getAtoms(ca1, res1, n, false);
260                Atom[] cod2 = getAtoms(ca2, res2, n, false);
261
262                // double *cod1 = pro1->Cod4Res(n, res1);
263                // double *cod2 = pro2->Cod4Res(n, res2);
264
265                Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(cod1),
266                                Calc.atomsToPoints(cod2));
267
268                Matrix r = Matrices.getRotationJAMA(transform);
269                Atom t = Calc.getTranslationVector(transform);
270
271                logger.debug("transPdb: transforming orig coordinates with matrix: {}",
272                                r);
273
274                if (afpChain != null) {
275                        Matrix[] ms = afpChain.getBlockRotationMatrix();
276                        if (ms == null)
277                                ms = new Matrix[afpChain.getBlockNum()];
278
279                        ms[blockNr] = r;
280
281                        Atom[] shifts = afpChain.getBlockShiftVector();
282                        if (shifts == null)
283                                shifts = new Atom[afpChain.getBlockNum()];
284                        shifts[blockNr] = t;
285
286                        afpChain.setBlockRotationMatrix(ms);
287                        afpChain.setBlockShiftVector(shifts);
288                }
289
290                for (Atom a : ca2)
291                        Calc.transform(a.getGroup(), transform);
292
293        }
294
295        // like Cod4Res
296        // most likely the clone flag is not needed
297        private static Atom[] getAtoms(Atom[] ca, int[] positions, int length,
298                        boolean clone) {
299
300                List<Atom> atoms = new ArrayList<Atom>();
301                for (int i = 0; i < length; i++) {
302                        int p = positions[i];
303                        Atom a;
304                        if (clone) {
305                                a = (Atom) ca[p].clone();
306                                a.setGroup((Group) ca[p].getGroup().clone());
307                        } else {
308                                a = ca[p];
309                        }
310                        atoms.add(a);
311                }
312                return atoms.toArray(new Atom[atoms.size()]);
313        }
314
315        /**
316         * Clones and copies the Atoms from p2 into p1 range is between r1 and r2
317         *
318         * @param p1
319         * @param p2
320         * @param r1
321         * @param r2
322         */
323        // orig name: modifyCod
324        private static void cloneAtomRange(Atom[] p1, Atom[] p2, int r1, int r2)
325                        throws StructureException {
326
327                logger.debug("modifyCod from: {} to: {}", r1, r2);
328
329                // special clone method, can;t use StructureTools.cloneCAArray, since we
330                // access the data
331                // slightly differently here.
332                List<Chain> model = new ArrayList<Chain>();
333                for (int i = r1; i < r2; i++) {
334
335                        Group g = p2[i].getGroup();
336                        Group newG = (Group) g.clone();
337
338                        p1[i] = newG.getAtom(p2[i].getName());
339                        Chain parentC = g.getChain();
340
341                        Chain newChain = null;
342
343                        for (Chain c : model) {
344                                if (c.getName().equals(parentC.getName())) {
345                                        newChain = c;
346                                        break;
347                                }
348                        }
349
350                        if (newChain == null) {
351                                newChain = new ChainImpl();
352                                newChain.setName(parentC.getName());
353                                model.add(newChain);
354                        }
355
356                        newChain.addGroup(newG);
357
358                } // modify caCod
359
360        }
361
362        /**
363         * Return the rmsd of the CAs between the input pro and this protein, at
364         * given positions. quite similar to transPdb but while that one transforms
365         * the whole ca2, this one only works on the res1 and res2 positions.
366         *
367         * Modifies the coordinates in the second set of Atoms (pro).
368         *
369         * @return rmsd of CAs
370         */
371        private static double calCaRmsd(Atom[] ca1, Atom[] pro, int resn,
372                        int[] res1, int[] res2) throws StructureException {
373
374                Atom[] cod1 = getAtoms(ca1, res1, resn, false);
375                Atom[] cod2 = getAtoms(pro, res2, resn, false);
376
377                if (cod1.length == 0 || cod2.length == 0) {
378                        logger.info("length of atoms  == 0!");
379                        return 99;
380                }
381
382                Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(cod1),
383                                Calc.atomsToPoints(cod2));
384
385                for (Atom a : cod2)
386                        Calc.transform(a.getGroup(), transform);
387
388                return Calc.rmsd(cod1, cod2);
389        }
390
391        /**
392         * Set the list of equivalent residues in the two proteins given a list of
393         * AFPs
394         *
395         * WARNING: changes the values for FocusRes1, focusRes2 and FocusResn in
396         * afpChain!
397         *
398         * @param afpChain
399         *            the AFPChain to store resuts
400         * @param afpn
401         *            nr of afp
402         * @param afpPositions
403         * @param listStart
404         * @return nr of eq residues
405         */
406
407        public static int afp2Res(AFPChain afpChain, int afpn, int[] afpPositions,
408                        int listStart) {
409                int[] res1 = afpChain.getFocusRes1();
410                int[] res2 = afpChain.getFocusRes2();
411                int minLen = afpChain.getMinLen();
412
413                int n = 0;
414
415                List<AFP> afpSet = afpChain.getAfpSet();
416
417                for (int i = listStart; i < listStart + afpn; i++) {
418                        int a = afpPositions[i];
419                        for (int j = 0; j < afpSet.get(a).getFragLen(); j++) {
420                                if (n >= minLen) {
421                                        throw new RuntimeException(
422                                                        "Error: too many residues in AFPChainer.afp2Res!");
423                                }
424                                res1[n] = afpSet.get(a).getP1() + j;
425                                res2[n] = afpSet.get(a).getP2() + j;
426                                n++;
427                        }
428                }
429
430                afpChain.setFocusRes1(res1);
431                afpChain.setFocusRes2(res2);
432                afpChain.setFocusResn(n);
433
434                if (n == 0) {
435                        logger.warn("n=0!!! + " + afpn + " " + listStart + " "
436                                        + afpPositions.length);
437                }
438                return n;
439        }
440
441}