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; 028 029import org.biojava.nbio.structure.*; 030import org.biojava.nbio.structure.align.model.AFP; 031import org.biojava.nbio.structure.align.model.AFPChain; 032import org.biojava.nbio.structure.geometry.Matrices; 033import org.biojava.nbio.structure.geometry.SuperPositions; 034import org.biojava.nbio.structure.jama.Matrix; 035import org.slf4j.Logger; 036import org.slf4j.LoggerFactory; 037 038import java.util.ArrayList; 039import java.util.List; 040 041import javax.vecmath.Matrix4d; 042 043//import org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol; 044 045public class AFPTwister { 046 private final static Logger logger = LoggerFactory 047 .getLogger(AFPTwister.class); 048 049 // private static final boolean showAlignmentSteps = false; 050 051 /** 052 * calculate the total rmsd of the blocks output a merged pdb file for both 053 * proteins protein 1, in chain A protein 2 is twisted according to the 054 * twists detected, in chain B 055 * 056 * @return twisted Groups 057 */ 058 public static Group[] twistPDB(AFPChain afpChain, Atom[] ca1, Atom[] ca2) 059 throws StructureException { 060 // -------------------------------------------------------- 061 062 if (afpChain.isShortAlign()) 063 return new Group[0]; 064 065 List<AFP> afpSet = afpChain.getAfpSet(); 066 067 int blockNum = afpChain.getBlockNum(); 068 069 int i, b2, e2; 070 071 // superimposing according to the initial AFP-chaining 072 Atom[] origCA = StructureTools.cloneAtomArray(ca2); 073 Atom[] iniTwistPdb = StructureTools.cloneAtomArray(ca2); 074 075 int[] blockResSize = afpChain.getBlockResSize(); 076 077 int[][][] blockResList = afpChain.getBlockResList(); 078 079 int[] afpChainList = afpChain.getAfpChainList(); 080 int[] block2Afp = afpChain.getBlock2Afp(); 081 int[] blockSize = afpChain.getBlockSize(); 082 int[] focusAfpList = afpChain.getFocusAfpList(); 083 084 if (focusAfpList == null) { 085 focusAfpList = new int[afpChain.getMinLen()]; 086 afpChain.setFocusAfpList(focusAfpList); 087 } 088 089 int focusAfpn = 0; 090 e2 = 0; 091 b2 = 0; 092 093 logger.debug("blockNUm at twister: ", blockNum); 094 095 for (int bk = 0; bk < blockNum; bk++) { 096 097 // THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING... 098 // copies the atoms over to iniTwistPdb later on in modifyCod 099 100 transformOrigPDB(blockResSize[bk], blockResList[bk][0], 101 blockResList[bk][1], ca1, ca2, null, -1); 102 103 // transform pro2 according to comparison of pro1 and pro2 at give 104 // residues 105 if (bk > 0) { 106 b2 = e2; 107 } 108 109 logger.debug("b2 is {} before modifyCon", b2); 110 111 if (bk < (blockNum - 1)) { 112 113 // bend at the middle of two consecutive AFPs 114 int afpPos = afpChainList[block2Afp[bk] + blockSize[bk] - 1]; 115 AFP a1 = afpSet.get(afpPos); 116 e2 = a1.getP2(); 117 118 int afpPos2 = afpChainList[block2Afp[bk + 1]]; 119 AFP a2 = afpSet.get(afpPos2); 120 e2 = (a2.getP2() - e2) / 2 + e2; 121 logger.debug("e2 : {}", e2); 122 } else { 123 // last one is until the end... 124 e2 = ca2.length; 125 } 126 127 // this copies the coordinates over into iniTwistPdb 128 cloneAtomRange(iniTwistPdb, ca2, b2, e2); 129 130 // bound[bk] = e2; 131 for (i = 0; i < blockSize[bk]; i++) { 132 focusAfpList[focusAfpn] = afpChainList[block2Afp[bk] + i]; 133 focusAfpn++; 134 } 135 } 136 137 int focusResn = afp2Res(afpChain, focusAfpn, focusAfpList, 0); 138 139 afpChain.setTotalLenIni(focusResn); 140 141 logger.debug(String.format("calrmsdini for %d residues", focusResn)); 142 143 double totalRmsdIni = calCaRmsd(ca1, iniTwistPdb, focusResn, 144 afpChain.getFocusRes1(), afpChain.getFocusRes2()); 145 146 afpChain.setTotalRmsdIni(totalRmsdIni); 147 logger.debug("got iniRMSD: {}", totalRmsdIni); 148 if (totalRmsdIni == 5.76611141613097) { 149 logger.debug("{}", afpChain.getAlnseq1()); 150 logger.debug("{}", afpChain.getAlnsymb()); 151 logger.debug("{}", afpChain.getAlnseq2()); 152 } 153 154 afpChain.setFocusAfpList(focusAfpList); 155 afpChain.setBlock2Afp(block2Afp); 156 afpChain.setAfpChainList(afpChainList); 157 158 return twistOptimized(afpChain, ca1, origCA); 159 } 160 161 /** 162 * superimposing according to the optimized alignment 163 * 164 * @param afpChain 165 * @param ca1 166 * @param ca2 167 * @return Group array twisted. 168 * @throws StructureException 169 */ 170 public static Group[] twistOptimized(AFPChain afpChain, Atom[] ca1, 171 Atom[] ca2) throws StructureException { 172 173 Atom[] optTwistPdb = new Atom[ca2.length]; 174 175 int gPos = -1; 176 for (Atom a : ca2) { 177 gPos++; 178 optTwistPdb[gPos] = a; 179 } 180 181 int blockNum = afpChain.getBlockNum(); 182 183 int b2 = 0; 184 int e2 = 0; 185 int focusResn = 0; 186 int[] focusRes1 = afpChain.getFocusRes1(); 187 int[] focusRes2 = afpChain.getFocusRes2(); 188 189 if (focusRes1 == null) { 190 focusRes1 = new int[afpChain.getCa1Length()]; 191 afpChain.setFocusRes1(focusRes1); 192 } 193 if (focusRes2 == null) { 194 focusRes2 = new int[afpChain.getCa2Length()]; 195 afpChain.setFocusRes2(focusRes2); 196 } 197 198 int[] optLen = afpChain.getOptLen(); 199 int[][][] optAln = afpChain.getOptAln(); 200 201 for (int bk = 0; bk < blockNum; bk++) { 202 // THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING... 203 // copies the atoms over to iniTwistPdb later on in modifyCod 204 transformOrigPDB(optLen[bk], optAln[bk][0], optAln[bk][1], ca1, 205 ca2, afpChain, bk); 206 207 // transform pro2 according to comparison of pro1 and pro2 at give 208 // residues 209 if (bk > 0) { 210 b2 = e2; 211 } 212 if (bk < blockNum - 1) { // bend at the middle of two consecutive 213 // blocks 214 e2 = optAln[bk][1][optLen[bk] - 1]; 215 e2 = (optAln[bk + 1][1][0] - e2) / 2 + e2; 216 } else { 217 e2 = ca2.length; 218 } 219 cloneAtomRange(optTwistPdb, ca2, b2, e2); 220 for (int i = 0; i < optLen[bk]; i++) { 221 focusRes1[focusResn] = optAln[bk][0][i]; 222 focusRes2[focusResn] = optAln[bk][1][i]; 223 focusResn++; 224 } 225 } 226 int totalLenOpt = focusResn; 227 logger.debug("calrmsdopt for {} residues", focusResn); 228 double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1, 229 focusRes2); 230 logger.debug("got opt RMSD: {}", totalRmsdOpt); 231 int optLength = afpChain.getOptLength(); 232 233 if (totalLenOpt != optLength) { 234 logger.warn("Final alignment length is different {} {}", 235 totalLenOpt, optLength); 236 } 237 logger.debug("final alignment length {}, rmsd {}", focusResn, 238 totalRmsdOpt); 239 240 afpChain.setTotalLenOpt(totalLenOpt); 241 afpChain.setTotalRmsdOpt(totalRmsdOpt); 242 243 return StructureTools.cloneGroups(optTwistPdb); 244 245 } 246 247 /** 248 * transform the coordinates in the ca2 according to the superimposing of 249 * the given position pairs. No Cloning, transforms input atoms. 250 */ 251 // orig name: transPdb 252 private static void transformOrigPDB(int n, int[] res1, int[] res2, 253 Atom[] ca1, Atom[] ca2, AFPChain afpChain, int blockNr) 254 throws StructureException { 255 logger.debug( 256 "transforming original coordinates {} len1: {} res1: {} len2: {} res2: {}", 257 n, ca1.length, res1.length, ca2.length, res2.length); 258 259 Atom[] cod1 = getAtoms(ca1, res1, n, false); 260 Atom[] cod2 = getAtoms(ca2, res2, n, false); 261 262 // double *cod1 = pro1->Cod4Res(n, res1); 263 // double *cod2 = pro2->Cod4Res(n, res2); 264 265 Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(cod1), 266 Calc.atomsToPoints(cod2)); 267 268 Matrix r = Matrices.getRotationJAMA(transform); 269 Atom t = Calc.getTranslationVector(transform); 270 271 logger.debug("transPdb: transforming orig coordinates with matrix: {}", 272 r); 273 274 if (afpChain != null) { 275 Matrix[] ms = afpChain.getBlockRotationMatrix(); 276 if (ms == null) 277 ms = new Matrix[afpChain.getBlockNum()]; 278 279 ms[blockNr] = r; 280 281 Atom[] shifts = afpChain.getBlockShiftVector(); 282 if (shifts == null) 283 shifts = new Atom[afpChain.getBlockNum()]; 284 shifts[blockNr] = t; 285 286 afpChain.setBlockRotationMatrix(ms); 287 afpChain.setBlockShiftVector(shifts); 288 } 289 290 for (Atom a : ca2) 291 Calc.transform(a.getGroup(), transform); 292 293 } 294 295 // like Cod4Res 296 // most likely the clone flag is not needed 297 private static Atom[] getAtoms(Atom[] ca, int[] positions, int length, 298 boolean clone) { 299 300 List<Atom> atoms = new ArrayList<Atom>(); 301 for (int i = 0; i < length; i++) { 302 int p = positions[i]; 303 Atom a; 304 if (clone) { 305 a = (Atom) ca[p].clone(); 306 a.setGroup((Group) ca[p].getGroup().clone()); 307 } else { 308 a = ca[p]; 309 } 310 atoms.add(a); 311 } 312 return atoms.toArray(new Atom[atoms.size()]); 313 } 314 315 /** 316 * Clones and copies the Atoms from p2 into p1 range is between r1 and r2 317 * 318 * @param p1 319 * @param p2 320 * @param r1 321 * @param r2 322 */ 323 // orig name: modifyCod 324 private static void cloneAtomRange(Atom[] p1, Atom[] p2, int r1, int r2) 325 throws StructureException { 326 327 logger.debug("modifyCod from: {} to: {}", r1, r2); 328 329 // special clone method, can;t use StructureTools.cloneCAArray, since we 330 // access the data 331 // slightly differently here. 332 List<Chain> model = new ArrayList<Chain>(); 333 for (int i = r1; i < r2; i++) { 334 335 Group g = p2[i].getGroup(); 336 Group newG = (Group) g.clone(); 337 338 p1[i] = newG.getAtom(p2[i].getName()); 339 Chain parentC = g.getChain(); 340 341 Chain newChain = null; 342 343 for (Chain c : model) { 344 if (c.getName().equals(parentC.getName())) { 345 newChain = c; 346 break; 347 } 348 } 349 350 if (newChain == null) { 351 newChain = new ChainImpl(); 352 newChain.setName(parentC.getName()); 353 model.add(newChain); 354 } 355 356 newChain.addGroup(newG); 357 358 } // modify caCod 359 360 } 361 362 /** 363 * Return the rmsd of the CAs between the input pro and this protein, at 364 * given positions. quite similar to transPdb but while that one transforms 365 * the whole ca2, this one only works on the res1 and res2 positions. 366 * 367 * Modifies the coordinates in the second set of Atoms (pro). 368 * 369 * @return rmsd of CAs 370 */ 371 private static double calCaRmsd(Atom[] ca1, Atom[] pro, int resn, 372 int[] res1, int[] res2) throws StructureException { 373 374 Atom[] cod1 = getAtoms(ca1, res1, resn, false); 375 Atom[] cod2 = getAtoms(pro, res2, resn, false); 376 377 if (cod1.length == 0 || cod2.length == 0) { 378 logger.info("length of atoms == 0!"); 379 return 99; 380 } 381 382 Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(cod1), 383 Calc.atomsToPoints(cod2)); 384 385 for (Atom a : cod2) 386 Calc.transform(a.getGroup(), transform); 387 388 return Calc.rmsd(cod1, cod2); 389 } 390 391 /** 392 * Set the list of equivalent residues in the two proteins given a list of 393 * AFPs 394 * 395 * WARNING: changes the values for FocusRes1, focusRes2 and FocusResn in 396 * afpChain! 397 * 398 * @param afpChain 399 * the AFPChain to store resuts 400 * @param afpn 401 * nr of afp 402 * @param afpPositions 403 * @param listStart 404 * @return nr of eq residues 405 */ 406 407 public static int afp2Res(AFPChain afpChain, int afpn, int[] afpPositions, 408 int listStart) { 409 int[] res1 = afpChain.getFocusRes1(); 410 int[] res2 = afpChain.getFocusRes2(); 411 int minLen = afpChain.getMinLen(); 412 413 int n = 0; 414 415 List<AFP> afpSet = afpChain.getAfpSet(); 416 417 for (int i = listStart; i < listStart + afpn; i++) { 418 int a = afpPositions[i]; 419 for (int j = 0; j < afpSet.get(a).getFragLen(); j++) { 420 if (n >= minLen) { 421 throw new RuntimeException( 422 "Error: too many residues in AFPChainer.afp2Res!"); 423 } 424 res1[n] = afpSet.get(a).getP1() + j; 425 res2[n] = afpSet.get(a).getP2() + j; 426 n++; 427 } 428 } 429 430 afpChain.setFocusRes1(res1); 431 afpChain.setFocusRes2(res2); 432 afpChain.setFocusResn(n); 433 434 if (n == 0) { 435 logger.warn("n=0!!! + " + afpn + " " + listStart + " " 436 + afpPositions.length); 437 } 438 return n; 439 } 440 441}