Skip to content

Commit cc79186

Browse files
lafitalafita
committed
Merge branch 'svdsuperimpIssue' into superpose
2 parents 85b4b69 + a8828ca commit cc79186

5 files changed

Lines changed: 120 additions & 41 deletions

File tree

‎biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/geometry/TestProteinSuperposition.java‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
package org.biojava.nbio.structure.test.geometry;
2+
3+
import static org.junit.Assert.*;
4+
5+
import java.io.IOException;
6+
7+
import javax.vecmath.Point3d;
8+
9+
import org.biojava.nbio.structure.Calc;
10+
import org.biojava.nbio.structure.Chain;
11+
import org.biojava.nbio.structure.Structure;
12+
import org.biojava.nbio.structure.StructureException;
13+
import org.biojava.nbio.structure.StructureIO;
14+
import org.biojava.nbio.structure.StructureTools;
15+
import org.biojava.nbio.structure.geometry.SuperPosition;
16+
import org.biojava.nbio.structure.geometry.SuperPositionQCP;
17+
import org.biojava.nbio.structure.geometry.SuperPositionQuat;
18+
import org.biojava.nbio.structure.geometry.SuperPositionSVD;
19+
import org.junit.BeforeClass;
20+
import org.junit.Test;
21+
22+
public class TestProteinSuperposition {
23+
24+
private static Point3d[] chain1;
25+
private static Point3d[] chain2;
26+
27+
@BeforeClass
28+
public static void setUpBeforeClass() throws StructureException, IOException {
29+
Structure s = StructureIO.getStructure("1smt");
30+
Chain origChainA = s.getPolyChainByPDB("A");
31+
Chain clonedChainA = (Chain) origChainA.clone();
32+
33+
chain1 = Calc.atomsToPoints(StructureTools.getAtomCAArray(origChainA));
34+
chain2 = Calc.atomsToPoints(StructureTools.getAtomCAArray(clonedChainA));
35+
36+
}
37+
38+
@Test
39+
public void testSuperpositionSVD() {
40+
41+
SuperPosition sup = new SuperPositionSVD(false);
42+
43+
double rmsd = sup.getRmsd(chain1, chain2);
44+
45+
assertEquals(0.0, rmsd, 0.0001);
46+
}
47+
48+
@Test
49+
public void testSuperpositionQCP() {
50+
51+
SuperPosition sup = new SuperPositionQCP(false);
52+
53+
double rmsd = sup.getRmsd(chain1, chain2);
54+
55+
assertEquals(0.0, rmsd, 0.0001);
56+
}
57+
58+
@Test
59+
public void testSuperpositionQuat() {
60+
61+
SuperPosition sup = new SuperPositionQuat(false);
62+
63+
double rmsd = sup.getRmsd(chain1, chain2);
64+
65+
assertEquals(0.0, rmsd, 0.0001);
66+
}
67+
}

‎biojava-structure/src/main/java/org/biojava/nbio/structure/Calc.java‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -998,7 +998,7 @@ public static final Atom createVirtualCBAtom(AminoAcid amino)
998998

999999
SuperPositionSVD svd = new SuperPositionSVD(false);
10001000

1001-
Matrix4d transform = svd.superpose(arr1, arr2);
1001+
Matrix4d transform = svd.superpose(Calc.atomsToPoints(arr1), Calc.atomsToPoints(arr2));
10021002
Matrix rotMatrix = Matrices.getRotationJAMA(transform);
10031003
Atom tranMatrix = getTranslationVector(transform);
10041004

‎biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java‎

Lines changed: 48 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,6 @@
2828
import org.biojava.nbio.structure.align.model.AFPChain;
2929
import org.biojava.nbio.structure.align.xml.AFPChainXMLParser;
3030
import org.biojava.nbio.structure.geometry.Matrices;
31-
import org.biojava.nbio.structure.geometry.SuperPositionSVD;
3231
import org.biojava.nbio.structure.geometry.SuperPositions;
3332
import org.biojava.nbio.structure.jama.Matrix;
3433
import org.slf4j.Logger;
@@ -855,7 +854,8 @@ public static AFPChain replaceOptAln(AFPChain afpChain, Atom[] ca1, Atom[] ca2,
855854
* contains much of the same code, but stores results in a CECalculator
856855
* instance rather than an AFPChain
857856
*/
858-
public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException {
857+
public static void updateSuperposition(AFPChain afpChain, Atom[] ca1,
858+
Atom[] ca2) throws StructureException {
859859

860860
//Update ca information, because the atom array might also be changed
861861
afpChain.setCa1Length(ca1.length);
@@ -878,30 +878,8 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
878878
// create new arrays for the subset of atoms in the alignment.
879879
Atom[] ca1aligned = new Atom[afpChain.getOptLength()];
880880
Atom[] ca2aligned = new Atom[afpChain.getOptLength()];
881-
int pos=0;
882-
int[] blockLens = afpChain.getOptLen();
883-
int[][][] optAln = afpChain.getOptAln();
884-
assert(afpChain.getBlockNum() <= optAln.length);
885-
886-
for (int block=0; block < afpChain.getBlockNum(); block++) {
887-
for(int i=0;i<blockLens[block];i++) {
888-
int pos1 = optAln[block][0][i];
889-
int pos2 = optAln[block][1][i];
890-
Atom a1 = ca1[pos1];
891-
Atom a2 = (Atom) ca2[pos2].clone();
892-
ca1aligned[pos] = a1;
893-
ca2aligned[pos] = a2;
894-
pos++;
895-
}
896-
}
897-
898-
// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
899-
if (pos != afpChain.getOptLength()){
900-
logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLength());
901-
// we need to resize the array, because we allocated too many atoms earlier on.
902-
ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
903-
ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
904-
}
881+
882+
fillAlignedAtomArrays(afpChain, ca1, ca2, ca1aligned, ca2aligned);
905883

906884
//Superimpose the two structures in correspondance to the new alignment
907885
Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
@@ -910,7 +888,6 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
910888
Matrix matrix = Matrices.getRotationJAMA(trans);
911889
Atom shift = Calc.getTranslationVector(trans);
912890

913-
914891
Matrix[] blockMxs = new Matrix[afpChain.getBlockNum()];
915892
Arrays.fill(blockMxs, matrix);
916893
afpChain.setBlockRotationMatrix(blockMxs);
@@ -928,6 +905,9 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
928905
double tmScore = Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
929906
afpChain.setTotalRmsdOpt(rmsd);
930907
afpChain.setTMScore(tmScore);
908+
909+
int[] blockLens = afpChain.getOptLen();
910+
int[][][] optAln = afpChain.getOptAln();
931911

932912
//Calculate the RMSD and TM score for every block of the new alignment
933913
double[] blockRMSD = new double[afpChain.getBlockNum()];
@@ -947,7 +927,7 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
947927
position++;
948928
}
949929
if (position != afpChain.getOptLen()[k]){
950-
logger.warn("AFPChainScorer getTMScore: Problems reconstructing block alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLen()[k]);
930+
logger.warn("AFPChainScorer getTMScore: Problems reconstructing block alignment! nr of loaded atoms is " + position + " but should be " + afpChain.getOptLen()[k]);
951931
// we need to resize the array, because we allocated too many atoms earlier on.
952932
ca1block = (Atom[]) resizeArray(ca1block, position);
953933
ca2block = (Atom[]) resizeArray(ca2block, position);
@@ -956,8 +936,8 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
956936
Matrix4d transb = SuperPositions.superpose(Calc.atomsToPoints(ca1block),
957937
Calc.atomsToPoints(ca2block));
958938

959-
Matrix matrixb = Matrices.getRotationJAMA(trans);
960-
Atom shiftb = Calc.getTranslationVector(trans);
939+
blockMxs[k] = Matrices.getRotationJAMA(trans);
940+
blockShifts[k] = Calc.getTranslationVector(trans);
961941

962942
Calc.transform(ca2block, transb);
963943

@@ -1380,4 +1360,42 @@ public static void shiftCA2(AFPChain afpChain, Atom[] ca2, Matrix m, Atom shift
13801360
twistedGroups[i]=g;
13811361
}
13821362
}
1363+
1364+
/**
1365+
* Fill the aligned Atom arrays with the equivalent residues in the afpChain.
1366+
* @param afpChain
1367+
* @param ca1
1368+
* @param ca2
1369+
* @param ca1aligned
1370+
* @param ca2aligned
1371+
*/
1372+
public static void fillAlignedAtomArrays(AFPChain afpChain, Atom[] ca1,
1373+
Atom[] ca2, Atom[] ca1aligned, Atom[] ca2aligned) {
1374+
1375+
int pos=0;
1376+
int[] blockLens = afpChain.getOptLen();
1377+
int[][][] optAln = afpChain.getOptAln();
1378+
assert(afpChain.getBlockNum() <= optAln.length);
1379+
1380+
for (int block=0; block < afpChain.getBlockNum(); block++) {
1381+
for(int i=0;i<blockLens[block];i++) {
1382+
int pos1 = optAln[block][0][i];
1383+
int pos2 = optAln[block][1][i];
1384+
Atom a1 = ca1[pos1];
1385+
Atom a2 = (Atom) ca2[pos2].clone();
1386+
ca1aligned[pos] = a1;
1387+
ca2aligned[pos] = a2;
1388+
pos++;
1389+
}
1390+
}
1391+
1392+
// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
1393+
if (pos != afpChain.getOptLength()){
1394+
logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLength());
1395+
// we need to resize the array, because we allocated too many atoms earlier on.
1396+
ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
1397+
ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
1398+
}
1399+
1400+
}
13831401
}

‎biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionSVD.java‎

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,6 @@
2525
import javax.vecmath.Matrix4d;
2626
import javax.vecmath.Point3d;
2727

28-
import org.biojava.nbio.structure.Atom;
29-
import org.biojava.nbio.structure.Calc;
3028
import org.biojava.nbio.structure.jama.Matrix;
3129
import org.biojava.nbio.structure.jama.SingularValueDecomposition;
3230

@@ -55,10 +53,6 @@ public SuperPositionSVD(boolean centered) {
5553
super(centered);
5654
}
5755

58-
public Matrix4d superpose(Atom[] fixed, Atom[] moved) {
59-
return superpose(Calc.atomsToPoints(fixed), Calc.atomsToPoints(moved));
60-
}
61-
6256
@Override
6357
public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
6458

@@ -70,7 +64,7 @@ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
7064
double[][] centAcoords = new double[][] { { cena.x, cena.y, cena.z } };
7165
Matrix centroidA = new Matrix(centAcoords);
7266

73-
double[][] centBcoords = new double[][] { { cenb.x, cenb.x, cenb.x } };
67+
double[][] centBcoords = new double[][] { { cenb.x, cenb.y, cenb.z } };
7468
Matrix centroidB = new Matrix(centBcoords);
7569

7670
// center at centroid

‎biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPosition.java‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -121,9 +121,9 @@ public void testSuperposition() {
121121
cloud1.get(c).length, svdTime, quatTime, qcpTime));
122122

123123
// Check that the transformation matrix was recovered
124-
assertTrue(transform.epsilonEquals(svdTransform, 0.01));
125-
assertTrue(transform.epsilonEquals(quatTransform, 0.01));
126-
assertTrue(transform.epsilonEquals(qcpTransform, 0.01));
124+
assertTrue(transform.epsilonEquals(svdTransform, 0.05));
125+
assertTrue(transform.epsilonEquals(quatTransform, 0.05));
126+
assertTrue(transform.epsilonEquals(qcpTransform, 0.05));
127127
}
128128

129129
}

0 commit comments

Comments
 (0)