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