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