Source code

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.fatcat.calc;
028
029import org.biojava.nbio.structure.Atom;
030import org.biojava.nbio.structure.StructureException;
031import org.biojava.nbio.structure.align.model.AFP;
032import org.biojava.nbio.structure.align.model.AFPChain;
033
034import java.util.List;
035
036public class AFPOptimizer
037{
038
039        public static final boolean debug = FatCatAligner.debug;
040
041        /**
042         * optimize the alignment by dynamic programming
043         */
044        public static void optimizeAln(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws StructureException
045        {
046
047                int minLen = afpChain.getMinLen();
048                int fragLen = params.getFragLen();
049
050
051                long optStart = System.currentTimeMillis();
052                int     i, a, k, p1, p2, bk, b1, b2, e1, e2, a1, a2;
053                int     iniLen;
054                int[][]     iniSet = new int[2][minLen];
055                int     maxi = 100;
056
057                int[][][] optAln = afpChain.getOptAln();
058                int[] optLen = afpChain.getOptLen();
059
060                int maxTra = params.getMaxTra();
061                double[] optRmsd = afpChain.getOptRmsd();
062
063                int blockNum = afpChain.getBlockNum();
064
065                if(optAln == null)      {
066                        optAln     = new int[maxTra+1][2][minLen];
067                        optLen     = new int[maxTra+1];
068                        afpChain.setOptLen(optLen);
069                        optRmsd    = new double[maxTra+1];
070                        afpChain.setOptRmsd(optRmsd);
071                }
072
073                List<AFP> afpSet = afpChain.getAfpSet();
074
075                int optLength         = afpChain.getOptLength();
076                int[] afpChainList    = afpChain.getAfpChainList();
077                int[] block2Afp       = afpChain.getBlock2Afp();
078                int[] blockSize       = afpChain.getBlockSize();
079
080
081                if (debug)
082                        System.out.println("AFPOptimizer got blockNum: " +blockNum);
083                //optimize each alignment defined by a block
084                b1 = b2 = e1 = e2 = optLength = 0;
085                for(bk = 0; bk < blockNum; bk ++)       {
086                        //initial aligned position
087                        iniLen = 0;
088                        if(bk > 0)      {
089                                b1 = e1;
090                                b2 = e2;
091                        }
092                        if(bk < blockNum - 1)   {
093                                a1 = afpChainList[block2Afp[bk] + blockSize[bk] - 1]; //the last AFP in current block
094                                a2 = afpChainList[block2Afp[bk + 1]];  //the first AFP in next block
095                                e1 = (afpSet.get(a1).getP1() + fragLen +  afpSet.get(a2).getP1()) / 2;
096                                e2 = (afpSet.get(a1).getP2() + fragLen +  afpSet.get(a2).getP2()) / 2;
097                        } //use the middle point of the current and next AFPs. old (starting point of next AFP)
098                        else    {
099                                e1 = ca1.length;
100                                e2 = ca2.length;
101                        }
102
103
104                        for(i = block2Afp[bk]; i < block2Afp[bk] + blockSize[bk]; i ++) {
105                                a = afpChainList[i];
106                                p1 = afpSet.get(a).getP1();
107                                p2 = afpSet.get(a).getP2();
108                                for(k = 0; k < afpSet.get(a).getFragLen(); k ++)     {
109                                        iniSet[0][iniLen] = p1 + k - b1; //note -b1
110                                        iniSet[1][iniLen] = p2 + k - b2; //note -b2
111                                        iniLen ++;
112
113                                }
114                        }
115                        //optimize the align by dynamic programming & constraint the optimization region
116                        if(debug) {
117                                System.err.println(String.format("optimize block %d (%d afp), region %d-%d(len %d), %d-%d(len %d)\n",
118                                                bk, blockSize[bk], b1, e1, e1-b1, b2, e2, e2-b2));
119
120                                System.err.println(" initial alignment Length: " + iniLen );
121                        }
122
123                        StructureAlignmentOptimizer opt = new StructureAlignmentOptimizer(b1,e1, ca1, b2,e2, ca2, iniLen, iniSet);
124                        opt.runOptimization(maxi);
125                        optRmsd[bk] = opt.optimizeResult(optLen,bk,optAln[bk]);
126
127                        //System.out.println(optRmsd[bk]);
128                        // SALNOPT *opt = new SALNOPT(e1-b1, &pro1->caCod[3 * b1], e2-b2, &pro2->caCod[3 * b2], iniLen, iniSet, maxi);
129                        // optRmsd[bk] = opt->OptimizeResult(&optLen[bk], optAln[bk]);
130
131                        if(debug)
132                                System.out.println(String.format(" optimized len=%d, rmsd %f\n", optLen[bk], optRmsd[bk]));
133
134                        for(i = 0; i < optLen[bk]; i ++)        {
135                                optAln[bk][0][i] += b1; //restore the position
136                                optAln[bk][1][i] += b2; //restore the position
137                        }
138                        optLength += optLen[bk];
139
140
141                }
142
143                long optEnd = System.currentTimeMillis();
144                if(debug)       System.out.println("complete AlignOpt " + (optEnd-optStart) +"\n");
145
146                if(optLength < minLen) {
147                        int[][][] optAln_trim = new int[maxTra + 1][2][optLength];
148                        for (i = 0; i < maxTra + 1; i ++) {
149                                System.arraycopy(optAln[i][0], 0, optAln_trim[i][0], 0, optLength);
150                                System.arraycopy(optAln[i][1], 0, optAln_trim[i][1], 0, optLength);
151                        }
152                        afpChain.setOptAln(optAln_trim);
153                } else {
154                        afpChain.setOptAln(optAln);
155                }
156
157                afpChain.setBlockNum(blockNum);
158                afpChain.setOptLength(optLength);
159                afpChain.setAfpChainList(afpChainList);
160                afpChain.setBlock2Afp(block2Afp);
161                afpChain.setBlockSize(blockSize);
162
163
164        }
165
166        /**
167         * get the afp list and residue list for each block
168         */
169
170        public static void blockInfo(AFPChain afpChain)
171        {
172                int     i, j, k, a, n;
173
174                int blockNum = afpChain.getBlockNum();
175
176                int[] blockSize =afpChain.getBlockSize();
177                int[] afpChainList = afpChain.getAfpChainList();
178                int[] block2Afp = afpChain.getBlock2Afp();
179                int[][][]blockResList = afpChain.getBlockResList();
180
181                List<AFP>afpSet = afpChain.getAfpSet();
182                int[] blockResSize = afpChain.getBlockResSize();
183
184                for(i = 0; i < blockNum; i ++)  {
185                        n = 0;
186                        for(j = 0; j < blockSize[i]; j ++)      {
187                                //the index in afpChainList, not in the whole afp set
188                                a = afpChainList[block2Afp[i] + j];
189                                for(k = 0; k < afpSet.get(a).getFragLen(); k ++)     {
190                                        blockResList[i][0][n] = afpSet.get(a).getP1() + k;
191                                        blockResList[i][1][n] = afpSet.get(a).getP2() + k;
192                                        n ++;
193                                }
194                        }
195                        blockResSize[i] = n;
196                }
197
198                afpChain.setBlockResSize(blockResSize);
199                afpChain.setBlockSize(blockSize);
200                afpChain.setAfpChainList(afpChainList);
201                afpChain.setBlock2Afp(block2Afp);
202                afpChain.setBlockResList(blockResList);
203        }
204
205        /**
206         * to update the chaining score after block delete and merge processed
207         * the blockScore value is important for significance evaluation
208         */
209        public static void updateScore(FatCatParameters params, AFPChain afpChain)
210        {
211                int     i, j, bknow, bkold, g1, g2;
212
213
214                afpChain.setConn(0d);
215                afpChain.setDVar(0d);
216
217                int blockNum = afpChain.getBlockNum();
218                int alignScoreUpdate = 0;
219                double[] blockScore = afpChain.getBlockScore();
220                int[] blockGap = afpChain.getBlockGap();
221                int[] blockSize =afpChain.getBlockSize();
222                int[] afpChainList = afpChain.getAfpChainList();
223                List<AFP>afpSet = afpChain.getAfpSet();
224                int[] block2Afp = afpChain.getBlock2Afp();
225
226                double torsionPenalty = params.getTorsionPenalty();
227
228
229                bkold = 0;
230                for(i = 0; i < blockNum; i ++)  {
231                        blockScore[i] = 0;
232                        blockGap[i] = 0;
233                        for(j = 0; j < blockSize[i]; j ++)      {
234                                bknow = afpChainList[block2Afp[i] + j];
235                                if(j == 0)      {
236                                        blockScore[i] = afpSet.get(bknow).getScore();
237                                }
238                                else    {
239                                        AFPChainer.afpPairConn(bkold, bknow, params, afpChain); //note: j, i
240                                        Double conn = afpChain.getConn();
241                                        blockScore[i] += afpSet.get(bknow).getScore() + conn;
242                                        g1 = afpSet.get(bknow).getP1() - afpSet.get(bkold).getP1() - afpSet.get(bkold).getFragLen();
243                                        g2 = afpSet.get(bknow).getP2() - afpSet.get(bkold).getP2() - afpSet.get(bkold).getFragLen();
244                                        blockGap[i] += (g1 > g2)?g1:g2;
245                                }
246                                bkold = bknow;
247                        }
248                        alignScoreUpdate += blockScore[i];
249                }
250                if(blockNum >= 2)       {
251                        alignScoreUpdate += (blockNum - 1) * torsionPenalty;
252                }
253
254                afpChain.setBlockGap(blockGap);
255                afpChain.setAlignScoreUpdate(alignScoreUpdate);
256                afpChain.setBlockScore(blockScore);
257                afpChain.setBlockSize(blockSize);
258                afpChain.setAfpChainList(afpChainList);
259                afpChain.setBlock2Afp(block2Afp);
260        }
261
262
263
264}