org.biojava.nbio.structure.align.AFPTwister Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of biojava-structure Show documentation
Show all versions of biojava-structure Show documentation
The protein structure modules of BioJava.
/* This class is based on the original FATCAT implementation by
*
* Yuzhen Ye & Adam Godzik (2003)
* Flexible structure alignment by chaining aligned fragment pairs allowing twists.
* Bioinformatics vol.19 suppl. 2. ii246-ii255.
* http://www.ncbi.nlm.nih.gov/pubmed/14534198
*
*
* Thanks to Yuzhen Ye and A. Godzik for granting permission to freely use and redistribute this code.
*
* This code may be freely distributed and modified under the
* terms of the GNU Lesser General Public Licence. This should
* be distributed with the code. If you do not have a copy,
* see:
*
* http://www.gnu.org/copyleft/lesser.html
*
* Copyright for this code is held jointly by the individual
* authors. These should be listed in @author doc comments.
*
*
* Created on Jun 17, 2009
* Created by Andreas Prlic - RCSB PDB
*
*/
package org.biojava.nbio.structure.align;
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.model.AFP;
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.ArrayList;
import java.util.List;
//import org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol;
public class AFPTwister
{
private final static Logger logger = LoggerFactory.getLogger(AFPTwister.class);
// private static final boolean showAlignmentSteps = false;
/**
* calculate the total rmsd of the blocks
* output a merged pdb file for both proteins
protein 1, in chain A
protein 2 is twisted according to the twists detected, in chain B
* @return twisted Groups
*/
public static Group[] twistPDB(AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws StructureException {
//--------------------------------------------------------
if ( afpChain.isShortAlign())
return new Group[0];
List afpSet = afpChain.getAfpSet();
int blockNum = afpChain.getBlockNum();
int i, b2, e2;
//superimposing according to the initial AFP-chaining
Atom[] origCA = StructureTools.cloneAtomArray(ca2);
Atom[] iniTwistPdb = StructureTools.cloneAtomArray(ca2);
int[] blockResSize = afpChain.getBlockResSize();
int[][][] blockResList = afpChain.getBlockResList();
int[] afpChainList = afpChain.getAfpChainList();
int[] block2Afp = afpChain.getBlock2Afp();
int[] blockSize = afpChain.getBlockSize();
int[] focusAfpList = afpChain.getFocusAfpList();
if ( focusAfpList == null){
focusAfpList = new int[afpChain.getMinLen()];
afpChain.setFocusAfpList(focusAfpList);
}
int focusAfpn = 0;
e2 = 0;
b2 = 0;
logger.debug("blockNUm at twister: ", blockNum);
for(int bk = 0; bk < blockNum; bk ++) {
// THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
// copies the atoms over to iniTwistPdb later on in modifyCod
transformOrigPDB(blockResSize[bk], blockResList[bk][0], blockResList[bk][1], ca1, ca2, null,-1);
//transform pro2 according to comparison of pro1 and pro2 at give residues
if(bk > 0) {
b2 = e2;
}
logger.debug("b2 is {} before modifyCon", b2);
if(bk < (blockNum - 1)) {
//bend at the middle of two consecutive AFPs
int afpPos = afpChainList[block2Afp[bk] + blockSize[bk] - 1];
AFP a1 = afpSet.get(afpPos);
e2 = a1.getP2() ;
int afpPos2 = afpChainList[block2Afp[bk + 1]];
AFP a2 = afpSet.get(afpPos2);
e2 = (a2.getP2() - e2)/ 2 + e2;
logger.debug("e2 : {}", e2);
} else{
// last one is until the end...
e2 = ca2.length;
}
// this copies the coordinates over into iniTwistPdb
cloneAtomRange(iniTwistPdb, ca2, b2, e2);
//bound[bk] = e2;
for(i = 0; i < blockSize[bk]; i ++) {
focusAfpList[focusAfpn] = afpChainList[block2Afp[bk] + i];
focusAfpn++;
}
}
int focusResn = afp2Res(afpChain, focusAfpn, focusAfpList, 0);
afpChain.setTotalLenIni(focusResn);
logger.debug(String.format("calrmsdini for %d residues", focusResn));
double totalRmsdIni = calCaRmsd(ca1,iniTwistPdb, focusResn, afpChain.getFocusRes1() , afpChain.getFocusRes2() );
afpChain.setTotalRmsdIni(totalRmsdIni);
logger.debug("got iniRMSD: {}", totalRmsdIni);
if ( totalRmsdIni == 5.76611141613097) {
logger.debug("{}", afpChain.getAlnseq1());
logger.debug("{}", afpChain.getAlnsymb());
logger.debug("{}", afpChain.getAlnseq2());
}
afpChain.setFocusAfpList(focusAfpList);
afpChain.setBlock2Afp(block2Afp);
afpChain.setAfpChainList(afpChainList);
return twistOptimized(afpChain,ca1,origCA);
}
/** superimposing according to the optimized alignment
*
* @param afpChain
* @param ca1
* @param ca2
* @return Group array twisted.
* @throws StructureException
*/
public static Group[] twistOptimized(AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws StructureException {
Atom[] optTwistPdb = new Atom[ca2.length];
int gPos = -1;
for (Atom a : ca2){
gPos++;
optTwistPdb[gPos] = a;
}
int blockNum = afpChain.getBlockNum();
int b2 = 0;
int e2 = 0;
int focusResn = 0;
int[] focusRes1 = afpChain.getFocusRes1();
int[] focusRes2 = afpChain.getFocusRes2();
if(focusRes1 == null) {
focusRes1 = new int[afpChain.getCa1Length()];
afpChain.setFocusRes1(focusRes1);
}
if(focusRes2 == null) {
focusRes2 = new int[afpChain.getCa2Length()];
afpChain.setFocusRes2(focusRes2);
}
int[] optLen = afpChain.getOptLen();
int[][][] optAln = afpChain.getOptAln();
for(int bk = 0; bk < blockNum; bk ++) {
// THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
// copies the atoms over to iniTwistPdb later on in modifyCod
transformOrigPDB(optLen[bk], optAln[bk][0], optAln[bk][1], ca1, ca2,afpChain,bk);
//transform pro2 according to comparison of pro1 and pro2 at give residues
if(bk > 0) { b2 = e2; }
if(bk < blockNum - 1) { //bend at the middle of two consecutive blocks
e2 = optAln[bk][1][optLen[bk] - 1];
e2 = (optAln[bk + 1][1][0] - e2)/ 2 + e2;
}
else { e2 = ca2.length; }
cloneAtomRange(optTwistPdb, ca2, b2, e2);
for(int i = 0; i < optLen[bk]; i ++) {
focusRes1[focusResn] = optAln[bk][0][i];
focusRes2[focusResn] = optAln[bk][1][i];
focusResn ++;
}
}
int totalLenOpt = focusResn;
logger.debug("calrmsdopt for {} residues", focusResn);
double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1, focusRes2);
logger.debug("got opt RMSD: {}", totalRmsdOpt);
int optLength = afpChain.getOptLength();
if(totalLenOpt != optLength) {
logger.warn("Final alignment length is different {} {}", totalLenOpt, optLength);
}
logger.debug("final alignment length {}, rmsd {}", focusResn, totalRmsdOpt);
afpChain.setTotalLenOpt(totalLenOpt);
afpChain.setTotalRmsdOpt(totalRmsdOpt);
return StructureTools.cloneGroups(optTwistPdb);
}
/**
* transform the coordinates in the ca2 according to the superimposing of the given position pairs.
* No Cloning, transforms input atoms.
*/
// orig name: transPdb
private static void transformOrigPDB(int n, int[] res1, int[] res2, Atom[] ca1, Atom[] ca2, AFPChain afpChain, int blockNr)
throws StructureException
{
logger.debug("transforming original coordinates {} len1: {} res1: {} len2: {} res2: {}", n, ca1.length, res1.length, ca2.length, res2.length);
Atom[] cod1 = getAtoms(ca1, res1, n,false);
Atom[] cod2 = getAtoms(ca2, res2, n,false);
//double *cod1 = pro1->Cod4Res(n, res1);
//double *cod2 = pro2->Cod4Res(n, res2);
Matrix r;
Atom t;
SVDSuperimposer svd = new SVDSuperimposer(cod1, cod2);
r = svd.getRotation();
t = svd.getTranslation();
logger.debug("transPdb: transforming orig coordinates with matrix: {}", r);
if ( afpChain != null ){
Matrix[] ms = afpChain.getBlockRotationMatrix();
if ( ms == null)
ms = new Matrix[afpChain.getBlockNum()];
ms[blockNr]= r;
Atom[] shifts = afpChain.getBlockShiftVector();
if ( shifts == null)
shifts = new Atom[afpChain.getBlockNum()];
shifts[blockNr] = t;
afpChain.setBlockRotationMatrix(ms);
afpChain.setBlockShiftVector(shifts);
}
for (Atom a : ca2){
Calc.rotate(a.getGroup(),r);
Calc.shift(a.getGroup(),t);
}
}
// like Cod4Res
// most likely the clone flag is not needed
private static Atom[] getAtoms(Atom[] ca, int[] positions, int length, boolean clone){
List atoms = new ArrayList();
for ( int i = 0 ; i < length ; i++){
int p = positions[i];
Atom a;
if ( clone ){
a = (Atom)ca[p].clone();
a.setGroup((Group)ca[p].getGroup().clone());
}
else {
a = ca[p];
}
atoms.add(a);
}
return atoms.toArray(new Atom[atoms.size()]);
}
/** Clones and copies the Atoms from p2 into p1 range is between r1 and r2
*
* @param p1
* @param p2
* @param r1
* @param r2
*/
// orig name: modifyCod
private static void cloneAtomRange(Atom[] p1, Atom[] p2, int r1, int r2)
throws StructureException
{
logger.debug("modifyCod from: {} to: {}", r1, r2);
// special clone method, can;t use StructureTools.cloneCAArray, since we access the data
// slightly differently here.
List model = new ArrayList();
for(int i = r1; i < r2; i ++) {
Group g = p2[i].getGroup();
Group newG = (Group)g.clone();
p1[i] = newG.getAtom(p2[i].getName());
Chain parentC = g.getChain();
Chain newChain = null;
for( Chain c: model){
if (c.getChainID().equals(parentC.getChainID())){
newChain = c;
break;
}
}
if ( newChain == null){
newChain = new ChainImpl();
newChain.setChainID(parentC.getChainID());
model.add(newChain);
}
newChain.addGroup(newG);
} //modify caCod
}
/**
* Return the rmsd of the CAs between the input pro and this protein, at given positions.
* quite similar to transPdb but while that one transforms the whole ca2, this one only works on the res1 and res2 positions.
*
* Modifies the coordinates in the second set of Atoms (pro).
*
* @return rmsd of CAs
*/
private static double calCaRmsd(Atom[] ca1, Atom[] pro, int resn, int[] res1, int[] res2) throws StructureException
{
Atom[] cod1 = getAtoms(ca1, res1,resn,false);
Atom[] cod2 = getAtoms(pro, res2,resn,false);
Matrix r;
Atom t;
if ( cod1.length == 0 || cod2.length == 0){
logger.info("length of atoms == 0!");
return 99;
}
SVDSuperimposer svd = new SVDSuperimposer(cod1, cod2);
r = svd.getRotation();
t = svd.getTranslation();
for (Atom a : cod2){
Calc.rotate(a.getGroup(), r);
Calc.shift(a.getGroup(), t);
}
return SVDSuperimposer.getRMS(cod1, cod2);
}
/**
* Set the list of equivalent residues in the two proteins given a list of AFPs
*
* WARNING: changes the values for FocusRes1, focusRes2 and FocusResn in afpChain!
*
* @param afpChain the AFPChain to store resuts
* @param afpn nr of afp
* @param afpPositions
* @param listStart
* @return nr of eq residues
*/
public static int afp2Res(AFPChain afpChain, int afpn, int[] afpPositions , int listStart )
{
int[] res1 = afpChain.getFocusRes1();
int[] res2 = afpChain.getFocusRes2();
int minLen = afpChain.getMinLen();
int n = 0;
List afpSet =afpChain.getAfpSet();
for(int i = listStart; i < listStart+afpn; i ++) {
int a = afpPositions[i];
for(int j = 0; j < afpSet.get(a).getFragLen(); j ++) {
if(n >= minLen) {
throw new RuntimeException("Error: too many residues in AFPChainer.afp2Res!");
}
res1[n] = afpSet.get(a).getP1() + j;
res2[n] = afpSet.get(a).getP2() + j;
n ++;
}
}
afpChain.setFocusRes1(res1);
afpChain.setFocusRes2(res2);
afpChain.setFocusResn(n);
if ( n == 0 ){
logger.warn("n=0!!! + " + afpn + " " + listStart + " " + afpPositions.length);
}
return n;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy