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