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