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.align.model.AFP;
031import org.biojava.nbio.structure.align.model.AFPChain;
032
033import java.util.List;
034
035/** does post processing after alignment chaingin
036 *
037 * @author Andreas Prlic
038 *
039 */
040public class AFPPostProcessor
041{
042
043        public static final boolean debug = FatCatAligner.debug;
044
045        public static void postProcess(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2){
046
047                int blockNum = afpChain.getBlockNum();
048                afpChain.setBlockNumIni(blockNum);
049                //PostProcess of chaining result
050
051                afpChain.setBlockNumIni(blockNum);
052
053                //split blocks (introduce twists) with high RMSD
054                splitBlock(params,afpChain, ca1,ca2);
055                blockNum = afpChain.getBlockNum();
056                afpChain.setBlockNumSpt( blockNum);
057
058                if ( debug){
059                        System.err.println("AFPPOstProcessor: postProcess blocknum = blocknumSpt:" + blockNum);
060                }
061
062                //redo: merge blocks with similar transformations & remove small blocks
063                //if(blockNum >= 2)     ClustBlock();
064
065                deleteBlock(params,afpChain,ca1,ca2);
066                mergeBlock(params,afpChain,ca1,ca2);
067
068                afpChain.setBlockNumClu(afpChain.getBlockNum());
069
070
071        }
072
073        /**
074         * in some special cases, there is no maginificent twists in the
075        final chaining result; however, their rmsd (original and after
076        optimizing) are very large. Therefore, a post-process is taken
077        to split the blocks further at the ralative bad connections (
078        with relative high distance variation)
079        to be tested:
080          split or not according to optimized or initial chaining???
081         */
082
083        private static void splitBlock(FatCatParameters params, AFPChain afpChain, Atom[] ca1, Atom[] ca2)
084        {
085                if ( debug)
086                        System.err.println("AFPPostProcessor: splitBlock");
087                int     i, a, bk, cut;
088                double  maxs, maxt;
089                int blockNum = afpChain.getBlockNum();
090                int maxTra = params.getMaxTra();
091                double badRmsd = params.getBadRmsd();
092
093                int     blockNum0 = blockNum;
094
095                double[] blockRmsd = afpChain.getBlockRmsd();
096                int[] blockSize = afpChain.getBlockSize();
097                int[] block2Afp = afpChain.getBlock2Afp();
098                double[] afpChainTwiList = afpChain.getAfpChainTwiList();
099
100                bk = 0;
101                while(blockNum < maxTra + 1)    {
102                        maxs = 0;
103                        for(i = 0; i < blockNum; i ++)   {
104                                if(blockRmsd[i] > maxs && blockSize[i] > 2) { //according to the optimized alignment
105                                        maxs = blockRmsd[i];
106                                        bk = i;
107                                } //!(Note: optRmsd, not blockRmsd, according to the optimized alignment
108                        }
109                        if(maxs < badRmsd)      break;
110                        maxt = 0;
111                        cut = 0;
112                        for(i = 1; i < blockSize[bk]; i ++)     {
113                                a = i + block2Afp[bk];
114                                if(afpChainTwiList[a] > maxt)   {
115                                        maxt = afpChainTwiList[a];
116                                        cut = i;
117
118                                }
119                        }
120                        if(debug)
121                                System.out.println(String.format("block %d original size %d rmsd %.3f maxt %.2f cut at %d\n", bk, blockSize[bk], maxs, maxt, cut));
122                        for(i = blockNum - 1; i > bk; i --)     {
123                                block2Afp[i + 1] = block2Afp[i];
124                                blockSize[i + 1] = blockSize[i];
125                                blockRmsd[i + 1] = blockRmsd[i];
126                        } //update block information
127                        block2Afp[bk + 1] = cut + block2Afp[bk];
128                        blockSize[bk + 1] = blockSize[bk] - cut;
129                        blockSize[bk] = cut;
130
131                        if(debug)
132                                System.out.println(String.format("  split into %d and %d sizes\n", blockSize[bk], blockSize[bk + 1]));
133
134
135                        int[] afpChainList = afpChain.getAfpChainList();
136                        //int[] subrange1    = getSubrange(afpChainList, block2Afp[bk + 1] );
137                        blockRmsd[bk + 1]  = AFPChainer.calAfpRmsd(blockSize[bk + 1],  afpChainList, block2Afp[bk + 1] , afpChain, ca1, ca2);
138
139                        //int[] subrange2    = getSubrange(afpChainList, block2Afp[bk] );
140                        blockRmsd[bk]      = AFPChainer.calAfpRmsd(blockSize[bk],      afpChainList, block2Afp[bk], afpChain, ca1, ca2);
141
142                        //split a block at the biggest position
143                        blockNum ++;
144                        afpChain.setAfpChainList(afpChainList);
145                }
146                if(blockNum - blockNum0 > 0)    {
147                        if(debug)
148                                System.out.println(String.format("Split %d times:\n", blockNum - blockNum0));
149                        for(i = 0; i < blockNum; i ++)  {
150                                if(debug)
151                                        System.out.println(String.format("  block %d size %d from %d rmsd %.3f\n", i, blockSize[i], block2Afp[i], blockRmsd[i]));
152                        }
153                }
154
155
156                afpChain.setBlockNum(blockNum);
157                afpChain.setBlockSize(blockSize);
158                afpChain.setBlockRmsd(blockRmsd);
159                afpChain.setBlock2Afp(block2Afp);
160
161
162        }
163
164        /**
165         * remove the artifical small rigid-body superimpose in the middle
166         clust the similar superimpositions (caused by the small flexible
167         region, which is detected as a seperate rigid superimposing region by adding
168         two twists before and after it(artifically!)
169         one possible solution: allowing long enough loops in the chaining process,
170         which however increase the calculation complexity
171         */
172        private static void deleteBlock(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2)
173        {
174                int blockNum = afpChain.getBlockNum();
175                List<AFP> afpSet = afpChain.getAfpSet();
176
177                int[] afpChainList =afpChain.getAfpChainList();
178
179
180
181                int[] block2Afp = afpChain.getBlock2Afp();
182                int[] blockSize = afpChain.getBlockSize();
183
184                double[] blockRmsd = afpChain.getBlockRmsd();
185
186                int fragLen = params.getFragLen();
187
188                //remove those blocks (both in terminals and in the middle) with only a AFP
189                //but still keep those small blocks spaning large regions
190                if(blockNum <= 1)       return;
191                int     blockNumOld = blockNum;
192                int     i, j, b1, b2, e1, e2, len;
193                e1 = e2 = 0;
194                for(i = 0; i < blockNum; i ++) {
195                        b1 = e1;
196                        b2 = e2;
197                        if(i < blockNum - 1)    {
198                                e1 = afpSet.get(afpChainList[block2Afp[i + 1]]).getP1();
199                                e2 = afpSet.get(afpChainList[block2Afp[i + 1]]).getP2();
200                        }
201                        else    {
202                                e1 = ca1.length;
203                                e2 = ca2.length;
204                        }
205                        if(blockSize[i] > 1)    continue;
206                        len = (e1 - b1) < (e2 - b2)?(e1 - b1):(e2 - b2);
207                        //if(i == blockNum - 1) blockNum --;
208                        if(len < 2 * fragLen)   {
209                                for(j = i; j < blockNum - 1; j ++)      {
210                                        blockRmsd[j] = blockRmsd[j + 1];
211                                        blockSize[j] = blockSize[j + 1];
212                                        block2Afp[j] = block2Afp[j + 1];
213                                }
214                                blockNum --;
215                                i --;
216                        } //delete a block
217                }
218                if(blockNumOld > blockNum)
219                        if(debug)
220                                System.out.println(
221                                                String.format("Delete %d small blocks\n", blockNumOld - blockNum)
222                                );
223
224
225                if (debug)
226                        System.err.println("deleteBlock: end blockNum:"+ blockNum);
227                afpChain.setBlock2Afp(block2Afp);
228                afpChain.setBlockSize(blockSize);
229                afpChain.setAfpChainList(afpChainList);
230                afpChain.setBlockNum(blockNum);
231                afpChain.setBlockRmsd(blockRmsd);
232        }
233
234
235        /**
236         * Merge consecutive blocks with similar transformation
237         */
238        private static  void mergeBlock(FatCatParameters params, AFPChain afpChain,Atom[] ca1,Atom[] ca2)
239        {
240
241                int blockNum = afpChain.getBlockNum();
242                double badRmsd = params.getBadRmsd();
243
244                int[] block2Afp = afpChain.getBlock2Afp();
245                int[] blockSize = afpChain.getBlockSize();
246
247                double[] blockRmsd = afpChain.getBlockRmsd();
248
249                int afpChainTwiNum = afpChain.getAfpChainTwiNum();
250
251                //clustering the neighbor blocks if their transformations are similar
252                int     i, j, b1, b2, minb1, minb2;
253                double  minrmsd;
254                int     merge = 0;
255                int     blockNumOld = blockNum;
256                double[][]  rmsdlist = null;
257                if(blockNum > 1)        {
258                        rmsdlist = new double[blockNumOld][blockNumOld];
259                        for(b1 = 0; b1 < blockNum - 1; b1 ++)   {
260                                for(b2 = b1 + 1; b2 < blockNum; b2 ++)  {
261                                        rmsdlist[b1][b2] = combineRmsd(b1, b2,afpChain,ca1,ca2);
262                                }
263                        }
264                }
265                minb1 = 0;
266                while(blockNum > 1)     {
267                        minrmsd = 1000;
268                        for(i = 0; i < blockNum - 1; i ++)      {
269                                j = i + 1; //only consider neighbor blocks
270                                if(minrmsd > rmsdlist[i][j])    {
271                                        minrmsd = rmsdlist[i][j];
272                                        minb1 = i;
273                                }
274                        }
275                        minb2 = minb1 + 1; //merge those most similar blocks
276                        //maxrmsd = (blockRmsd[minb1] > blockRmsd[minb2])?blockRmsd[minb1]:blockRmsd[minb2];
277                        if(minrmsd < badRmsd)   {
278                                if(debug)
279                                        System.out.println(String.format("merge block %d (rmsd %.3f) and %d (rmsd %.3f), total rmsd %.3f\n",
280                                                        minb1, blockRmsd[minb1], minb2, blockRmsd[minb2], minrmsd));
281                                blockSize[minb1] += blockSize[minb2];
282                                blockRmsd[minb1] = minrmsd;
283                                for(i = minb2; i < blockNum - 1; i ++)  {
284                                        block2Afp[i] = block2Afp[i + 1];
285                                        blockSize[i] = blockSize[i + 1];
286                                        blockRmsd[i] = blockRmsd[i + 1];
287                                } //update block information
288                                afpChainTwiNum --;
289                                blockNum --;
290                                for(b1 = 0; b1 < blockNum - 1; b1 ++)   {
291                                        for(b2 = b1 + 1; b2 < blockNum; b2 ++) {
292                                                if(b1 == minb1 || b2 == minb1)  {
293                                                        rmsdlist[b1][b2] = combineRmsd(b1, b2, afpChain,ca1,ca2);
294                                                }
295                                                else if(b2 < minb1)     continue;
296                                                else if(b1 < minb1)     {
297                                                        rmsdlist[b1][b2] = rmsdlist[b1][b2 + 1];
298                                                }
299                                                else    {
300                                                        rmsdlist[b1][b2] = rmsdlist[b1 + 1][b2 + 1];
301                                                }
302                                        }
303                                } //update the rmsdlist
304                                merge ++;
305                        } //merge two blocks
306                        else if(minrmsd >= 100) break;
307                        else    {
308                                rmsdlist[minb1][minb2] += 100;
309                        } //not merge, modify the rmsdlist so that this combination is not considered in next iteration
310                }
311
312                if(merge > 0)       {
313                        if(debug)
314                                System.out.println(String.format("Merge %d blocks, remaining %d blocks\n", merge, blockNum));
315                }
316
317                if (debug){
318                        System.err.println("AFPPostProcessor: mergeBlock end blocknum:" + blockNum);
319                }
320                afpChain.setBlock2Afp(block2Afp);
321                afpChain.setBlockSize(blockSize);
322                afpChain.setBlockNum(blockNum);
323                afpChain.setBlockRmsd(blockRmsd);
324                afpChain.setAfpChainTwiNum(afpChainTwiNum);
325        }
326
327
328        /**
329        return the rmsd of two blocks
330         */
331        private static double combineRmsd(int b1, int b2, AFPChain afpChain,Atom[] ca1,Atom[] ca2)
332        {
333                int     i;
334                int     afpn = 0;
335
336                int[] afpChainList =afpChain.getAfpChainList();
337
338                int[] block2Afp = afpChain.getBlock2Afp();
339                int[] blockSize = afpChain.getBlockSize();
340
341
342                int[]   list = new int[blockSize[b1]+blockSize[b2]];
343                for(i = block2Afp[b1]; i < block2Afp[b1] + blockSize[b1]; i ++) {
344                        list[afpn ++] = afpChainList[i];
345                }
346                for(i = block2Afp[b2]; i < block2Afp[b2] + blockSize[b2]; i ++) {
347                        list[afpn ++] = afpChainList[i];
348                }
349                double  rmsd = AFPChainer.calAfpRmsd(afpn, list,0, afpChain,ca1,ca2);
350
351                afpChain.setBlock2Afp(block2Afp);
352                afpChain.setBlockSize(blockSize);
353                afpChain.setAfpChainList(afpChainList);
354
355                return rmsd;
356        }
357
358
359}