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.Calc; 031import org.biojava.nbio.structure.StructureException; 032import org.biojava.nbio.structure.align.AFPTwister; 033import org.biojava.nbio.structure.align.model.AFP; 034import org.biojava.nbio.structure.align.model.AFPChain; 035import org.biojava.nbio.structure.geometry.SuperPositions; 036import org.biojava.nbio.structure.jama.Matrix; 037 038import java.util.List; 039 040import javax.vecmath.Matrix4d; 041 042/** a class to chain AFPs to an alignment 043 * 044 * @author Andreas Prlic 045 * 046 */ 047public class AFPChainer 048{ 049 public static final boolean debug = FatCatAligner.debug; 050 // private static final boolean showAlig = false; 051 052 /** 053 // Key function: chain (assembly) the AFPs 054 // a AFP (k) is defined as (i, j, k), with i and j are staring points 055 // AFP extension (eg. AFP(k-1) -> AFP(k) ) requiring 056 // AFP(k-1) < AFP(k)(refer AFP.h definition), 057 // ie i(k-1) < i(k) and j(k-1) < j(k) 058 // in the figure, only (2) AFP can extend to AFP(k) 059 // Key features: a coordination transformation is allowed in the AFP extension 060 // gap penalties are also considered 061 // 062 // protein1 063 // --------------------------- 064 // | \ | 065 // | \(1) | 066 // | \ \ | 067 // | \(2) \ | 068 // p | \ | 069 // r | \ | 070 // o | \(3) \(i,j, k) | 071 // t | \ \ | 072 // e | \ | 073 // i | | 074 // n | | 075 // 2 --------------------------- 076 // schematic of AFP chaining 077 */ 078 public static void doChainAfp(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2){ 079 List<AFP> afpSet = afpChain.getAfpSet(); 080 081 afpChain.setConn(0d); 082 afpChain.setDVar(0d); 083 084 int afpNum = afpSet.size(); 085 if(afpNum <= 0) return; 086 087 088 int[] twi = afpChain.getTwi(); 089 if(twi == null) { 090 twi = new int[afpNum]; 091 afpChain.setTwi(twi); 092 } 093 094 //transformation, calculated at DoChainAfp, be used in List extraction 095 096 //forward: calculate the score matrix 097 boolean isConnected = false; 098 int i, j, j0, n; 099 double stmp; 100 101 afpChain.setConn(0d); 102 afpChain.setDVar(0d); 103 104 double[] sco = new double[afpNum]; //the score ending at an AFP 105 int[] pre = new int[afpNum]; //the previous AFP 106 double maxsco = 0; 107 int maxafp = 0; 108 int[] list = new int[afpNum]; 109 110 int maxGap = params.getMaxGap(); 111 int fragLen = params.getFragLen(); 112 int maxTra = params.getMaxTra(); 113 114 115 Matrix disTable1 = getDisTable(maxGap + 2 * fragLen + 1,ca1); 116 Matrix disTable2 = getDisTable(maxGap + 2 * fragLen + 1,ca2); 117 118 afpChain.setDisTable1(disTable1); 119 afpChain.setDisTable2(disTable2); 120 121 for(i = 0; i < afpNum; i ++) { 122 sco[i] = afpSet.get(i).getScore(); //start from itself 123 pre[i] = -1; 124 twi[i] = 0; 125 if ( afpSet.get(i).getP1() < fragLen || afpSet.get(i).getP2() < fragLen) 126 n = 0; 127 else 128 n = getCompatibleAfps(i, list, params, afpChain); //get a compatible list 129 //printf("afp %d, compatible %d\n", i, n); 130 for(j0 = 0; j0 < n; j0 ++) { 131 j = list[j0]; 132 isConnected = afpPairConn(j, i, params,afpChain); //note: j, i 133 Double conn = afpChain.getConn(); 134 int t = 0; 135 if ( isConnected) 136 t=1; 137 if(twi[j] + t > maxTra) continue; 138 //two many transformation are disfavored 139 stmp = sco[j] + afpSet.get(i).getScore() + conn; 140 if(stmp > sco[i]) { //considered all previous compatible AFPs 141 sco[i] = stmp; 142 twi[i] = twi[j] + t; 143 pre[i] = j; 144 } 145 } 146 if(maxsco < sco[i]) { 147 maxsco = sco[i]; 148 maxafp = i; 149 } 150 } 151 152 int currafp = maxafp; 153 if(debug) 154 System.out.printf("maximum score %f, %d%n%n", maxsco, twi[currafp]); 155 156 //trace-back from maxafp (maxsco) 157 158 afpChain.setAlignScore(maxsco); 159 afpChain.setAlignScoreUpdate(maxsco); 160 afpChain.setAfpChainTwiNum(0); 161 162 traceBack(pre, currafp, twi[currafp],params,afpChain,ca1,ca2); 163 164 165 } 166 167 private static Matrix getDisTable(int maxlen, Atom[]ca) 168 169 { 170 int length = ca.length; 171 Matrix dis = new Matrix(length,length); 172 173 int i, j; 174 for(i = 0; i < length; i ++) { 175 dis.set(i,i,0); 176 for(j = i + 1;( j < length) && (j <= i + maxlen); j ++) { 177 dis.set(i,j,0); 178 179 double val = dis.get(i,j) + (Calc.getDistance(ca[i],ca[j])) * (Calc.getDistance(ca[i],ca[j])); 180 dis.set(i,j,val); 181 182 183 dis.set(i,j,Math.sqrt(dis.get(i,j))); 184 dis.set(j,i,dis.get(i,j)); 185 } 186 } 187 return dis; 188 189 } 190 191 /* 192 193 derive the compabitle AFP lists for AFP-chaining 194 this is important for speeding up the process 195 for a given AFP(i1,j1), there are three regions that could be the starting 196 point for the compabitle AFPs of AFP(i1,j1) 197 // a1 a2 a3 198 // i1-G i1-f-c i1-f+1 i1 199 // | | | | 200 // ---------------------------- 201 // | B | | 202 //b1 j1-G -| ---------------| | 203 // | | | | | 204 // | | C | 3 | | 205 // | | | | | 206 //b2 j1-f-c -| |--------------| | 207 // | | 2 | 1 | | 208 //b3 j1-f+1 -|------------------ | 209 // | A | 210 // j1 -|---------------------\ 211 // | \ (AFP(i1,j1)) 212 // ----------------------------- 213 // 214 f: the length of AFPs (we use segments of same length) 215 G: g + f, where g is the maximum allowed gaps 216 c: the maximum allowed cross-over in AFP-connection, 217 here we use c = f, and j1-f-c = j1-2f 218 incompatible region A: its AFPs overlap with given AFP(i1,j1) 219 incompatible region B: the gaps between its AFP with given AFP is larger than g 220 incompatible region C: its AFPs connect with given AFP but cross a given threshold. 221 compatible region 1: [i1-f-c,i1-f+1>,[j1-f-c,j1-f+1> or [a2,a3],[b2,b3] 222 compatible region 2: [i1-G,i1-f-c],[j1-f-c,j1-f] or [a1,a2],[b2,b3] 223 combine 1 and 2 : [i1-G,i1-f],[j1-f-c,j1-f] or [a1,a3],[b2,b3] 224 compatible region 3: [i1-f-c,i1-f],[j1-G,j1-f-c] or [a2,a3],[b1,b2] 225 c->misCut 226 f->fragLen 227 G->fragLen+maxGap->maxGapFrag 228 * 229 * 230 */ 231 private static int getCompatibleAfps(int afp, int[] list, FatCatParameters params, AFPChain afpChain){ 232 233 int i, j, i1, j1, f, G, c, a1, a2, a3, b1, b2, b3, s1, s2; 234 235 int fragLen = params.getFragLen(); 236 int maxGapFrag = params.getMaxGapFrag(); 237 int misCut = params.getMisCut(); 238 int maxTra = params.getMaxTra(); 239 List<AFP> afpSet = afpChain.getAfpSet(); 240 241 f = fragLen; 242 G = maxGapFrag; 243 c = misCut; 244 245 i1 = afpSet.get(afp).getP1(); 246 j1 = afpSet.get(afp).getP2(); 247 a3 = i1 - f; 248 a2 = a3 - c; 249 a1 = i1 - G; 250 a2 = a2>0?a2:0; 251 a1 = a1>0?a1:0; 252 253 b3 = j1 - f; 254 b2 = b3 - c; 255 b1 = j1 - G; 256 b2 = (b2 > 0)?b2:0; 257 b1 = (b1 > 0)?b1:0; 258 259 int[][] afpAftIndex = afpChain.getAfpAftIndex(); 260 int[][] afpBefIndex = afpChain.getAfpBefIndex(); 261 int[] twi = afpChain.getTwi(); 262 263 264 265 int n = 0; 266 //compatible region 1-2, [a1,a3][b2,b3] 267 for(i = a1; i <= a3; i ++) {//i <= a3 instead of i < a3 268 s1 = afpAftIndex[i][b2]; //note afpAftIndex, not afpIndex 269 if(s1 < 0) continue;//no AFP for the given i with j > b2 270 s2 = afpBefIndex[i][b3]; //afps is sorted by j given a i,it's sparse matrix 271 if(s2 < 0) continue;//no AFP for the given i with j < b3 272 for(j = s1; j <= s2; j ++) { //j <= s2 instead of j < s2 273 if(twi[j] <= maxTra) { 274 list[n ++] = j; 275 } 276 } 277 } 278 279 //compatible region 3 [a2,a3][b1,b2] 280 for(i = a2; i <= a3; i ++) { 281 s1 = afpAftIndex[i][b1]; 282 if(s1 < 0) continue; 283 s2 = afpBefIndex[i][b2]; //afps is sorted by j given a i 284 if(s2 < 0) continue; 285 //note j < s2, as the cases of j == s2 is alread considered in previous region 286 for(j = s1; j < s2; j ++) { 287 if(twi[j] <= maxTra) { 288 list[n ++] = j; 289 } 290 } 291 } 292 293 return n; 294 295 } 296 297 298 /** 299 //Key function: calculate the connectivity of AFP pairs 300 //no compatibility criteria is executed 301 //note: afp1 is previous to afp2 in terms of the position 302 //this module must be optimized 303 * 304 * @param afp1 305 * @param afp2 306 * @return flag if they are connected 307 */ 308 public static boolean afpPairConn(int afp1, int afp2, FatCatParameters params, AFPChain afpChain) 309 310 { 311 312 Double conn = afpChain.getConn(); 313 Double dvar = afpChain.getDVar(); 314 315 double misScore = params.getMisScore(); 316 double maxPenalty = params.getMaxPenalty(); 317 double disCut = params.getDisCut(); 318 double gapExtend = params.getGapExtend(); 319 double torsionPenalty = params.getTorsionPenalty(); 320 double disSmooth = params.getDisSmooth(); 321 322 List<AFP> afpSet = afpChain.getAfpSet(); 323 324 int m = calcGap(afpSet.get(afp2),afpSet.get(afp1)); 325 int g = calcMismatch(afpSet.get(afp2),afpSet.get(afp1)); 326 327 328 double gp = misScore * m; //on average, penalty for a mismatch is misScore, no modification on score 329 if(g > 0) { 330 gp += gapExtend * g; 331 } 332 if(gp < maxPenalty) gp = maxPenalty; //penalty cut-off 333 //note: use < (smaller) instead of >, because maxPenalty is a negative number 334 335 double d; 336 d = calAfpDis(afp1, afp2,params, afpChain); 337 //note: the 'dis' value is numerically equivalent to the 'rms' with exceptions 338 339 boolean ch = false; 340 double tp = 0.0; 341 if(d >= disCut) { 342 tp = torsionPenalty; 343 ch = true; 344 } //use the variation of the distances between AFPs 345 else if(d > disCut - disSmooth) { 346 double wt = Math.sqrt((d - disCut + disSmooth) / disSmooth); 347 //using sqrt: penalty increase with dis more quicker than linear function 348 tp = torsionPenalty * wt; 349 } 350 351 dvar = d; 352 conn = tp + gp; 353 354 afpChain.setConn(conn); 355 afpChain.setDVar(dvar); 356 return ch; 357 } 358 359 /** 360 * return the gaps between this and afp 361 * requiring afp1 > afp2 362 * ( operator % in C) 363 */ 364 365 private static int calcGap(AFP afp1 , AFP afp2) 366 { 367 int g = (afp1.getP1() - afp2.getP1()) - (afp1.getP2() - afp2.getP2()); 368 if(g < 0) g = -g; 369 return g; 370 } 371 372 /** 373 * return the mis-matched between this and afp 374 * requiring AFP1 > afp2 375 * (operator / in C) 376 * @param afp1 377 * @param afp2 378 * @return 379 */ 380 381 //-------------------------------------------- 382 private static int calcMismatch(AFP afp1, AFP afp2) 383 { 384 int l1 = afp1.getP1() - afp2.getP1() - afp2.getFragLen(); 385 int l2 = afp1.getP2() - afp2.getP2() - afp2.getFragLen(); 386 return (l1 > l2?l2:l1); 387 } 388 389 /** 390 //return the root mean square of the distance matrix between the residues 391 //from the segments that form the given AFP list 392 //this value can be a measurement (2) for the connectivity of the AFPs 393 //and its calculation is quicker than the measurement (1), rmsd 394 //currently only deal with AFP pair 395// 396// |-d1--| 397// |--d2---| 398// |---d3----| 399 //----------------------------------------------------------------------- 400 //this module is optimized 401 * 402 * @param afp1 403 * @param afp2 404 * @return 405 */ 406 private static double calAfpDis(int afp1, int afp2, FatCatParameters params, AFPChain afpChain) 407 { 408 409 List<AFP> afpSet = afpChain.getAfpSet(); 410 411 Matrix disTable1 = afpChain.getDisTable1(); 412 Matrix disTable2 = afpChain.getDisTable2(); 413 414 int fragLen = params.getFragLen(); 415 double afpDisCut = params.getAfpDisCut(); 416 double disCut = params.getDisCut(); 417 double fragLenSq = params.getFragLenSq(); 418 419 int i, j, ai, bi, aj, bj; 420 double d; 421 double rms = 0; 422 for(i = 0; i < fragLen; i ++) { 423 ai = afpSet.get(afp1).getP1() + i; 424 bi = afpSet.get(afp1).getP2() + i; 425 for(j = 0; j < fragLen; j ++) { 426 aj = afpSet.get(afp2).getP1() + j; 427 bj = afpSet.get(afp2).getP2() + j; 428 d = disTable1.get(aj,ai) - disTable2.get(bj,bi); 429 rms += d * d; 430 if(rms > afpDisCut) { return (disCut); } 431 } 432 } 433 return (Math.sqrt(rms / fragLenSq)); 434 } 435 436 437 /** 438 * derive the optimal chaining of AFPs by trace-back 439 */ 440 private static void traceBack(int[] pre, int currafp0, int twist, FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[]ca2) 441 { 442 443 afpChain.setConn(0d); 444 afpChain.setDVar(0d); 445 446 int minLen = afpChain.getMinLen(); 447 List<AFP> afpSet = afpChain.getAfpSet(); 448 449 int afpChainLen = 0; 450 451 //trace-back from currafp (maxsco) 452 int[] afpchain = new int[minLen]; 453 int[] afptwibin = new int[minLen]; 454 double[] afptwilist = new double[minLen]; 455 int currafp = currafp0; 456 int s = 0; 457 458 afptwibin[s] = 0; 459 afpchain[s ++] = currafp; 460 461 boolean isConnected = false; 462 int prevafp; 463 // int alnlen = afpSet.get(afpchain[s]).getFragLen(); 464 465 //Double conn = new Double(0) ; 466 //Double dvar = new Double(0); 467 468 while((prevafp = pre[currafp]) != -1) { 469 470 isConnected = afpPairConn(prevafp, currafp, params,afpChain); 471 472 if ( isConnected ) 473 afptwibin[s - 1] = 1; 474 else 475 afptwibin[s - 1] = 0; 476 Double dvar = afpChain.getDVar(); 477 afptwilist[s - 1] = dvar; 478 //note s - 1: the transformation of prevafp-currafp is recorded in currafp 479 currafp = prevafp; 480 // alnlen += afpSet.get(currafp).getFragLen(); 481 afpchain[s ++] = currafp; 482 } 483 484 afpChainLen = s; 485 afpChain.setAfpChainLen(afpChainLen); 486 487 //first afp without transformation 488 if ( isConnected) 489 afptwibin[s - 1] = 1; 490 else 491 afptwibin[s - 1] = 0; 492 493 //if(debug) 494 // System.out.println(String.format("including %d AFPs, %d residues\n", afpChainLen, alnlen)); 495 496 //record the optimal alignment in afpChainList (afpChainLen) 497 498 int[] afpChainList = afpChain.getAfpChainList(); 499 double[] afpChainTwiBin = afpChain.getAfpChainTwiBin(); 500 double[] afpChainTwiList = afpChain.getAfpChainTwiList(); 501 502 if(afpChainList == null) { 503 afpChainList = new int[s]; 504 afpChain.setAfpChainList(afpChainList); 505 afpChainTwiBin = new double[s]; 506 afpChain.setAfpChainTwiBin(afpChainTwiBin); 507 afpChainTwiList = new double[s]; 508 afpChain.setAfpChainTwiList(afpChainTwiList); 509 } 510 511 int afpChainTwiNum = afpChain.getAfpChainTwiNum(); 512 513 int i; 514 for(i = 0; i < s; i ++) { 515 afpChainList[i] = afpchain[s - 1 - i]; 516 afpChainTwiBin[i] = afptwibin[s - 1 - i]; 517 afpChainTwiList[i] = afptwilist[s - 1 - i]; 518 afpChainTwiNum += afptwibin[s - 1 - i]; 519 } 520 521 if(afpChainTwiNum != twist) { 522 System.err.println(String.format("AFPChainer Warning: the twists number is not consistent %d %d\n", afpChainTwiNum, twist)); 523 } 524 525 double alignScore = afpChain.getAlignScore(); 526 527 double checkscore = afpSet.get(afpChainList[0]).getScore(); 528 for(i = 1; i < afpChainLen; i ++) { 529 isConnected = afpPairConn(afpChainList[i - 1], afpChainList[i], params,afpChain); 530 checkscore = checkscore + afpSet.get(afpChainList[i]).getScore() + afpChain.getConn(); 531 } 532 if(Math.abs(checkscore - alignScore) > 1e-4) { 533 System.err.println(String.format("AFPChainer Warning: fail in alignment score checking %.4f %.4f\n", alignScore, checkscore)); 534 } 535 536 if ( debug ) 537 System.out.println("traceBack:" + afpChainLen + " " + afpChainList.length); 538 double rmsd = calAfpRmsd(afpChainLen, afpChainList,0, afpChain,ca1,ca2); 539 afpChain.setChainRmsd(rmsd); 540 if (debug ) 541 System.out.println("Chain RMSD: " + rmsd); 542 int b1 = 0; 543 int bk = 0; 544 int a, b; 545 afpChain.setChainLen( 0); 546 int chainLen = afpChain.getChainLen(); 547 int[] block2Afp = afpChain.getBlock2Afp(); 548 549 550 double[] blockRmsd = afpChain.getBlockRmsd(); 551 int[] blockSize = afpChain.getBlockSize(); 552 553 block2Afp[0] = 0; 554 for(i = 0; i < afpChainLen; i ++) { 555 a = afpChainList[i]; 556 chainLen += afpSet.get(a).getFragLen(); 557 if(i > 0) { 558 b = afpChainList[i - 1]; 559 int misLen = afpChain.getMisLen(); 560 misLen += calcMismatch(afpSet.get(a),afpSet.get(b)); 561 afpChain.setMisLen(misLen); 562 int gapLen = afpChain.getGapLen(); 563 gapLen += calcGap(afpSet.get(a),afpSet.get(b)); 564 afpChain.setGapLen(gapLen); 565 566 } 567 568 if(afpChainTwiBin[i] == 1) { 569 if (debug) 570 System.out.println(" ** calAfpTmsd : afpChainWtiBin == 1 : i: "+i+" i-b1: " + (i-b1) + " b1: " + b1 + " afpChainList.len: " + afpChainList.length); 571 572 //int len = afpChainList.length - b1 +1; 573 //int[] fakeList = new int[len]; 574 //int pos = -1; 575 //for ( int fPos = b1 ; fPos< afpChainList.length ; fPos++){ 576 // pos++; 577 // fakeList[pos] = afpChainList[fPos]; 578 //} 579 if (debug ) 580 System.err.println("calculation calAfpRmsd " + i + " " + b1 + " " ); 581 582 583 rmsd = calAfpRmsd(i - b1, afpChainList,b1, afpChain, ca1, ca2); 584 585 586 blockRmsd[bk] = rmsd; 587 blockSize[bk] = i - b1; 588 b1 = i; 589 590 //System.out.println("block2Afp.length:"+ block2Afp.length + " " + bk + " " + i + " " + afpChain.getMaxTra() ); 591 block2Afp[++bk] = i; //next-block 592 } 593 } 594 afpChain.setBlock2Afp(block2Afp); 595 afpChain.setChainLen(chainLen); 596 597 if (debug) 598 System.out.println("after loop over all afpChainList " + (i-b1) + " " + b1); 599 600 rmsd = calAfpRmsd(i - b1, afpChainList, b1, afpChain,ca1,ca2); 601 if (debug) 602 System.out.println("*** i:" + i + " b1: " + b1 + " i-b1 " + (i-b1)); 603 604 605 606 // argh this is the block RMSD, not the chain RMSD! 607 //afpChain.setChainRmsd(rmsd); 608 //rmsd = calAfpRmsd(i - b1, afpChainList[b1]); 609 blockSize[bk] = i - b1; 610 blockRmsd[bk] = rmsd; 611 612 afpChain.setBlockSize(blockSize); 613 afpChain.setBlockRmsd(blockRmsd); 614 int blockNum = afpChain.getBlockNum(); 615 blockNum = ++bk; 616 if ( debug) 617 System.err.println("AFPChainser setBlockNUm:" + blockNum); 618 afpChain.setBlockNum(blockNum); 619 afpChain.setAfpChainList(afpChainList); 620 afpChain.setAfpChainTwiList(afpChainTwiList); 621 622 } 623 624 /** 625 //return the rmsd of the residues from the segments that form the given AFP list 626 //this value can be a measurement (1) for the connectivity of the AFPs 627 * 628 * @param afpn 629 * @param afpPositions the positions of AFPs to work on. 630 * @param listStart the starting position in the list of AFPs 631 * @param afpChain 632 * @param ca1 633 * @param ca2 634 * @return rmsd 635 */ 636 637 protected static double calAfpRmsd(int afpn, int[] afpPositions, int listStart, AFPChain afpChain,Atom[] ca1, Atom[] ca2) 638 { 639 640 641 if (debug) 642 System.out.println("XXX calling afp2res "+ afpn + " " + afpPositions.length); 643 644 int focusResn = AFPTwister.afp2Res(afpChain,afpn, afpPositions, listStart); 645 646 647 int[] focusRes1 = afpChain.getFocusRes1(); 648 int[] focusRes2 = afpChain.getFocusRes2(); 649 650 if (debug) 651 System.out.println("XXX calculating calAfpRmsd: " + focusResn + " " + focusRes1.length + " " + focusRes2.length + " " + afpChain.getMinLen() + " " ); 652 double rmsd = getRmsd(focusResn, focusRes1, focusRes2 , afpChain,ca1, ca2); 653 654 return rmsd; 655 } 656 657 658 659 /** the the RMSD for the residues defined in the two arrays 660 * 661 * @param focusResn 662 * @param focusRes1 663 * @param focusRes2 664 * @return 665 */ 666 private static double getRmsd(int focusResn, int[] focusRes1, int[] focusRes2, AFPChain afpChain, Atom[] ca1, Atom[] ca2){ 667 668 669 670 Atom[] tmp1 = new Atom[focusResn]; 671 Atom[] tmp2 = new Atom[focusResn]; 672 673 for ( int i =0 ; i< focusResn;i++){ 674 tmp1[i] = ca1[focusRes1[i]]; 675 tmp2[i] = (Atom)ca2[focusRes2[i]].clone(); 676 if (tmp1[i].getCoords() == null){ 677 System.err.println("tmp1 got null: " +i + " pos: " + focusRes1[i]); 678 } 679 if (tmp2[i].getCoords() == null){ 680 System.err.println("tmp1 got null: " +i + " pos: " + focusRes2[i]); 681 } 682 //XX 683 //tmp2[i].setParent((Group) ca2[focusRes2[i]].getParent().clone()); 684 } 685 double rmsd = 99; 686 try { 687 rmsd = getRmsd(tmp1,tmp2); 688 689 } catch (Exception e){ 690 e.printStackTrace(); 691 } 692 return rmsd; 693 694 } 695 /** Calculate the RMSD for two sets of atoms. Rotates the 2nd atom set so make sure this does not cause problems later 696 * 697 * 698 * @param catmp1 699 * @return 700 */ 701 private static double getRmsd(Atom[] catmp1, Atom[] catmp2) throws StructureException{ 702 703 Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(catmp1), 704 Calc.atomsToPoints(catmp2)); 705 706 Calc.transform(catmp2, trans); 707 708 // if ( showAlig) { 709 // StructureAlignmentJmol jmol = new StructureAlignmentJmol(); 710 // jmol.setTitle("AFPCHainer: getRmsd" + rmsd); 711 // 712 // Chain c1 = new ChainImpl(); 713 // c1.setName("A"); 714 // for ( Atom a : catmp1){ 715 // c1.addGroup(a.getParent()); 716 // } 717 // 718 // Chain c2 = new ChainImpl(); 719 // c2.setName("B"); 720 // for ( Atom a : catmp2){ 721 // c2.addGroup(a.getParent()); 722 // } 723 // 724 // Structure fake = new StructureImpl(); 725 // fake.setPDBCode("AFPCHainer: getRmsd" + rmsd); 726 // List<Chain> model1 = new ArrayList<Chain>(); 727 // model1.add(c1); 728 // List<Chain> model2 = new ArrayList<Chain>(); 729 // model2.add(c2); 730 // fake.addModel(model1); 731 // fake.addModel(model2); 732 // fake.setNmr(true); 733 // 734 // jmol.setStructure(fake); 735 // jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " + 736 // "select not protein and not solvent; spacefill on;"); 737 // jmol.evalString("select */1 ; color red; model 1; "); 738 // 739 // // now show both models again. 740 // jmol.evalString("model 0;"); 741 // } 742 743 return Calc.rmsd(catmp1,catmp2); 744 } 745 746}