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}