diff --git a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/StructureTest.java b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/StructureTest.java
index ec2f589929..76edc8dcd4 100644
--- a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/StructureTest.java
+++ b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/StructureTest.java
@@ -21,16 +21,17 @@
package org.biojava.nbio.structure.test;
import org.biojava.nbio.structure.*;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.io.FileParsingParameters;
import org.biojava.nbio.structure.io.PDBFileParser;
import org.biojava.nbio.structure.io.SSBondImpl;
-import org.biojava.nbio.structure.jama.Matrix;
-
import java.io.IOException;
import java.io.InputStream;
import java.util.List;
import java.util.Set;
+import javax.vecmath.Matrix4d;
+
import org.junit.BeforeClass;
import org.junit.Test;
@@ -41,8 +42,6 @@
* @author Andreas Prlic
* @since 1.5
*/
-
-
public class StructureTest {
private static Structure structure;
@@ -266,18 +265,12 @@ public void testMutation() throws Exception {
atoms2[1] = g2.getAtom("CA");
atoms2[2] = g2.getAtom("CB");
+ Matrix4d transform = SuperPositions.superpose(
+ Calc.atomsToPoints(atoms1), Calc.atomsToPoints(atoms2));
- SVDSuperimposer svds = new SVDSuperimposer(atoms1,atoms2);
-
-
- Matrix rotMatrix = svds.getRotation();
- Atom tran = svds.getTranslation();
-
- Group newGroup = (Group)g2.clone();
-
- Calc.rotate(newGroup,rotMatrix);
+ Group newGroup = (Group) g2.clone();
- Calc.shift(newGroup,tran);
+ Calc.transform(newGroup, transform);
Atom ca1 = g1.getAtom("CA");
Atom oldca2 = g2.getAtom("CA");
diff --git a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/FlipAFPChainTest.java b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/FlipAFPChainTest.java
index 9f62b85778..f08d8c3a5a 100644
--- a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/FlipAFPChainTest.java
+++ b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/FlipAFPChainTest.java
@@ -25,6 +25,7 @@
package org.biojava.nbio.structure.test.align;
import junit.framework.TestCase;
+
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.StructureAlignment;
import org.biojava.nbio.structure.align.StructureAlignmentFactory;
@@ -150,7 +151,7 @@ private double getRMSD(AFPChain afpChain, Atom[] ca1, Atom[] ca2)
assertTrue(catmp1.length == afpChain.getNrEQR());
- return SVDSuperimposer.getRMS(catmp1,catmp2);
+ return Calc.rmsd(catmp1,catmp2);
}
public static void rotateAtoms2(AFPChain afpChain,Atom[] ca2){
diff --git a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/FlipAFPChainTest.java b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/FlipAFPChainTest.java
index 8aa21461c4..a7b0c15777 100644
--- a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/FlipAFPChainTest.java
+++ b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/FlipAFPChainTest.java
@@ -25,6 +25,7 @@
package org.biojava.nbio.structure.test.align.fatcat;
import junit.framework.TestCase;
+
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.StructureAlignment;
import org.biojava.nbio.structure.align.StructureAlignmentFactory;
@@ -171,7 +172,7 @@ private double getRMSD(AFPChain afpChain, Atom[] ca1, Atom[] ca2)
assertTrue(catmp1.length == afpChain.getNrEQR());
- return SVDSuperimposer.getRMS(catmp1,catmp2);
+ return Calc.rmsd(catmp1,catmp2);
}
public static void rotateAtoms2(AFPChain afpChain,Atom[] ca2){
diff --git a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/TestFlexibleRotationMatrices.java b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/TestFlexibleRotationMatrices.java
index 551a6d5288..b1ff3ddb8b 100644
--- a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/TestFlexibleRotationMatrices.java
+++ b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/align/fatcat/TestFlexibleRotationMatrices.java
@@ -25,6 +25,7 @@
package org.biojava.nbio.structure.test.align.fatcat;
import junit.framework.TestCase;
+
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.StructureAlignment;
import org.biojava.nbio.structure.align.fatcat.FatCatFlexible;
@@ -34,19 +35,18 @@
import org.biojava.nbio.structure.align.util.AtomCache;
import org.biojava.nbio.structure.align.xml.AFPChainXMLConverter;
import org.biojava.nbio.structure.align.xml.AFPChainXMLParser;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
-
+import javax.vecmath.Matrix4d;
public class TestFlexibleRotationMatrices extends TestCase{
-
-
-
public void testFlexibleRotationMatrices(){
String name1 = "1a21.A";
@@ -200,7 +200,7 @@ private void compareBlock(int blockNr, AFPChain afpChain, AFPChain newChain,
// calc RMSD
- double rmsdFile = SVDSuperimposer.getRMS(blockSet1, blockSet2);
+ double rmsdFile = Calc.rmsd(blockSet1, blockSet2);
// this is the value from the file. it never seems to match precisely, probably is calculated from initial block.
// we can't reproduce the initial block, since we don;t serialize it.
@@ -209,10 +209,12 @@ private void compareBlock(int blockNr, AFPChain afpChain, AFPChain newChain,
// THIS IS CALCULATING THE "correct" rotation matrix, that should be in the file
- SVDSuperimposer svd = new SVDSuperimposer(blockSet1, blockSet2copy);
+ Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(blockSet1),
+ Calc.atomsToPoints(blockSet2copy));
+
//double rmsdForce = SVDSuperimposer.getRMS(atomSet1, atomSet2);
- Matrix m = svd.getRotation();
- Atom s = svd.getTranslation();
+ Matrix m = Matrices.getRotationJAMA(transform);
+ Atom s = Calc.getTranslationVector(transform);
Matrix max = maxs2[blockNr];
Atom shift = shifts2[blockNr];
@@ -231,7 +233,7 @@ private void compareBlock(int blockNr, AFPChain afpChain, AFPChain newChain,
Calc.rotate(a, m);
Calc.shift( a, s);
}
- double rmsd3 = SVDSuperimposer.getRMS(blockSet1,blockSet2copy);
+ double rmsd3 = Calc.rmsd(blockSet1,blockSet2copy);
assertTrue("The RMSD values don;t match after rotation / shift for block " + blockNr + "! should be: " + rmsd3 + " but found: " +rmsdFile, compareRmsd(rmsd3, rmsdFile));
diff --git a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/symmetry/TestQuatSymmetryDetectorExamples.java b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/symmetry/TestQuatSymmetryDetectorExamples.java
index 9ec0b15b1a..fc3e5d009b 100644
--- a/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/symmetry/TestQuatSymmetryDetectorExamples.java
+++ b/biojava-integrationtest/src/test/java/org/biojava/nbio/structure/test/symmetry/TestQuatSymmetryDetectorExamples.java
@@ -68,9 +68,30 @@ public void testNMR() throws IOException, StructureException {
assertFalse(symmetry.isPseudoStoichiometric());
}
+
+ /**
+ * Test a dihedral symmetry: 2VML
+ */
+ @Test
+ public void testDihedral() throws IOException, StructureException {
+
+ Structure pdb = StructureIO.getStructure("BIO:2vml:1");
+
+ SubunitClustererParameters clusterParams = new SubunitClustererParameters();
+ clusterParams.setClustererMethod(SubunitClustererMethod.SEQUENCE);
+ QuatSymmetryParameters symmParams = new QuatSymmetryParameters();
+ QuatSymmetryResults symmetry = QuatSymmetryDetector.calcGlobalSymmetry(
+ pdb, symmParams, clusterParams);
+
+ // D3 symmetry non pseudosymmetric
+ assertEquals("D3", symmetry.getSymmetry());
+ assertEquals("A6B6", symmetry.getStoichiometry());
+
+ }
/**
* Hemoglobin has both symmetry and pseudosymmetry: 4HHB
+ * TODO pseudosymmetry in hemoglobin fails
*
* @throws StructureException
* @throws IOException
diff --git a/biojava-structure-gui/src/main/java/demo/DemoQuatSymmetryJmol.java b/biojava-structure-gui/src/main/java/demo/DemoQuatSymmetryJmol.java
index 909240eb97..47865e88d8 100644
--- a/biojava-structure-gui/src/main/java/demo/DemoQuatSymmetryJmol.java
+++ b/biojava-structure-gui/src/main/java/demo/DemoQuatSymmetryJmol.java
@@ -61,7 +61,7 @@ public class DemoQuatSymmetryJmol {
public static void main(String[] args) throws IOException,
StructureException {
- String name = "1VYM";
+ String name = "2vml";
// Download the biological assembly
AtomCache cache = new AtomCache();
@@ -72,7 +72,6 @@ public static void main(String[] args) throws IOException,
SubunitClustererParameters cp = new SubunitClustererParameters();
cp.setClustererMethod(SubunitClustererMethod.SEQUENCE); // normal
// cp.setClustererMethod(SubunitClustererMethod.STRUCTURE); // pseudo
- // cp.setClustererMethod(SubunitClustererMethod.INTERNAL_SYMMETRY);
cp.setCoverageThreshold(0.9);
// Calculate and display the global symmetry
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/Calc.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/Calc.java
index 92d90088d6..829cb47589 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/Calc.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/Calc.java
@@ -23,6 +23,8 @@
package org.biojava.nbio.structure;
import org.biojava.nbio.structure.geometry.CalcPoint;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositionSVD;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -992,13 +994,14 @@ public static final Atom createVirtualCBAtom(AminoAcid amino)
arr2[1] = amino.getCA();
arr2[2] = amino.getC();
- // ok now we got the two arrays, do a SVD:
+ // ok now we got the two arrays, do a Superposition:
- SVDSuperimposer svd = new SVDSuperimposer(arr2, arr1);
-
- Matrix rotMatrix = svd.getRotation();
- Atom tranMatrix = svd.getTranslation();
+ SuperPositionSVD svd = new SuperPositionSVD(false);
+ Matrix4d transform = svd.superpose(arr1, arr2);
+ Matrix rotMatrix = Matrices.getRotationJAMA(transform);
+ Atom tranMatrix = getTranslationVector(transform);
+
Calc.rotate(aCB, rotMatrix);
Atom virtualCB = Calc.add(aCB, tranMatrix);
@@ -1180,23 +1183,6 @@ public static void shift(Atom[] ca, Atom b) {
Calc.shift(atom, b);
}
- /**
- * Convert JAMA rotation and translation to a Vecmath transformation matrix.
- * Because the JAMA matrix is a pre-multiplication matrix and the Vecmath
- * matrix is a post-multiplication one, the rotation matrix is transposed to
- * ensure that the transformation they produce is the same.
- *
- * @param rot
- * 3x3 Rotation matrix
- * @param trans
- * 3x1 Translation matrix
- * @return 4x4 transformation matrix
- */
- public static Matrix4d getTransformation(Matrix rot, Matrix trans) {
- return new Matrix4d(new Matrix3d(rot.getColumnPackedCopy()),
- new Vector3d(trans.getColumnPackedCopy()), 1.0);
- }
-
/**
* Convert JAMA rotation and translation to a Vecmath transformation matrix.
* Because the JAMA matrix is a pre-multiplication matrix and the Vecmath
@@ -1213,30 +1199,9 @@ public static Matrix4d getTransformation(Matrix rot, Atom trans) {
return new Matrix4d(new Matrix3d(rot.getColumnPackedCopy()),
new Vector3d(trans.getCoords()), 1.0);
}
-
- /**
- * Convert Vecmath transformation into a JAMA rotation matrix. Because the
- * JAMA matrix is a pre-multiplication matrix and the Vecmath matrix is a
- * post-multiplication one, the rotation matrix is transposed to ensure that
- * the transformation they produce is the same.
- *
- * @param transform
- * Matrix4d with transposed rotation matrix
- * @return
- */
- public static Matrix getRotationMatrix(Matrix4d transform) {
-
- Matrix rot = new Matrix(3, 3);
- for (int i = 0; i < 3; i++) {
- for (int j = 0; j < 3; j++) {
- rot.set(j, i, transform.getElement(i, j)); // transposed
- }
- }
- return rot;
- }
-
+
/**
- * Extract the translational vector of a Vecmath transformation.
+ * Extract the translational vector as an Atom of a transformation matrix.
*
* @param transform
* Matrix4d
@@ -1264,4 +1229,127 @@ public static Point3d[] atomsToPoints(Atom[] atoms) {
}
return points;
}
+
+ /**
+ * Calculate the RMSD of two Atom arrays, already superposed.
+ *
+ * @param x
+ * array of Atoms superposed to y
+ * @param y
+ * array of Atoms superposed to x
+ * @return RMSD
+ */
+ public static double rmsd(Atom[] x, Atom[] y) {
+ return CalcPoint.rmsd(atomsToPoints(x), atomsToPoints(y));
+ }
+
+ /**
+ * Calculate the TM-Score for the superposition.
+ *
+ * Normalizes by the maximum-length structure (that is, {@code max\{len1,len2\}}) rather than the minimum.
+ *
+ * Atom sets must be superposed.
+ *
+ *
+ * Citation:
+ * Zhang Y and Skolnick J (2004). "Scoring function for automated
+ * assessment of protein structure template quality". Proteins 57: 702 -
+ * 710.
+ *
+ * @param atomSet1
+ * atom array 1
+ * @param atomSet2
+ * atom array 2
+ * @param len1
+ * The full length of the protein supplying atomSet1
+ * @param len2
+ * The full length of the protein supplying atomSet2
+ * @return The TM-Score
+ * @throws StructureException
+ * @see {@link #getTMScore(Atom[], Atom[], int, int)}, which normalizes by
+ * the minimum length
+ */
+ public static double getTMScoreAlternate(Atom[] atomSet1, Atom[] atomSet2,
+ int len1, int len2) throws StructureException {
+ if (atomSet1.length != atomSet2.length) {
+ throw new StructureException(
+ "The two atom sets are not of same length!");
+ }
+ if (atomSet1.length > len1) {
+ throw new StructureException(
+ "len1 must be greater or equal to the alignment length!");
+ }
+ if (atomSet2.length > len2) {
+ throw new StructureException(
+ "len2 must be greater or equal to the alignment length!");
+ }
+
+ int Lmax = Math.max(len1, len2);
+ int Laln = atomSet1.length;
+
+ double d0 = 1.24 * Math.cbrt(Lmax - 15.) - 1.8;
+ double d0sq = d0 * d0;
+
+ double sum = 0;
+ for (int i = 0; i < Laln; i++) {
+ double d = Calc.getDistance(atomSet1[i], atomSet2[i]);
+ sum += 1. / (1 + d * d / d0sq);
+ }
+
+ return sum / Lmax;
+ }
+
+ /**
+ * Calculate the TM-Score for the superposition.
+ *
+ * Normalizes by the minimum-length structure (that is, {@code min\{len1,len2\}}).
+ *
+ * Atom sets must be pre-rotated.
+ *
+ *
+ * Citation:
+ * Zhang Y and Skolnick J (2004). "Scoring function for automated
+ * assessment of protein structure template quality". Proteins 57: 702 -
+ * 710.
+ *
+ * @param atomSet1
+ * atom array 1
+ * @param atomSet2
+ * atom array 2
+ * @param len1
+ * The full length of the protein supplying atomSet1
+ * @param len2
+ * The full length of the protein supplying atomSet2
+ * @return The TM-Score
+ * @throws StructureException
+ */
+ public static double getTMScore(Atom[] atomSet1, Atom[] atomSet2, int len1,
+ int len2) throws StructureException {
+ if (atomSet1.length != atomSet2.length) {
+ throw new StructureException(
+ "The two atom sets are not of same length!");
+ }
+ if (atomSet1.length > len1) {
+ throw new StructureException(
+ "len1 must be greater or equal to the alignment length!");
+ }
+ if (atomSet2.length > len2) {
+ throw new StructureException(
+ "len2 must be greater or equal to the alignment length!");
+ }
+
+ int Lmin = Math.min(len1, len2);
+ int Laln = atomSet1.length;
+
+ double d0 = 1.24 * Math.cbrt(Lmin - 15.) - 1.8;
+ double d0sq = d0 * d0;
+
+ double sum = 0;
+ for (int i = 0; i < Laln; i++) {
+ double d = Calc.getDistance(atomSet1[i], atomSet2[i]);
+ sum += 1. / (1 + d * d / d0sq);
+ }
+
+ return sum / Lmin;
+ }
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/SVDSuperimposer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/SVDSuperimposer.java
deleted file mode 100644
index 6bd953dc19..0000000000
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/SVDSuperimposer.java
+++ /dev/null
@@ -1,428 +0,0 @@
-/*
- * BioJava development 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.
- *
- * For more information on the BioJava project and its aims,
- * or to join the biojava-l mailing list, visit the home page
- * at:
- *
- * http://www.biojava.org/
- *
- * Created on Dec 4, 2005
- *
- */
-package org.biojava.nbio.structure;
-
-import javax.vecmath.Matrix4d;
-import javax.vecmath.Point3d;
-
-import org.biojava.nbio.structure.geometry.CalcPoint;
-import org.biojava.nbio.structure.jama.Matrix;
-import org.biojava.nbio.structure.jama.SingularValueDecomposition;
-
-
-/** A class that calculates the superimposition between two sets of atoms
- * inspired by the biopython SVDSuperimposer class...
- *
- *
- * example usage:
- *
- // get some arbitrary amino acids from somewhere
- String filename = "/Users/ap3/WORK/PDB/5pti.pdb" ;
-
- PDBFileReader pdbreader = new PDBFileReader();
- Structure struc = pdbreader.getStructure(filename);
- Group g1 = (Group)struc.getChain(0).getGroup(21).clone();
- Group g2 = (Group)struc.getChain(0).getGroup(53).clone();
-
- if ( g1.getPDBName().equals("GLY")){
- if ( g1 instanceof AminoAcid){
- Atom cb = Calc.createVirtualCBAtom((AminoAcid)g1);
- g1.addAtom(cb);
- }
- }
-
- if ( g2.getPDBName().equals("GLY")){
- if ( g2 instanceof AminoAcid){
- Atom cb = Calc.createVirtualCBAtom((AminoAcid)g2);
- g2.addAtom(cb);
- }
- }
-
- Structure struc2 = new StructureImpl((Group)g2.clone());
-
- System.out.println(g1);
- System.out.println(g2);
-
-
- Atom[] atoms1 = new Atom[3];
- Atom[] atoms2 = new Atom[3];
-
- atoms1[0] = g1.getAtom("N");
- atoms1[1] = g1.getAtom("CA");
- atoms1[2] = g1.getAtom("CB");
-
-
- atoms2[0] = g2.getAtom("N");
- atoms2[1] = g2.getAtom("CA");
- atoms2[2] = g2.getAtom("CB");
-
-
- SVDSuperimposer svds = new SVDSuperimposer(atoms1,atoms2);
-
-
- Matrix rotMatrix = svds.getRotation();
- Atom tranMatrix = svds.getTranslation();
-
-
- // now we have all the info to perform the rotations ...
-
- Calc.rotate(struc2,rotMatrix);
-
- // shift structure 2 onto structure one ...
- Calc.shift(struc2,tranMatrix);
-
- //
- // write the whole thing to a file to view in a viewer
-
- String outputfile = "/Users/ap3/WORK/PDB/rotated.pdb";
-
- FileOutputStream out= new FileOutputStream(outputfile);
- PrintStream p = new PrintStream( out );
-
- Structure newstruc = new StructureImpl();
-
- Chain c1 = new ChainImpl();
- c1.setName("A");
- c1.addGroup(g1);
- newstruc.addChain(c1);
-
- Chain c2 = struc2.getChain(0);
- c2.setName("B");
- newstruc.addChain(c2);
-
- // show where the group was originally ...
- Chain c3 = new ChainImpl();
- c3.setName("C");
- //c3.addGroup(g1);
- c3.addGroup(g2);
-
- newstruc.addChain(c3);
- p.println(newstruc.toPDB());
-
- p.close();
-
- System.out.println("wrote to file " + outputfile);
-
-
-
- * TODO the documentaton should contain a reference and explanation
- * of the algorithm, rather than example code.
- * The example should go to demo (Aleix 08/2016)
- *
- * @author Andreas Prlic
- * @since 1.5
- * @version %I% %G%
-
- */
-public class SVDSuperimposer {
-
- private Matrix rot;
- private Matrix tran;
-
- private Matrix centroidA;
- private Matrix centroidB;
-
- /** Create a SVDSuperimposer object and calculate a SVD
- * superimposition of two sets of atoms.
- *
- * @param atomSet1 Atom array 1, fixed
- * @param atomSet2 Atom array 2, moved/transformed
- * @throws StructureException
- */
- public SVDSuperimposer(Atom[] atomSet1,Atom[]atomSet2)
- throws StructureException{
-
- if ( atomSet1.length != atomSet2.length ){
- throw new StructureException("The two atom sets are not of same length!");
- }
-
- Atom cena = Calc.getCentroid(atomSet1);
- Atom cenb = Calc.getCentroid(atomSet2);
-
- double[][] centAcoords = new double[][]{{cena.getX(),cena.getY(),cena.getZ()}};
- centroidA = new Matrix(centAcoords);
-
- double[][] centBcoords = new double[][]{{cenb.getX(),cenb.getY(),cenb.getZ()}};
- centroidB = new Matrix(centBcoords);
-
- // center at centroid
-
- Atom[] ats1 = Calc.centerAtoms(atomSet1,cena);
- Atom[] ats2 = Calc.centerAtoms(atomSet2,cenb);
-
- double[][] coordSet1 = new double[ats1.length][3];
- double[][] coordSet2 = new double[ats2.length][3];
-
- // copy the atoms into the internal coords;
- for (int i =0 ; i< ats1.length;i++) {
- coordSet1[i] = ats1[i].getCoords();
- coordSet2[i] = ats2[i].getCoords();
- }
-
- calculate(coordSet1,coordSet2);
-
- }
-
- /** Create a SVDSuperimposer object and calculate a SVD
- * superimposition of two sets of points.
- *
- * @param pointSet1 Point3d array 1, fixed
- * @param pointSet2 Point3d array 2, moved/ transformed
- * @throws StructureException
- */
- public SVDSuperimposer(Point3d[] pointSet1, Point3d[] pointSet2)
- throws StructureException{
-
- if ( pointSet1.length != pointSet2.length ){
- throw new StructureException("The two point sets are not of same length!");
- }
-
- Point3d cena = CalcPoint.centroid(pointSet1);
- Point3d cenb = CalcPoint.centroid(pointSet2);
-
- double[][] centAcoords = new double[][]{{cena.x,cena.y,cena.z}};
- centroidA = new Matrix(centAcoords);
-
- double[][] centBcoords = new double[][]{{cenb.x,cenb.x,cenb.x}};
- centroidB = new Matrix(centBcoords);
-
- // center at centroid
-
- cena.negate();
- cenb.negate();
-
- Point3d[] ats1 = CalcPoint.clonePoint3dArray(pointSet1);
- CalcPoint.translate(cena, ats1);
-
- Point3d[] ats2 = CalcPoint.clonePoint3dArray(pointSet2);
- CalcPoint.translate(cenb, ats2);
-
- double[][] coordSet1 = new double[ats1.length][3];
- double[][] coordSet2 = new double[ats2.length][3];
-
- // copy the atoms into the internal coords;
- for (int i =0 ; i< ats1.length;i++) {
- coordSet1[i] = new double[3];
- ats1[i].get(coordSet1[i]);
- coordSet2[i] = new double[3];
- ats2[i].get(coordSet2[i]);
- }
-
- calculate(coordSet1,coordSet2);
-
- }
-
-
- /** Do the actual calculation.
- *
- * @param coordSet1 coordinates for atom array 1
- * @param coordSet2 coordiantes for atom array 2
- */
- private void calculate(double[][] coordSet1, double[][]coordSet2){
- // now this is the bridge to the Jama package:
- Matrix a = new Matrix(coordSet1);
- Matrix b = new Matrix(coordSet2);
-
-
-// # correlation matrix
-
- Matrix b_trans = b.transpose();
- Matrix corr = b_trans.times(a);
-
-
- SingularValueDecomposition svd = corr.svd();
-
- Matrix u = svd.getU();
- // v is alreaady transposed ! difference to numermic python ...
- Matrix vt =svd.getV();
-
- Matrix vt_orig = (Matrix) vt.clone();
- Matrix u_transp = u.transpose();
-
- Matrix rot_nottrans = vt.times(u_transp);
- rot = rot_nottrans.transpose();
-
- // check if we have found a reflection
-
- double det = rot.det();
-
- if (det<0) {
- vt = vt_orig.transpose();
- vt.set(2,0,(0 - vt.get(2,0)));
- vt.set(2,1,(0 - vt.get(2,1)));
- vt.set(2,2,(0 - vt.get(2,2)));
-
- Matrix nv_transp = vt.transpose();
- rot_nottrans = nv_transp.times(u_transp);
- rot = rot_nottrans.transpose();
-
- }
-
- Matrix cb_tmp = centroidB.times(rot);
- tran = centroidA.minus(cb_tmp);
-
-
- }
-
- /** Calculate the RMS (root mean square) deviation of two sets of atoms.
- *
- * Atom sets must be pre-rotated.
- *
- * @param atomSet1 atom array 1
- * @param atomSet2 atom array 2
- * @return the RMS of two atom sets
- * @throws StructureException
- */
- public static double getRMS(Atom[] atomSet1, Atom[] atomSet2) throws StructureException {
- if ( atomSet1.length != atomSet2.length ){
- throw new StructureException("The two atom sets are not of same length!");
- }
-
- double sum = 0.0;
- for ( int i =0 ; i < atomSet1.length;i++){
- double d = Calc.getDistance(atomSet1[i],atomSet2[i]);
- sum += (d*d);
-
- }
-
- double avd = ( sum/ atomSet1.length);
- //System.out.println("av dist" + avd);
- return Math.sqrt(avd);
-
- }
-
- /**
- * Calculate the TM-Score for the superposition.
- *
- * Normalizes by the minimum-length structure (that is, {@code min\{len1,len2\}}).
- *
- * Atom sets must be pre-rotated.
- *
- * Citation:
- * Zhang Y and Skolnick J (2004). "Scoring function for automated assessment
- * of protein structure template quality". Proteins 57: 702 - 710.
- *
- * @param atomSet1 atom array 1
- * @param atomSet2 atom array 2
- * @param len1 The full length of the protein supplying atomSet1
- * @param len2 The full length of the protein supplying atomSet2
- * @return The TM-Score
- * @throws StructureException
- */
- public static double getTMScore(Atom[] atomSet1, Atom[] atomSet2, int len1, int len2) throws StructureException {
- if ( atomSet1.length != atomSet2.length ){
- throw new StructureException("The two atom sets are not of same length!");
- }
- if ( atomSet1.length > len1 ){
- throw new StructureException("len1 must be greater or equal to the alignment length!");
- }
- if ( atomSet2.length > len2 ){
- throw new StructureException("len2 must be greater or equal to the alignment length!");
- }
-
- int Lmin = Math.min(len1,len2);
- int Laln = atomSet1.length;
-
- double d0 = 1.24 * Math.cbrt(Lmin - 15.) - 1.8;
- double d0sq = d0*d0;
-
- double sum = 0;
- for(int i=0;iNormalizes by the maximum-length structure (that is, {@code max\{len1,len2\}}) rather than the minimum.
- *
- * Atom sets must be pre-rotated.
- *
- * Citation:
- * Zhang Y and Skolnick J (2004). "Scoring function for automated assessment
- * of protein structure template quality". Proteins 57: 702 - 710.
- *
- * @param atomSet1 atom array 1
- * @param atomSet2 atom array 2
- * @param len1 The full length of the protein supplying atomSet1
- * @param len2 The full length of the protein supplying atomSet2
- * @return The TM-Score
- * @throws StructureException
- * @see {@link #getTMScore(Atom[], Atom[], int, int)}, which normalizes by the minimum length
- */
- public static double getTMScoreAlternate(Atom[] atomSet1, Atom[] atomSet2, int len1, int len2) throws StructureException {
- if ( atomSet1.length != atomSet2.length ){
- throw new StructureException("The two atom sets are not of same length!");
- }
- if ( atomSet1.length > len1 ){
- throw new StructureException("len1 must be greater or equal to the alignment length!");
- }
- if ( atomSet2.length > len2 ){
- throw new StructureException("len2 must be greater or equal to the alignment length!");
- }
-
- int Lmax = Math.max(len1,len2);
- int Laln = atomSet1.length;
-
- double d0 = 1.24 * Math.cbrt(Lmax - 15.) - 1.8;
- double d0sq = d0*d0;
-
- double sum = 0;
- for(int i=0;i afpSet = afpChain.getAfpSet();
int blockNum = afpChain.getBlockNum();
- int i, b2, e2;
-
+ int i, b2, e2;
- //superimposing according to the initial AFP-chaining
- Atom[] origCA = StructureTools.cloneAtomArray(ca2);
+ // 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();
+ int[] afpChainList = afpChain.getAfpChainList();
+ int[] block2Afp = afpChain.getBlock2Afp();
+ int[] blockSize = afpChain.getBlockSize();
+ int[] focusAfpList = afpChain.getFocusAfpList();
- if ( focusAfpList == null){
+ if (focusAfpList == null) {
focusAfpList = new int[afpChain.getMinLen()];
afpChain.setFocusAfpList(focusAfpList);
}
- int focusAfpn = 0;
+ int focusAfpn = 0;
e2 = 0;
b2 = 0;
-
logger.debug("blockNUm at twister: ", blockNum);
- for(int bk = 0; bk < blockNum; bk ++) {
+ 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);
+ 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) {
+ // 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)) {
+ if (bk < (blockNum - 1)) {
- //bend at the middle of two consecutive AFPs
+ // bend at the middle of two consecutive AFPs
int afpPos = afpChainList[block2Afp[bk] + blockSize[bk] - 1];
AFP a1 = afpSet.get(afpPos);
- e2 = a1.getP2() ;
+ e2 = a1.getP2();
int afpPos2 = afpChainList[block2Afp[bk + 1]];
AFP a2 = afpSet.get(afpPos2);
- e2 = (a2.getP2() - e2)/ 2 + e2;
+ e2 = (a2.getP2() - e2) / 2 + e2;
logger.debug("e2 : {}", e2);
- } else{
+ } else {
// last one is until the end...
e2 = ca2.length;
}
@@ -122,26 +127,25 @@ public static Group[] twistPDB(AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws
// this copies the coordinates over into iniTwistPdb
cloneAtomRange(iniTwistPdb, ca2, b2, e2);
- //bound[bk] = e2;
- for(i = 0; i < blockSize[bk]; i ++) {
+ // bound[bk] = e2;
+ for (i = 0; i < blockSize[bk]; i++) {
focusAfpList[focusAfpn] = afpChainList[block2Afp[bk] + i];
focusAfpn++;
}
}
-
- int focusResn = afp2Res(afpChain, focusAfpn, focusAfpList, 0);
+ 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() );
-
+ double totalRmsdIni = calCaRmsd(ca1, iniTwistPdb, focusResn,
+ afpChain.getFocusRes1(), afpChain.getFocusRes2());
afpChain.setTotalRmsdIni(totalRmsdIni);
logger.debug("got iniRMSD: {}", totalRmsdIni);
- if ( totalRmsdIni == 5.76611141613097) {
+ if (totalRmsdIni == 5.76611141613097) {
logger.debug("{}", afpChain.getAlnseq1());
logger.debug("{}", afpChain.getAlnsymb());
logger.debug("{}", afpChain.getAlnseq2());
@@ -151,10 +155,11 @@ public static Group[] twistPDB(AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws
afpChain.setBlock2Afp(block2Afp);
afpChain.setAfpChainList(afpChainList);
- return twistOptimized(afpChain,ca1,origCA);
+ return twistOptimized(afpChain, ca1, origCA);
}
- /** superimposing according to the optimized alignment
+ /**
+ * superimposing according to the optimized alignment
*
* @param afpChain
* @param ca1
@@ -162,65 +167,75 @@ public static Group[] twistPDB(AFPChain afpChain,Atom[] ca1, Atom[] ca2) throws
* @return Group array twisted.
* @throws StructureException
*/
- public static Group[] twistOptimized(AFPChain afpChain,Atom[] ca1, Atom[] ca2) 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){
+ for (Atom a : ca2) {
gPos++;
optTwistPdb[gPos] = a;
}
int blockNum = afpChain.getBlockNum();
- int b2 = 0;
- int e2 = 0;
+ 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 (focusRes1 == null) {
+ focusRes1 = new int[afpChain.getCa1Length()];
+ afpChain.setFocusRes1(focusRes1);
}
- if(focusRes2 == null) {
- focusRes2 = new int[afpChain.getCa2Length()];
- afpChain.setFocusRes2(focusRes2);
+ if (focusRes2 == null) {
+ focusRes2 = new int[afpChain.getCa2Length()];
+ afpChain.setFocusRes2(focusRes2);
}
- int[] optLen = afpChain.getOptLen();
- int[][][] optAln = afpChain.getOptAln();
+ int[] optLen = afpChain.getOptLen();
+ int[][][] optAln = afpChain.getOptAln();
- for(int bk = 0; bk < blockNum; bk ++) {
- // THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
+ 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);
+ 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
+ // 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;
+ e2 = (optAln[bk + 1][1][0] - e2) / 2 + e2;
+ } else {
+ e2 = ca2.length;
}
- else { e2 = ca2.length; }
cloneAtomRange(optTwistPdb, ca2, b2, e2);
- for(int i = 0; i < optLen[bk]; i ++) {
+ for (int i = 0; i < optLen[bk]; i++) {
focusRes1[focusResn] = optAln[bk][0][i];
focusRes2[focusResn] = optAln[bk][1][i];
- focusResn ++;
+ focusResn++;
}
}
int totalLenOpt = focusResn;
logger.debug("calrmsdopt for {} residues", focusResn);
- double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1, focusRes2);
+ 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);
+ if (totalLenOpt != optLength) {
+ logger.warn("Final alignment length is different {} {}",
+ totalLenOpt, optLength);
}
- logger.debug("final alignment length {}, rmsd {}", focusResn, totalRmsdOpt);
+ logger.debug("final alignment length {}, rmsd {}", focusResn,
+ totalRmsdOpt);
afpChain.setTotalLenOpt(totalLenOpt);
afpChain.setTotalRmsdOpt(totalRmsdOpt);
@@ -229,42 +244,42 @@ public static Group[] twistOptimized(AFPChain afpChain,Atom[] ca1, Atom[] ca2) t
}
-
/**
- * transform the coordinates in the ca2 according to the superimposing of the given position pairs.
- * No Cloning, transforms input atoms.
+ * 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);
+ 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);
+ 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);
+ // double *cod1 = pro1->Cod4Res(n, res1);
+ // double *cod2 = pro2->Cod4Res(n, res2);
- Matrix r;
- Atom t;
+ Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(cod1),
+ Calc.atomsToPoints(cod2));
- SVDSuperimposer svd = new SVDSuperimposer(cod1, cod2);
+ Matrix r = Matrices.getRotationJAMA(transform);
+ Atom t = Calc.getTranslationVector(transform);
- r = svd.getRotation();
- t = svd.getTranslation();
+ logger.debug("transPdb: transforming orig coordinates with matrix: {}",
+ r);
- logger.debug("transPdb: transforming orig coordinates with matrix: {}", r);
-
- if ( afpChain != null ){
+ if (afpChain != null) {
Matrix[] ms = afpChain.getBlockRotationMatrix();
- if ( ms == null)
+ if (ms == null)
ms = new Matrix[afpChain.getBlockNum()];
- ms[blockNr]= r;
+ ms[blockNr] = r;
Atom[] shifts = afpChain.getBlockShiftVector();
- if ( shifts == null)
+ if (shifts == null)
shifts = new Atom[afpChain.getBlockNum()];
shifts[blockNr] = t;
@@ -272,27 +287,24 @@ private static void transformOrigPDB(int n, int[] res1, int[] res2, Atom[] ca1,
afpChain.setBlockShiftVector(shifts);
}
- for (Atom a : ca2){
- Calc.rotate(a.getGroup(),r);
- Calc.shift(a.getGroup(),t);
- }
-
+ for (Atom a : ca2)
+ Calc.transform(a.getGroup(), transform);
}
// like Cod4Res
// most likely the clone flag is not needed
- private static Atom[] getAtoms(Atom[] ca, int[] positions, int length, boolean clone){
+ private static Atom[] getAtoms(Atom[] ca, int[] positions, int length,
+ boolean clone) {
List atoms = new ArrayList();
- for ( int i = 0 ; i < length ; i++){
+ 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 {
+ if (clone) {
+ a = (Atom) ca[p].clone();
+ a.setGroup((Group) ca[p].getGroup().clone());
+ } else {
a = ca[p];
}
atoms.add(a);
@@ -300,9 +312,8 @@ private static Atom[] getAtoms(Atom[] ca, int[] positions, int length, boolean c
return atoms.toArray(new Atom[atoms.size()]);
}
-
-
- /** Clones and copies the Atoms from p2 into p1 range is between r1 and r2
+ /**
+ * Clones and copies the Atoms from p2 into p1 range is between r1 and r2
*
* @param p1
* @param p2
@@ -311,32 +322,32 @@ private static Atom[] getAtoms(Atom[] ca, int[] positions, int length, boolean c
*/
// orig name: modifyCod
private static void cloneAtomRange(Atom[] p1, Atom[] p2, int r1, int r2)
- throws StructureException
- {
+ throws StructureException {
logger.debug("modifyCod from: {} to: {}", r1, r2);
- // special clone method, can;t use StructureTools.cloneCAArray, since we access the data
+ // 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 ++) {
+ for (int i = r1; i < r2; i++) {
Group g = p2[i].getGroup();
- Group newG = (Group)g.clone();
+ 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.getName().equals(parentC.getName())){
+ for (Chain c : model) {
+ if (c.getName().equals(parentC.getName())) {
newChain = c;
break;
}
}
- if ( newChain == null){
+ if (newChain == null) {
newChain = new ChainImpl();
newChain.setName(parentC.getName());
model.add(newChain);
@@ -344,81 +355,75 @@ private static void cloneAtomRange(Atom[] p1, Atom[] p2, int r1, int r2)
newChain.addGroup(newG);
-
- } //modify caCod
+ } // 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.
+ * 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).
+ * 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
- {
+ 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);
-
- 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){
+ 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();
+ Matrix4d transform = SuperPositions.superpose(Calc.atomsToPoints(cod1),
+ Calc.atomsToPoints(cod2));
- for (Atom a : cod2){
+ for (Atom a : cod2)
+ Calc.transform(a.getGroup(), transform);
- Calc.rotate(a.getGroup(), r);
- Calc.shift(a.getGroup(), t);
- }
-
- return SVDSuperimposer.getRMS(cod1, cod2);
+ return Calc.rmsd(cod1, cod2);
}
/**
- * Set the list of equivalent residues in the two proteins given a list of AFPs
+ * 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!
+ * WARNING: changes the values for FocusRes1, focusRes2 and FocusResn in
+ * afpChain!
*
- * @param afpChain the AFPChain to store resuts
- * @param afpn nr of afp
+ * @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 )
- {
+ 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();
+ List afpSet = afpChain.getAfpSet();
- for(int i = listStart; i < listStart+afpn; i ++) {
+ 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!");
+ 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 ++;
+ n++;
}
}
@@ -426,8 +431,9 @@ public static int afp2Res(AFPChain afpChain, int afpn, int[] afpPositions , int
afpChain.setFocusRes2(res2);
afpChain.setFocusResn(n);
- if ( n == 0 ){
- logger.warn("n=0!!! + " + afpn + " " + listStart + " " + afpPositions.length);
+ if (n == 0) {
+ logger.warn("n=0!!! + " + afpn + " " + listStart + " "
+ + afpPositions.length);
}
return n;
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/StructurePairAligner.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/StructurePairAligner.java
index ea7e238c8d..3abf06c4cf 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/StructurePairAligner.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/StructurePairAligner.java
@@ -27,6 +27,8 @@
import org.biojava.nbio.structure.align.helper.AlignTools;
import org.biojava.nbio.structure.align.helper.JointFragments;
import org.biojava.nbio.structure.align.pairwise.*;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.io.PDBFileParser;
import org.biojava.nbio.structure.io.PDBFileReader;
import org.biojava.nbio.structure.jama.Matrix;
@@ -41,86 +43,88 @@
import java.util.Comparator;
import java.util.List;
+import javax.vecmath.Matrix4d;
/**
* Perform a pairwise protein structure superimposition.
*
*
- * The algorithm is a distance matrix based, rigid body protein structure superimposition.
- * It is based on a variation of the PSC++ algorithm provided by Peter Lackner
- * (Peter.Lackner@sbg.ac.at, personal communication) .
+ * The algorithm is a distance matrix based, rigid body protein structure
+ * superimposition. It is based on a variation of the PSC++ algorithm provided
+ * by Peter Lackner (Peter.Lackner@sbg.ac.at, personal communication) .
*
*
*
*
* Example
- *
+ *
+ *
* public void run(){
-
- // first load two example structures
- {@link InputStream} inStream1 = this.getClass().getResourceAsStream("/files/5pti.pdb");
- {@link InputStream} inStream2 = this.getClass().getResourceAsStream("/files/1tap.pdb");
-
- {@link Structure} structure1 = null;
- {@link Structure} structure2 = null;
-
- {@link PDBFileParser} pdbpars = new {@link PDBFileParser}();
- structure1 = pdbpars.parsePDBFile(inStream1) ;
- structure2 = pdbpars.parsePDBFile(inStream2);
-
-
- // calculate structure superimposition for two complete structures
- {@link StructurePairAligner} aligner = new {@link StructurePairAligner}();
-
-
- // align the full 2 structures with default parameters.
- // see StructurePairAligner for more options and how to align
- // any set of Atoms
- aligner.align(structure1,structure2);
-
- {@link AlternativeAlignment}[] aligs = aligner.getAlignments();
- {@link AlternativeAlignment} a = aligs[0];
- System.out.println(a);
-
- //display the alignment in Jmol
-
- // first get an artificial structure for the alignment
- {@link Structure} artificial = a.getAlignedStructure(structure1, structure2);
-
-
- // and then send it to Jmol (only will work if Jmol is in the Classpath)
-
- BiojavaJmol jmol = new BiojavaJmol();
- jmol.setTitle(artificial.getName());
- jmol.setStructure(artificial);
-
- // color the two structures
-
-
- jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " +
- "select not protein and not solvent; spacefill on;");
- jmol.evalString("select *"+"/1 ; color red; model 1; ");
-
-
- // now color the equivalent residues ...
-
- String[] pdb1 = a.getPDBresnum1();
- for (String res : pdb1 ){
- jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;");
- }
-
- jmol.evalString("select *"+"/2; color blue; model 2;");
- String[] pdb2 = a.getPDBresnum2();
- for (String res :pdb2 ){
- jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;");
- }
-
-
- // now show both models again.
- jmol.evalString("model 0;");
-
- }
- *
+ *
+ * // first load two example structures
+ * {@link InputStream} inStream1 = this.getClass().getResourceAsStream("/files/5pti.pdb");
+ * {@link InputStream} inStream2 = this.getClass().getResourceAsStream("/files/1tap.pdb");
+ *
+ * {@link Structure} structure1 = null;
+ * {@link Structure} structure2 = null;
+ *
+ * {@link PDBFileParser} pdbpars = new {@link PDBFileParser}();
+ * structure1 = pdbpars.parsePDBFile(inStream1) ;
+ * structure2 = pdbpars.parsePDBFile(inStream2);
+ *
+ *
+ * // calculate structure superimposition for two complete structures
+ * {@link StructurePairAligner} aligner = new {@link StructurePairAligner}();
+ *
+ *
+ * // align the full 2 structures with default parameters.
+ * // see StructurePairAligner for more options and how to align
+ * // any set of Atoms
+ * aligner.align(structure1,structure2);
+ *
+ * {@link AlternativeAlignment}[] aligs = aligner.getAlignments();
+ * {@link AlternativeAlignment} a = aligs[0];
+ * System.out.println(a);
+ *
+ * //display the alignment in Jmol
+ *
+ * // first get an artificial structure for the alignment
+ * {@link Structure} artificial = a.getAlignedStructure(structure1, structure2);
+ *
+ *
+ * // and then send it to Jmol (only will work if Jmol is in the Classpath)
+ *
+ * BiojavaJmol jmol = new BiojavaJmol();
+ * jmol.setTitle(artificial.getName());
+ * jmol.setStructure(artificial);
+ *
+ * // color the two structures
+ *
+ *
+ * jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " +
+ * "select not protein and not solvent; spacefill on;");
+ * jmol.evalString("select *"+"/1 ; color red; model 1; ");
+ *
+ *
+ * // now color the equivalent residues ...
+ *
+ * String[] pdb1 = a.getPDBresnum1();
+ * for (String res : pdb1 ){
+ * jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;");
+ * }
+ *
+ * jmol.evalString("select *"+"/2; color blue; model 2;");
+ * String[] pdb2 = a.getPDBresnum2();
+ * for (String res :pdb2 ){
+ * jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;");
+ * }
+ *
+ *
+ * // now show both models again.
+ * jmol.evalString("model 0;");
+ *
+ * }
+ *
*
*
*
@@ -131,12 +135,13 @@
*/
public class StructurePairAligner {
- private final static Logger logger = LoggerFactory.getLogger(StructurePairAligner.class);
+ private final static Logger logger = LoggerFactory
+ .getLogger(StructurePairAligner.class);
AlternativeAlignment[] alts;
Matrix distanceMatrix;
StrucAligParameters params;
- FragmentPair[] fragPairs;
+ FragmentPair[] fragPairs;
List listeners = new ArrayList();
@@ -145,19 +150,19 @@ public StructurePairAligner() {
params = StrucAligParameters.getDefaultParameters();
reset();
alts = new AlternativeAlignment[0];
- distanceMatrix = new Matrix(0,0);
+ distanceMatrix = new Matrix(0, 0);
}
- public void addProgressListener(AlignmentProgressListener li){
+ public void addProgressListener(AlignmentProgressListener li) {
listeners.add(li);
}
- public void clearListeners(){
+ public void clearListeners() {
listeners.clear();
}
-
- /** example usage of this class
+ /**
+ * example usage of this class
*
* @param args
*/
@@ -167,19 +172,17 @@ public static void main(String[] args) throws Exception {
PDBFileReader pdbr = new PDBFileReader();
pdbr.setPath("/Users/andreas/WORK/PDB/");
-
- //String pdb1 = "1crl";
- //String pdb2 = "1ede";
+ // String pdb1 = "1crl";
+ // String pdb2 = "1ede";
String pdb1 = "1buz";
String pdb2 = "1ali";
- String outputfile = "/tmp/alig_"+pdb1+"_"+pdb2+".pdb";
+ String outputfile = "/tmp/alig_" + pdb1 + "_" + pdb2 + ".pdb";
// NO NEED TO DO CHANGE ANYTHING BELOW HERE...
StructurePairAligner sc = new StructurePairAligner();
-
// step1 : read molecules
logger.info("aligning {} vs. {}", pdb1, pdb2);
@@ -188,65 +191,61 @@ public static void main(String[] args) throws Exception {
Structure s2 = pdbr.getStructureById(pdb2);
// step 2 : do the calculations
- sc.align(s1,s2);
-
+ sc.align(s1, s2);
AlternativeAlignment[] aligs = sc.getAlignments();
- //cluster similar results together
+ // cluster similar results together
ClusterAltAligs.cluster(aligs);
-
// print the result:
- // the AlternativeAlignment object gives also access to rotation matrices / shift vectors.
+ // the AlternativeAlignment object gives also access to rotation
+ // matrices / shift vectors.
for (AlternativeAlignment aa : aligs) {
logger.info("Alternative Alignment: ", aa);
}
+ // convert AlternativeAlignemnt 1 to PDB file, so it can be opened with
+ // a viewer (e.g. Jmol, Rasmol)
-
- // convert AlternativeAlignemnt 1 to PDB file, so it can be opened with a viewer (e.g. Jmol, Rasmol)
-
- if ( aligs.length > 0) {
- AlternativeAlignment aa1 =aligs[0];
- String pdbstr = aa1.toPDB(s1,s2);
+ if (aligs.length > 0) {
+ AlternativeAlignment aa1 = aligs[0];
+ String pdbstr = aa1.toPDB(s1, s2);
logger.info("writing alignment to {}", outputfile);
- FileOutputStream out= new FileOutputStream(outputfile);
- PrintStream p = new PrintStream( out );
+ FileOutputStream out = new FileOutputStream(outputfile);
+ PrintStream p = new PrintStream(out);
- p.println (pdbstr);
+ p.println(pdbstr);
p.close();
out.close();
}
-
// display the alignment in Jmol
// only will work if Jmol is in the Classpath
- if ( aligs.length > 0) {
+ if (aligs.length > 0) {
- if (! GuiWrapper.isGuiModuleInstalled()){
+ if (!GuiWrapper.isGuiModuleInstalled()) {
logger.error("Could not find structure-gui modules in classpath, please install first!");
}
}
-
}
-
-
- private void reset(){
+ private void reset() {
alts = new AlternativeAlignment[0];
- distanceMatrix = new Matrix(0,0);
+ distanceMatrix = new Matrix(0, 0);
fragPairs = new FragmentPair[0];
}
-
- /** get the results of step 1 - the FragmentPairs used for seeding the alignment
+ /**
+ * get the results of step 1 - the FragmentPairs used for seeding the
+ * alignment
+ *
* @return a FragmentPair[] array
*/
@@ -254,14 +253,13 @@ public FragmentPair[] getFragmentPairs() {
return fragPairs;
}
-
-
public void setFragmentPairs(FragmentPair[] fragPairs) {
this.fragPairs = fragPairs;
}
-
- /** return the alternative alignments that can be found for the two structures
+ /**
+ * return the alternative alignments that can be found for the two
+ * structures
*
* @return AlternativeAlignment[] array
*/
@@ -269,15 +267,17 @@ public AlternativeAlignment[] getAlignments() {
return alts;
}
- /** return the difference of distance matrix between the two structures
+ /**
+ * return the difference of distance matrix between the two structures
*
* @return a Matrix
*/
- public Matrix getDistMat(){
+ public Matrix getDistMat() {
return distanceMatrix;
}
- /** get the parameters.
+ /**
+ * get the parameters.
*
* @return the Parameters.
*/
@@ -285,16 +285,18 @@ public StrucAligParameters getParams() {
return params;
}
- /** set the parameters to be used for the algorithm
+ /**
+ * set the parameters to be used for the algorithm
*
- * @param params the Parameter object
+ * @param params
+ * the Parameter object
*/
public void setParams(StrucAligParameters params) {
this.params = params;
}
-
- /** check if debug mode is set on
+ /**
+ * check if debug mode is set on
*
* @return debug flag
*/
@@ -303,27 +305,32 @@ public boolean isDebug() {
return false;
}
- /** set the debug flag
+ /**
+ * set the debug flag
*
- * @param debug flag
+ * @param debug
+ * flag
*/
@Deprecated
public void setDebug(boolean debug) {
}
- /** Calculate the alignment between the two full structures with default parameters
+ /**
+ * Calculate the alignment between the two full structures with default
+ * parameters
*
* @param s1
* @param s2
* @throws StructureException
*/
- public void align(Structure s1, Structure s2)
- throws StructureException {
+ public void align(Structure s1, Structure s2) throws StructureException {
- align(s1,s2,params);
+ align(s1, s2, params);
}
- /** Calculate the alignment between the two full structures with user provided parameters
+ /**
+ * Calculate the alignment between the two full structures with user
+ * provided parameters
*
* @param s1
* @param s2
@@ -334,28 +341,28 @@ public void align(Structure s1, Structure s2, StrucAligParameters params)
throws StructureException {
// step 1 convert the structures to Atom Arrays
-
Atom[] ca1 = getAlignmentAtoms(s1);
Atom[] ca2 = getAlignmentAtoms(s2);
- notifyStartingAlignment(s1.getName(),ca1,s2.getName(),ca2);
- align(ca1,ca2,params);
+ notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
+ align(ca1, ca2, params);
}
-
- /** Align two chains from the structures. Uses default parameters.
+ /**
+ * Align two chains from the structures. Uses default parameters.
*
* @param s1
* @param chainId1
* @param s2
* @param chainId2
*/
- public void align(Structure s1, String chainId1, Structure s2, String chainId2)
- throws StructureException{
- align(s1,chainId1,s2,chainId2, params);
+ public void align(Structure s1, String chainId1, Structure s2,
+ String chainId2) throws StructureException {
+ align(s1, chainId1, s2, chainId2, params);
}
- /** Aligns two chains from the structures using user provided parameters.
+ /**
+ * Aligns two chains from the structures using user provided parameters.
*
* @param s1
* @param chainId1
@@ -364,8 +371,9 @@ public void align(Structure s1, String chainId1, Structure s2, String chainId2)
* @param params
* @throws StructureException
*/
- public void align(Structure s1, String chainId1, Structure s2, String chainId2, StrucAligParameters params)
- throws StructureException{
+ public void align(Structure s1, String chainId1, Structure s2,
+ String chainId2, StrucAligParameters params)
+ throws StructureException {
reset();
this.params = params;
@@ -381,128 +389,138 @@ public void align(Structure s1, String chainId1, Structure s2, String chainId2,
Atom[] ca1 = getAlignmentAtoms(s3);
Atom[] ca2 = getAlignmentAtoms(s4);
- notifyStartingAlignment(s1.getName(),ca1,s2.getName(),ca2);
- align(ca1,ca2,params);
+ notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
+ align(ca1, ca2, params);
}
- /** Returns the atoms that are being used for the alignment. (E.g. Calpha only, etc.)
+ /**
+ * Returns the atoms that are being used for the alignment. (E.g. Calpha
+ * only, etc.)
*
* @param s
* @return an array of Atoms objects
*/
- public Atom[] getAlignmentAtoms(Structure s){
+ public Atom[] getAlignmentAtoms(Structure s) {
String[] atomNames = params.getUsedAtomNames();
- return StructureTools.getAtomArray(s,atomNames);
+ return StructureTools.getAtomArray(s, atomNames);
}
- /** calculate the protein structure superimposition, between two sets of atoms.
+ /**
+ * calculate the protein structure superimposition, between two sets of
+ * atoms.
*
*
*
- * @param ca1 set of Atoms of structure 1
- * @param ca2 set of Atoms of structure 2
- * @param params the parameters to use for the alignment
+ * @param ca1
+ * set of Atoms of structure 1
+ * @param ca2
+ * set of Atoms of structure 2
+ * @param params
+ * the parameters to use for the alignment
* @throws StructureException
*/
public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
throws StructureException {
-
reset();
this.params = params;
long timeStart = System.currentTimeMillis();
-// step 1 get all Diagonals of length X that are similar between both structures
+ // step 1 get all Diagonals of length X that are similar between both
+ // structures
logger.debug(" length atoms1:" + ca1.length);
logger.debug(" length atoms2:" + ca2.length);
logger.debug("step 1 - get fragments with similar intramolecular distances ");
- int k = params.getDiagonalDistance();
+ int k = params.getDiagonalDistance();
int k2 = params.getDiagonalDistance2();
int fragmentLength = params.getFragmentLength();
- if ( ca1.length < (fragmentLength + 1) ) {
- throw new StructureException("structure 1 too short ("+ca1.length+"), can not align");
+ if (ca1.length < (fragmentLength + 1)) {
+ throw new StructureException("structure 1 too short (" + ca1.length
+ + "), can not align");
}
- if ( ca2.length < (fragmentLength + 1) ){
- throw new StructureException("structure 2 too short ("+ca2.length+"), can not align");
+ if (ca2.length < (fragmentLength + 1)) {
+ throw new StructureException("structure 2 too short (" + ca2.length
+ + "), can not align");
}
int rows = ca1.length - fragmentLength + 1;
int cols = ca2.length - fragmentLength + 1;
- distanceMatrix = new Matrix(rows,cols,0.0);
+ distanceMatrix = new Matrix(rows, cols, 0.0);
double[] dist1 = AlignTools.getDiagonalAtK(ca1, k);
double[] dist2 = AlignTools.getDiagonalAtK(ca2, k);
double[] dist3 = new double[0];
double[] dist4 = new double[0];
- if ( k2 > 0) {
+ if (k2 > 0) {
dist3 = AlignTools.getDiagonalAtK(ca1, k2);
dist4 = AlignTools.getDiagonalAtK(ca2, k2);
}
- double[][] utmp = new double[][] {{0,0,1}};
+ double[][] utmp = new double[][] { { 0, 0, 1 } };
Atom unitvector = new AtomImpl();
unitvector.setCoords(utmp[0]);
List fragments = new ArrayList();
- for ( int i = 0 ; i< rows; i++){
+ for (int i = 0; i < rows; i++) {
- Atom[] catmp1 = AlignTools.getFragment( ca1, i, fragmentLength);
- Atom center1 = AlignTools.getCenter( ca1, i, fragmentLength);
+ Atom[] catmp1 = AlignTools.getFragment(ca1, i, fragmentLength);
+ Atom center1 = AlignTools.getCenter(ca1, i, fragmentLength);
- for ( int j = 0 ; j < cols ; j++){
+ for (int j = 0; j < cols; j++) {
- double rdd1 = AlignTools.rms_dk_diag(dist1,dist2,i,j,fragmentLength,k);
+ double rdd1 = AlignTools.rms_dk_diag(dist1, dist2, i, j,
+ fragmentLength, k);
double rdd2 = 0;
- if ( k2 > 0)
- rdd2 = AlignTools.rms_dk_diag(dist3,dist4,i,j,fragmentLength,k2);
+ if (k2 > 0)
+ rdd2 = AlignTools.rms_dk_diag(dist3, dist4, i, j,
+ fragmentLength, k2);
double rdd = rdd1 + rdd2;
- distanceMatrix.set(i,j,rdd);
-
-
- if ( rdd < params.getFragmentMiniDistance()) {
- FragmentPair f = new FragmentPair(fragmentLength,i,j);
- try {
-
- Atom[] catmp2 = AlignTools.getFragment(ca2, j, fragmentLength);
- Atom center2 = AlignTools.getCenter(ca2,j,fragmentLength);
-
- f.setCenter1(center1);
- f.setCenter2(center2);
-
- SVDSuperimposer svd = new SVDSuperimposer(catmp1,catmp2);
- Matrix rotmat = svd.getRotation();
- f.setRot(rotmat);
-
- Atom aunitv = (Atom)unitvector.clone();
- Calc.rotate(aunitv,rotmat);
- f.setUnitv(aunitv);
+ distanceMatrix.set(i, j, rdd);
+
+ if (rdd < params.getFragmentMiniDistance()) {
+ FragmentPair f = new FragmentPair(fragmentLength, i, j);
+ Atom[] catmp2 = AlignTools.getFragment(ca2, j,
+ fragmentLength);
+ Atom center2 = AlignTools.getCenter(ca2, j,
+ fragmentLength);
+
+ f.setCenter1(center1);
+ f.setCenter2(center2);
+
+ Matrix4d t = SuperPositions.superpose(
+ Calc.atomsToPoints(catmp1),
+ Calc.atomsToPoints(catmp2));
+
+ Matrix rotmat = Matrices.getRotationJAMA(t);
+ f.setRot(rotmat);
+
+ Atom aunitv = (Atom) unitvector.clone();
+ Calc.rotate(aunitv, rotmat);
+ f.setUnitv(aunitv);
+
+ boolean doNotAdd = false;
+ if (params.reduceInitialFragments()) {
+ doNotAdd = FragmentJoiner.reduceFragments(
+ fragments, f, distanceMatrix);
- boolean doNotAdd = false;
- if ( params.reduceInitialFragments()) {
- doNotAdd = FragmentJoiner.reduceFragments(fragments,f, distanceMatrix);
-
- }
- if ( doNotAdd)
- continue;
-
- fragments.add(f);
-
- } catch (StructureException e){
- logger.error("Exception: ", e);
- break;
}
+ if (doNotAdd)
+ continue;
+
+ fragments.add(f);
}
}
}
notifyFragmentListeners(fragments);
- FragmentPair[] fp = fragments.toArray(new FragmentPair[fragments.size()]);
+ FragmentPair[] fp = fragments
+ .toArray(new FragmentPair[fragments.size()]);
setFragmentPairs(fp);
logger.debug(" got # fragment pairs: {}", fp.length);
@@ -512,28 +530,24 @@ public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
// step 2 combine them to possible models
FragmentJoiner joiner = new FragmentJoiner();
-
JointFragments[] frags;
- if ( params.isJoinFast()) {
+ if (params.isJoinFast()) {
// apply the quick alignment procedure.
// less quality in alignments, better for DB searches...
- frags = joiner.approach_ap3(ca1,ca2,fp,params);
+ frags = joiner.approach_ap3(ca1, ca2, fp, params);
- joiner.extendFragments(ca1,ca2,frags,params);
+ joiner.extendFragments(ca1, ca2, frags, params);
- } else if ( params.isJoinPlo()){
+ } else if (params.isJoinPlo()) {
// this approach by StrComPy (peter lackner):
- frags = joiner.frag_pairwise_compat(fp,
- params.getAngleDiff(),
- params.getFragCompat(),
- params.getMaxrefine());
+ frags = joiner.frag_pairwise_compat(fp, params.getAngleDiff(),
+ params.getFragCompat(), params.getMaxrefine());
} else {
// my first implementation
- frags = joiner.approach_ap3(
- ca1,ca2, fp, params);
+ frags = joiner.approach_ap3(ca1, ca2, fp, params);
}
notifyJointFragments(frags);
@@ -543,34 +557,32 @@ public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
logger.debug("step 3 - refine alignments");
List aas = new ArrayList();
- for ( int i = 0 ; i < frags.length;i++){
+ for (int i = 0; i < frags.length; i++) {
JointFragments f = frags[i];
AlternativeAlignment a = new AlternativeAlignment();
a.apairs_from_idxlst(f);
- a.setAltAligNumber(i+1);
+ a.setAltAligNumber(i + 1);
a.setDistanceMatrix(distanceMatrix);
try {
- if ( params.getMaxIter() > 0 ){
+ if (params.getMaxIter() > 0) {
- a.refine(params,ca1,ca2);
- }
- else {
+ a.refine(params, ca1, ca2);
+ } else {
- a.finish(params,ca1,ca2);
+ a.finish(params, ca1, ca2);
}
- } catch (StructureException e){
+ } catch (StructureException e) {
logger.error("Refinement of fragment {} failed", i, e);
}
- a.calcScores(ca1,ca2);
+ a.calcScores(ca1, ca2);
aas.add(a);
}
-
// sort the alternative alignments
Comparator comp = new AltAligComparator();
- Collections.sort(aas,comp);
+ Collections.sort(aas, comp);
Collections.reverse(aas);
alts = aas.toArray(new AlternativeAlignment[aas.size()]);
@@ -581,25 +593,27 @@ public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
a.setAltAligNumber(aanbr);
}
- logger.debug("total calculation time: {} ms.", (System.currentTimeMillis()-timeStart));
+ logger.debug("total calculation time: {} ms.",
+ (System.currentTimeMillis() - timeStart));
}
- private void notifyStartingAlignment(String name1, Atom[] ca1, String name2, Atom[] ca2){
- for (AlignmentProgressListener li : listeners){
+ private void notifyStartingAlignment(String name1, Atom[] ca1,
+ String name2, Atom[] ca2) {
+ for (AlignmentProgressListener li : listeners) {
li.startingAlignment(name1, ca1, name2, ca2);
}
}
- private void notifyFragmentListeners(List fragments){
+ private void notifyFragmentListeners(List fragments) {
- for (AlignmentProgressListener li : listeners){
+ for (AlignmentProgressListener li : listeners) {
li.calculatedFragmentPairs(fragments);
}
}
- private void notifyJointFragments(JointFragments[] fragments){
- for (AlignmentProgressListener li : listeners){
+ private void notifyJointFragments(JointFragments[] fragments) {
+ for (AlignmentProgressListener li : listeners) {
li.jointFragments(fragments);
}
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CECalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CECalculator.java
index e625c48a6f..763c53e8e2 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CECalculator.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CECalculator.java
@@ -29,12 +29,13 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
import org.biojava.nbio.structure.Group;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.StructureTools;
import org.biojava.nbio.structure.align.model.AFP;
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.align.util.AFPAlignmentDisplay;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import org.biojava.nbio.core.sequence.compound.AminoAcidCompound;
import org.biojava.nbio.core.sequence.compound.AminoAcidCompoundSet;
@@ -42,6 +43,8 @@
import java.util.ArrayList;
import java.util.List;
+import javax.vecmath.Matrix4d;
+
/** This is based on the original Combinatorial Extension (CE) source code from 2003 or 2004 (CE version 2.3),
@@ -1906,26 +1909,26 @@ public double calc_rmsd(Atom[] pro1, Atom[] pro2, int strLen, boolean storeTrans
* @return RMSD
* @throws StructureException
*/
- public double calc_rmsd(Atom[] pro1, Atom[] pro2, int strLen, boolean storeTransform) throws StructureException {
+ public double calc_rmsd(Atom[] pro1, Atom[] pro2, int strLen,
+ boolean storeTransform) throws StructureException {
Atom[] cod1 = getAtoms(pro1, strLen,false);
Atom[] cod2 = getAtoms(pro2, strLen,true);
- assert(cod1.length == cod2.length);
- SVDSuperimposer svd = new SVDSuperimposer(cod1, cod2);
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(cod1),
+ Calc.atomsToPoints(cod2));
- Matrix matrix = svd.getRotation();
- Atom shift = svd.getTranslation();
+ Matrix matrix = Matrices.getRotationJAMA(trans);
+ Atom shift = Calc.getTranslationVector(trans);
if ( storeTransform) {
- r=matrix;
+ r = matrix;
t = shift;
}
- for (Atom a : cod2){
- Calc.rotate(a.getGroup(), matrix);
- Calc.shift(a.getGroup(), shift);
- }
- return SVDSuperimposer.getRMS(cod1, cod2);
+ for (Atom a : cod2)
+ Calc.transform(a.getGroup(), trans);
+
+ return Calc.rmsd(cod1, cod2);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCPMain.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCPMain.java
index 0e2f0ceaf8..79cf8d6a2e 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCPMain.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCPMain.java
@@ -29,6 +29,8 @@
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.align.util.AFPAlignmentDisplay;
import org.biojava.nbio.structure.align.util.ConfigurationException;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import java.lang.reflect.InvocationTargetException;
@@ -36,6 +38,8 @@
import java.util.Arrays;
import java.util.List;
+import javax.vecmath.Matrix4d;
+
/**
* A wrapper for {@link CeMain} which sets default parameters to be appropriate for finding
* circular permutations.
@@ -491,20 +495,20 @@ public static AFPChain filterDuplicateAFPs(AFPChain afpChain, CECalculator ceCal
Atom[] blockShifts = new Atom[blocks.size()];
if(alignLen>0) {
+
// superimpose
- SVDSuperimposer svd = new SVDSuperimposer(atoms1, atoms2);
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(atoms1),
+ Calc.atomsToPoints(atoms2));
- Matrix matrix = svd.getRotation();
- Atom shift = svd.getTranslation();
+ Matrix matrix = Matrices.getRotationJAMA(trans);
+ Atom shift = Calc.getTranslationVector(trans);
- for( Atom a : atoms2 ) {
- Calc.rotate(a.getGroup(), matrix);
- Calc.shift(a, shift);
- }
+ for( Atom a : atoms2 )
+ Calc.transform(a.getGroup(), trans);
//and get overall rmsd
- rmsd = SVDSuperimposer.getRMS(atoms1, atoms2);
- tmScore = SVDSuperimposer.getTMScore(atoms1, atoms2, ca1.length, ca2len);
+ rmsd = Calc.rmsd(atoms1, atoms2);
+ tmScore = Calc.getTMScore(atoms1, atoms2, ca1.length, ca2len);
// set all block rotations to the overall rotation
// It's not well documented if this is the expected behavior, but
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCalculatorEnhanced.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCalculatorEnhanced.java
index d68a35615f..bea92975ee 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCalculatorEnhanced.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/ce/CeCalculatorEnhanced.java
@@ -30,6 +30,8 @@
import org.biojava.nbio.structure.align.model.AFP;
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.align.util.AFPAlignmentDisplay;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import org.biojava.nbio.core.sequence.compound.AminoAcidCompound;
import org.biojava.nbio.core.sequence.compound.AminoAcidCompoundSet;
@@ -37,6 +39,8 @@
import java.util.ArrayList;
import java.util.List;
+import javax.vecmath.Matrix4d;
+
/** This is based on the original Combinatorial Extension (CE) source code from 2003 or 2004 (CE version 2.3),
@@ -1965,28 +1969,25 @@ public double calc_rmsd(Atom[] pro1, Atom[] pro2, int strLen, boolean storeTrans
Atom[] cod1 = getAtoms(pro1, strLen,false);
Atom[] cod2 = getAtoms(pro2, strLen,true);
- assert(cod1.length == cod2.length);
- SVDSuperimposer svd = new SVDSuperimposer(cod1, cod2);
-
- Matrix matrix = svd.getRotation();
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(cod1),
+ Calc.atomsToPoints(cod2));
+ Matrix matrix = Matrices.getRotationJAMA(trans);
+ Atom shift = Calc.getTranslationVector(trans);
+
if ( debug){
matrix.print(3,3);
}
- Atom shift = svd.getTranslation();
-
if ( storeTransform) {
- r=matrix;
+ r = matrix;
t = shift;
}
- for (Atom a : cod2){
-
- Calc.rotate(a.getGroup(), matrix);
- Calc.shift(a.getGroup(), shift);
- //Calc.rotate(a,r);
- //Calc.shift(a,t);
- }
+ for (Atom a : cod2)
+ Calc.transform(a.getGroup(), trans);
+
+ //Calc.rotate(a,r);
+ //Calc.shift(a,t);
// if ( show){
// StructureAlignmentJmol jmol = new StructureAlignmentJmol();
//
@@ -1999,7 +2000,7 @@ public double calc_rmsd(Atom[] pro1, Atom[] pro2, int strLen, boolean storeTrans
// jmol2.evalString("select * ; wireframe off; spacefill off; backbone on; color chain; select ligand;color cpk; wireframe 40;spacefill 120; ");
// jmol2.setTitle("calCaRmsd - pdb2");
// }
- return SVDSuperimposer.getRMS(cod1, cod2);
+ return Calc.rmsd(cod1, cod2);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPCalculator.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPCalculator.java
index 51c51f8c17..c9c6cbe7aa 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPCalculator.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPCalculator.java
@@ -29,6 +29,7 @@
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.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import java.util.ArrayList;
@@ -180,39 +181,26 @@ private static final boolean filterTerminal(Atom[] ca1, Atom[] ca2, int p1b, in
}
- private static final double getRmsd(Atom[] ca1, Atom[] ca2, int fragLen, int p1, int p2, Matrix m, Atom t) throws StructureException {
+ private static final double getRmsd(Atom[] ca1, Atom[] ca2, int fragLen,
+ int p1, int p2, Matrix m, Atom t) throws StructureException {
double rmsd = 99.9;
- Atom[] catmp1 = getFragment(ca1, p1, fragLen,false);
- Atom[] catmp2 = getFragment(ca2, p2, fragLen,true); // clone the atoms for fragment 2 so we can manipulate them...
+ Atom[] catmp1 = getFragment(ca1, p1, fragLen,false);
+ Atom[] catmp2 = getFragment(ca2, p2, fragLen,false);
- if ( catmp1 == null) {
- System.err.println("could not get fragment for ca1 " + p1 + " " + fragLen );
- return rmsd;
-
- }
-
- if ( catmp2 == null) {
- System.err.println("could not get fragment for ca2 " + p2 + " " + fragLen );
- return rmsd;
-
- }
-
- SVDSuperimposer svd = new SVDSuperimposer(catmp1, catmp2);
-
- m = svd.getRotation();
- t = svd.getTranslation();
-
- for (Atom a : catmp2){
- Calc.rotate(a,m);
- Calc.shift(a,t);
-
- }
+ if ( catmp1 == null) {
+ System.err.println("could not get fragment for ca1 " + p1 + " " + fragLen );
+ return rmsd;
+ }
- rmsd = SVDSuperimposer.getRMS(catmp1,catmp2);
+ if ( catmp2 == null) {
+ System.err.println("could not get fragment for ca2 " + p2 + " " + fragLen );
+ return rmsd;
+ }
- return rmsd;
+ return SuperPositions.getRmsd(Calc.atomsToPoints(catmp1),
+ Calc.atomsToPoints(catmp2));
}
/** get a continue subset of Atoms based by the starting position and the length
@@ -223,7 +211,8 @@ private static final double getRmsd(Atom[] ca1, Atom[] ca2, int fragLen, int p1,
* @param clone: returns a copy of the atom (in case the coordinate get manipulated...)
* @return an Atom[] array
*/
- private static final Atom[] getFragment(Atom[] caall, int pos, int fragmentLength , boolean clone){
+ private static final Atom[] getFragment(Atom[] caall, int pos, int fragmentLength ,
+ boolean clone){
if ( pos+fragmentLength > caall.length)
return null;
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPChainer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPChainer.java
index a82b9ac67d..45775748c5 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPChainer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/AFPChainer.java
@@ -28,15 +28,17 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.AFPTwister;
import org.biojava.nbio.structure.align.model.AFP;
import org.biojava.nbio.structure.align.model.AFPChain;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import java.util.List;
+import javax.vecmath.Matrix4d;
+
/** a class to chain AFPs to an alignment
*
* @author Andreas Prlic
@@ -698,18 +700,10 @@ private static double getRmsd(int focusResn, int[] focusRes1, int[] focusRes2, A
*/
private static double getRmsd(Atom[] catmp1, Atom[] catmp2) throws StructureException{
- SVDSuperimposer svd = new SVDSuperimposer(catmp1, catmp2);
-
- Matrix m = svd.getRotation();
- Atom t = svd.getTranslation();
-
- for (Atom a : catmp2){
- Calc.rotate(a,m);
- Calc.shift(a,t);
-
- }
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(catmp1),
+ Calc.atomsToPoints(catmp2));
- double rmsd = SVDSuperimposer.getRMS(catmp1,catmp2);
+ Calc.transform(catmp2, trans);
// if ( showAlig) {
// StructureAlignmentJmol jmol = new StructureAlignmentJmol();
@@ -746,9 +740,7 @@ private static double getRmsd(Atom[] catmp1, Atom[] catmp2) throws StructureExce
// jmol.evalString("model 0;");
// }
-
- return rmsd;
-
+ return Calc.rmsd(catmp1,catmp2);
}
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/StructureAlignmentOptimizer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/StructureAlignmentOptimizer.java
index 2908f91e42..dc40b6e6a5 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/StructureAlignmentOptimizer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/fatcat/calc/StructureAlignmentOptimizer.java
@@ -27,8 +27,10 @@
package org.biojava.nbio.structure.align.fatcat.calc;
+import javax.vecmath.Matrix4d;
+
import org.biojava.nbio.structure.*;
-import org.biojava.nbio.structure.jama.Matrix;
+import org.biojava.nbio.structure.geometry.SuperPositions;
@@ -257,33 +259,19 @@ private void superimposeBySet ()
}
//superimpose the equivalent residues
- SVDSuperimposer svd = new SVDSuperimposer(tmp1, tmp2);
-
- Matrix m = svd.getRotation();
- Atom t = svd.getTranslation();
-
- for (Atom a: tmp2) {
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(tmp1),
+ Calc.atomsToPoints(tmp2));
- Calc.rotate(a,m);
- Calc.shift(a,t);
-
- }
+ Calc.transform(tmp2, trans);
// weird, why does it take the RMSD before the rotation?
// the rmsd is only for the subset contained in the tmp arrays.
- rmsd = SVDSuperimposer.getRMS(tmp1,tmp2);
+ rmsd = Calc.rmsd(tmp1,tmp2);
//System.err.println("rmsd after superimpose by set: " + rmsd);
//transform structure 2 according to the superimposition of the equivalent residues
- for (i = 0 ; i< cod2.length; i++) {
- Atom a = cod2[i];
- Calc.rotate(a,m);
- Calc.shift(a,t);
-
- }
-
-
+ Calc.transform(cod2, trans);
// for(i = 0; i < equLen; i ++) {
// try {
@@ -292,9 +280,7 @@ private void superimposeBySet ()
// e.printStackTrace();
// }
// }
-
-
-
+
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/CoreSuperimposer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/CoreSuperimposer.java
index bac6dd554d..0cd14e0769 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/CoreSuperimposer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/CoreSuperimposer.java
@@ -27,12 +27,13 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.StructureTools;
import org.biojava.nbio.structure.align.multiple.Block;
import org.biojava.nbio.structure.align.multiple.BlockSet;
import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
+import org.biojava.nbio.structure.geometry.SuperPositionSVD;
+import org.biojava.nbio.structure.geometry.SuperPositions;
/**
* Superimposes the core aligned residues of every structure in a
@@ -45,7 +46,7 @@
* there is only one {@link BlockSet}, and a superposition for every
* BlockSet in case there is more than one (flexible alignment).
*
- * This class uses the {@link SVDSuperimposer} algorithm.
+ * This class uses the {@link SuperPositionSVD} algorithm.
*
* @author Aleix Lafita
* @since 4.2.0
@@ -161,9 +162,9 @@ public void superimpose(MultipleAlignment alignment)
array2 = StructureTools.cloneAtomArray(array2);
//From the superimposer we obtain the rotation and translation
- SVDSuperimposer svd = new SVDSuperimposer(array1, array2);
- Calc.transform(array2, svd.getTransformation());
- transforms.add(svd.getTransformation());
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(array1),
+ Calc.atomsToPoints(array2));
+ transforms.add(trans);
}
//Set transformation of the BlockSet
bs.setTransformations(transforms);
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/MultipleAlignmentScorer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/MultipleAlignmentScorer.java
index ed8315eb66..1422348e6d 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/MultipleAlignmentScorer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/MultipleAlignmentScorer.java
@@ -26,7 +26,6 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
import org.biojava.nbio.structure.jama.Matrix;
@@ -289,7 +288,7 @@ public static double getAvgTMScore(List transformed,
ref = Arrays.copyOf(ref, nonNullLen);
aln = Arrays.copyOf(aln, nonNullLen);
}
- sumTM += SVDSuperimposer.getTMScore(ref, aln, lengths.get(r1),
+ sumTM += Calc.getTMScore(ref, aln, lengths.get(r1),
lengths.get(r2));
comparisons++;
}
@@ -368,7 +367,7 @@ public static double getRefTMScore(List transformed,
ref = Arrays.copyOf(ref, nonNullLen);
aln = Arrays.copyOf(aln, nonNullLen);
}
- sumTM += SVDSuperimposer.getTMScore(ref, aln,
+ sumTM += Calc.getTMScore(ref, aln,
lengths.get(reference), lengths.get(r));
comparisons++;
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/ReferenceSuperimposer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/ReferenceSuperimposer.java
index 858d18279a..ca4cff2cdf 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/ReferenceSuperimposer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/ReferenceSuperimposer.java
@@ -27,12 +27,13 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.StructureTools;
import org.biojava.nbio.structure.align.multiple.Block;
import org.biojava.nbio.structure.align.multiple.BlockSet;
import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
+import org.biojava.nbio.structure.geometry.SuperPositionSVD;
+import org.biojava.nbio.structure.geometry.SuperPositions;
/**
* Superimposes each structure in a {@link MultipleAlignment} onto a reference
@@ -42,7 +43,7 @@
* there is only one {@link BlockSet}, and a superposition for every BlockSet
* in case there is more than one (flexible alignment).
*
- * This class uses the {@link SVDSuperimposer} algorithm.
+ * This class uses the {@link SuperPositionSVD} algorithm.
*
* @author Spencer Bliven
* @author Aleix Lafita
@@ -150,9 +151,9 @@ public void superimpose(MultipleAlignment alignment)
array2 = StructureTools.cloneAtomArray(array2);
//From the superimposer we obtain the rotation and translation
- SVDSuperimposer svd = new SVDSuperimposer(array1, array2);
- Calc.transform(array2, svd.getTransformation());
- transforms.add(svd.getTransformation());
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(array1),
+ Calc.atomsToPoints(array2));
+ transforms.add(trans);
}
//Set transformation of the BlockSet
bs.setTransformations(transforms);
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/AlternativeAlignment.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/AlternativeAlignment.java
index f049aecd6b..62c5015e8b 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/AlternativeAlignment.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/AlternativeAlignment.java
@@ -28,6 +28,8 @@
import org.biojava.nbio.structure.align.helper.AligMatEl;
import org.biojava.nbio.structure.align.helper.IndexPair;
import org.biojava.nbio.structure.align.helper.JointFragments;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import java.io.Serializable;
@@ -35,6 +37,8 @@
import java.util.List;
import java.util.logging.Logger;
+import javax.vecmath.Matrix4d;
+
/**
* Implements a class which handles one possible (alternative) solution.
*
@@ -776,13 +780,14 @@ private void super_pos_alig(Atom[]ca1,Atom[]ca2,int[] idx1, int[] idx2, boolean
ca2subset[i] = (Atom) ca2[pos2].clone();
}
- SVDSuperimposer svd = new SVDSuperimposer(ca1subset,ca2subset);
- this.currentRotMatrix = svd.getRotation();
- this.currentTranMatrix = svd.getTranslation();
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1subset),
+ Calc.atomsToPoints(ca2subset));
+ this.currentRotMatrix = Matrices.getRotationJAMA(trans);
+ this.currentTranMatrix = Calc.getTranslationVector(trans);
//currentRotMatrix.print(3,3);
if ( getRMS) {
rotateShiftAtoms(ca2subset);
- this.rms = SVDSuperimposer.getRMS(ca1subset,ca2subset);
+ this.rms = Calc.rmsd(ca1subset,ca2subset);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/FragmentJoiner.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/FragmentJoiner.java
index 3a6bc141fa..47cac82c94 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/FragmentJoiner.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/pairwise/FragmentJoiner.java
@@ -24,7 +24,6 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.StrucAligParameters;
import org.biojava.nbio.structure.align.helper.AlignTools;
@@ -314,7 +313,7 @@ public static double getRMS(Atom[] ca1, Atom[]ca2,JointFragments frag) throws St
Calc.shift(a, atom);
}
- rms = SVDSuperimposer.getRMS(ca1subset,ca2subset);
+ rms = Calc.rmsd(ca1subset,ca2subset);
return rms;
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/quaternary/QsAlign.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/quaternary/QsAlign.java
index 928d343c72..83cd550e61 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/quaternary/QsAlign.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/quaternary/QsAlign.java
@@ -11,7 +11,6 @@
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.Structure;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.StructureTools;
@@ -30,6 +29,7 @@
import org.biojava.nbio.structure.cluster.SubunitClustererParameters;
import org.biojava.nbio.structure.cluster.SubunitExtractor;
import org.biojava.nbio.structure.contact.Pair;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.geometry.UnitQuaternions;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -180,7 +180,7 @@ public static QsAlignResult align(List s1, List s2,
}
// 3- Check the RMSD condition
- double rmsd = SVDSuperimposer.getRMS(atoms1, atoms2c);
+ double rmsd = Calc.rmsd(atoms1, atoms2c);
if (rmsd > aParams.getMaxRmsd()) {
logger.debug(String.format("Subunit matching %d "
+ "vs %d of cluster %d could not be "
@@ -354,9 +354,8 @@ private static Matrix4d getTransformForClusterSubunitMap(
Pair pair = getAlignedAtomsForClusterSubunitMap(clusters,
clusterSubunitMap);
- SVDSuperimposer svd = new SVDSuperimposer(pair.getFirst(),
- pair.getSecond());
- return svd.getTransformation();
+ return SuperPositions.superpose(Calc.atomsToPoints(pair.getFirst()),
+ Calc.atomsToPoints(pair.getSecond()));
}
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPAlignmentDisplay.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPAlignmentDisplay.java
index 286aa3d7f5..21298c906e 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPAlignmentDisplay.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPAlignmentDisplay.java
@@ -22,6 +22,8 @@
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.ce.GuiWrapper;
import org.biojava.nbio.structure.align.model.AFPChain;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -29,6 +31,8 @@
import java.lang.reflect.InvocationTargetException;
import java.util.*;
+import javax.vecmath.Matrix4d;
+
public class AFPAlignmentDisplay
{
@@ -65,9 +69,10 @@ public static Matrix getRotMax(AFPChain afpChain,Atom[] ca1,Atom[] ca2) throws S
Atom[] a1 = getAlignedAtoms1(afpChain,ca1);
Atom[] a2 = getAlignedAtoms2(afpChain,ca2);
- SVDSuperimposer svd = new SVDSuperimposer(a1,a2);
-
- return svd.getRotation();
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(a1),
+ Calc.atomsToPoints(a2));
+
+ return Matrices.getRotationJAMA(trans);
}
@@ -77,9 +82,10 @@ public static Atom getTranslation(AFPChain afpChain,Atom[] ca1,Atom[] ca2) throw
Atom[] a1 = getAlignedAtoms1(afpChain,ca1);
Atom[] a2 = getAlignedAtoms2(afpChain,ca2);
- SVDSuperimposer svd = new SVDSuperimposer(a1,a2);
-
- return svd.getTranslation();
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(a1),
+ Calc.atomsToPoints(a2));
+
+ return Calc.getTranslationVector(trans);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPChainScorer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPChainScorer.java
index c9491196b0..030510157a 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPChainScorer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AFPChainScorer.java
@@ -25,12 +25,13 @@
package org.biojava.nbio.structure.align.util;
+import javax.vecmath.Matrix4d;
+
import org.biojava.nbio.structure.Atom;
import org.biojava.nbio.structure.Calc;
-import org.biojava.nbio.structure.SVDSuperimposer;
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.model.AFPChain;
-import org.biojava.nbio.structure.jama.Matrix;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -79,16 +80,12 @@ public static double getTMScore(AFPChain align, Atom[] ca1, Atom[] ca2) throws
ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
}
//Superimpose
- SVDSuperimposer svd = new SVDSuperimposer(ca1aligned, ca2aligned);
- Matrix matrix = svd.getRotation();
- Atom shift = svd.getTranslation();
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
+ Calc.atomsToPoints(ca2aligned));
- for(Atom a : ca2aligned) {
- Calc.rotate(a, matrix);
- Calc.shift(a, shift);
- }
+ Calc.transform(ca2aligned, trans);
- return SVDSuperimposer.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
+ return Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java
index 65d77182d3..3912226463 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java
@@ -27,6 +27,9 @@
import org.biojava.nbio.structure.align.fatcat.FatCatRigid;
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.align.xml.AFPChainXMLParser;
+import org.biojava.nbio.structure.geometry.Matrices;
+import org.biojava.nbio.structure.geometry.SuperPositionSVD;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -38,6 +41,8 @@
import java.util.regex.Matcher;
import java.util.regex.Pattern;
+import javax.vecmath.Matrix4d;
+
/**
* Methods for analyzing and manipulating AFPChains and for
* other pairwise alignment utilities.
@@ -899,9 +904,13 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
}
//Superimpose the two structures in correspondance to the new alignment
- SVDSuperimposer svd = new SVDSuperimposer(ca1aligned, ca2aligned);
- Matrix matrix = svd.getRotation();
- Atom shift = svd.getTranslation();
+ Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
+ Calc.atomsToPoints(ca2aligned));
+
+ Matrix matrix = Matrices.getRotationJAMA(trans);
+ Atom shift = Calc.getTranslationVector(trans);
+
+
Matrix[] blockMxs = new Matrix[afpChain.getBlockNum()];
Arrays.fill(blockMxs, matrix);
afpChain.setBlockRotationMatrix(blockMxs);
@@ -915,8 +924,8 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
}
//Calculate the RMSD and TM score for the new alignment
- double rmsd = SVDSuperimposer.getRMS(ca1aligned, ca2aligned);
- double tmScore = SVDSuperimposer.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
+ double rmsd = Calc.rmsd(ca1aligned, ca2aligned);
+ double tmScore = Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
afpChain.setTotalRmsdOpt(rmsd);
afpChain.setTMScore(tmScore);
@@ -944,16 +953,17 @@ public static void updateSuperposition(AFPChain afpChain, Atom[] ca1, Atom[] ca2
ca2block = (Atom[]) resizeArray(ca2block, position);
}
//Superimpose the two block structures
- SVDSuperimposer svdb = new SVDSuperimposer(ca1block, ca2block);
- Matrix matrixb = svdb.getRotation();
- Atom shiftb = svdb.getTranslation();
- for (Atom a : ca2block) {
- Calc.rotate(a, matrixb);
- Calc.shift(a, shiftb);
- }
+ Matrix4d transb = SuperPositions.superpose(Calc.atomsToPoints(ca1block),
+ Calc.atomsToPoints(ca2block));
+
+ Matrix matrixb = Matrices.getRotationJAMA(trans);
+ Atom shiftb = Calc.getTranslationVector(trans);
+
+ Calc.transform(ca2block, transb);
+
//Calculate the RMSD and TM score for the block
- double rmsdb = SVDSuperimposer.getRMS(ca1block, ca2block);
- double tmScoreb = SVDSuperimposer.getTMScore(ca1block, ca2block, ca1.length, ca2.length);
+ double rmsdb = Calc.rmsd(ca1block, ca2block);
+ double tmScoreb = Calc.getTMScore(ca1block, ca2block, ca1.length, ca2.length);
blockRMSD[k] = rmsdb;
blockScore[k] = tmScoreb;
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AtomCache.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AtomCache.java
index 26d717e849..ecb0c8f665 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AtomCache.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AtomCache.java
@@ -72,7 +72,7 @@ public class AtomCache {
* The default output bioassembly style: if true the bioassemblies are multimodel,
* if false the bioassemblies are flat with renamed chains for symmetry-partners.
*/
- public static final boolean DEFAULT_BIOASSEMBLY_STYLE = true;
+ public static final boolean DEFAULT_BIOASSEMBLY_STYLE = false;
public static final String BIOL_ASSEMBLY_IDENTIFIER = "BIO:";
public static final String CHAIN_NR_SYMBOL = ":";
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/RotationAxis.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/RotationAxis.java
index 779e78be6c..91d87512e1 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/RotationAxis.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/RotationAxis.java
@@ -26,6 +26,7 @@
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.contact.Pair;
+import org.biojava.nbio.structure.geometry.Matrices;
import org.biojava.nbio.structure.jama.Matrix;
import javax.vecmath.AxisAngle4d;
@@ -175,7 +176,7 @@ public RotationAxis(Matrix rotation, Atom translation) {
*/
public RotationAxis(Matrix4d transform) {
- Matrix rot = Calc.getRotationMatrix(transform);
+ Matrix rot = Matrices.getRotationJAMA(transform);
Atom transl = Calc.getTranslationVector(transform);
init(rot,transl);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitExtractor.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitExtractor.java
index d5e943a6c8..bb628f37aa 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitExtractor.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitExtractor.java
@@ -68,31 +68,24 @@ public static List extractSubunits(Structure structure,
// The extracted subunit container
List subunits = new ArrayList();
- // Biological assemblies require all models
- int models = 1;
- if (structure.isBiologicalAssembly())
- models = structure.nrModels();
+
+ for (Chain c : structure.getChains()) {
- logger.info("Protein models used in calculation: " + models);
+ // Only take protein chains
+ if (c.isProtein()) {
+ Atom[] ca = StructureTools.getRepresentativeAtomArray(c);
+ logger.debug("Chain " + c.getId() + "; CA Atoms: "
+ + ca.length + "; SEQRES: " + c.getSeqResSequence());
+ subunits.add(new Subunit(ca, c.getName(), null, structure));
- for (int i = 0; i < models; i++) {
- for (Chain c : structure.getChains(i)) {
-
- // Only take protein chains
- if (StructureTools.isProtein(c)) {
- Atom[] ca = StructureTools.getRepresentativeAtomArray(c);
- logger.info("Chain " + c.getId() + "; CA Atoms: "
- + ca.length + "; SEQRES: " + c.getSeqResSequence());
- subunits.add(new Subunit(ca, c.getName(), null, structure));
-
- }
}
}
+
// Calculate the minimum length of a Subunit
int adjustedMinLen = calcAdjustedMinimumSequenceLength(subunits,
absMinLen, fraction, minLen);
- logger.info("Adjusted minimum sequence length: " + adjustedMinLen);
+ logger.debug("Adjusted minimum sequence length: " + adjustedMinLen);
// Filter out short Subunits
for (int s = subunits.size() - 1; s >= 0; s--) {
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/CalcPoint.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/CalcPoint.java
index dd282dcb30..bb5a63c50e 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/CalcPoint.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/CalcPoint.java
@@ -93,8 +93,8 @@ public static Point3d[] clonePoint3dArray(Point3d[] x) {
/**
* Peter can you document this method? TODO
*
- * @param a
- * @param b
+ * @param moved
+ * @param fixed
* @return
*/
public static Matrix formMatrix(Point3d[] a, Point3d[] b) {
@@ -135,4 +135,155 @@ public static Matrix formMatrix(Point3d[] a, Point3d[] b) {
return new Matrix(f);
}
+ /**
+ * Returns the TM-Score for two superimposed sets of coordinates Yang Zhang
+ * and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics
+ * 57:702–710 (2004)
+ *
+ * @param x
+ * coordinate set 1
+ * @param y
+ * coordinate set 2
+ * @param lengthNative
+ * total length of native sequence
+ * @return
+ */
+ public static double TMScore(Point3d[] x, Point3d[] y, int lengthNative) {
+
+ if (x.length != y.length) {
+ throw new IllegalArgumentException(
+ "Point arrays are not of the same length.");
+ }
+
+ double d0 = 1.24 * Math.cbrt(x.length - 15.0) - 1.8;
+ double d0Sq = d0 * d0;
+
+ double sum = 0;
+ for (int i = 0; i < x.length; i++) {
+ sum += 1.0 / (1.0 + x[i].distanceSquared(y[i]) / d0Sq);
+ }
+
+ return sum / lengthNative;
+ }
+
+ /*
+ * Needs documentation!
+ *
+ * @param x
+ *
+ * @param y
+ *
+ * @return
+ */
+ public static double GTSlikeScore(Point3d[] x, Point3d[] y) {
+
+ if (x.length != y.length) {
+ throw new IllegalArgumentException(
+ "Point arrays are not of the same length.");
+ }
+
+ int contacts = 0;
+
+ for (Point3d px : x) {
+ double minDist = Double.MAX_VALUE;
+
+ for (Point3d py : y) {
+ minDist = Math.min(minDist, px.distanceSquared(py));
+ }
+
+ if (minDist > 64)
+ continue;
+ contacts++;
+
+ if (minDist > 16)
+ continue;
+ contacts++;
+
+ if (minDist > 4)
+ continue;
+ contacts++;
+
+ if (minDist > 1)
+ continue;
+ contacts++;
+ }
+
+ return contacts * 25.0 / x.length;
+ }
+
+ /**
+ * Calculate the RMSD of two point arrays, already superposed.
+ *
+ * @param x
+ * array of points superposed to y
+ * @param y
+ * array of points superposed to x
+ * @return RMSD
+ */
+ public static double rmsd(Point3d[] x, Point3d[] y) {
+
+ if (x.length != y.length) {
+ throw new IllegalArgumentException(
+ "Point arrays are not of the same length.");
+ }
+
+ double sum = 0.0;
+ for (int i = 0; i < x.length; i++) {
+ sum += x[i].distanceSquared(y[i]);
+ }
+ return Math.sqrt(sum / x.length);
+ }
+
+ /*
+ * Needs documentation!
+ *
+ * @param x
+ *
+ * @param y
+ *
+ * @return
+ */
+ public static double rmsdMin(Point3d[] x, Point3d[] y) {
+
+ if (x.length != y.length) {
+ throw new IllegalArgumentException(
+ "Point arrays are not of the same length.");
+ }
+
+ double sum = 0.0;
+ for (int i = 0; i < x.length; i++) {
+ double minDist = Double.MAX_VALUE;
+ for (int j = 0; j < y.length; j++) {
+ minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
+ }
+ sum += minDist;
+ }
+ return Math.sqrt(sum / x.length);
+ }
+
+ /*
+ * Needs documentation!
+ *
+ * @param x
+ *
+ * @param y
+ *
+ * @param maxDistance
+ *
+ * @return
+ */
+ public static int contacts(Point3d[] x, Point3d[] y, double maxDistance) {
+ int contacts = 0;
+ for (int i = 0; i < x.length; i++) {
+ double minDist = Double.MAX_VALUE;
+ for (int j = 0; j < y.length; j++) {
+ minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
+ }
+ if (minDist < maxDistance * maxDistance) {
+ contacts++;
+ }
+ }
+ return contacts;
+ }
+
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/Matrices.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/Matrices.java
new file mode 100644
index 0000000000..222d8cf840
--- /dev/null
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/Matrices.java
@@ -0,0 +1,89 @@
+package org.biojava.nbio.structure.geometry;
+
+import javax.vecmath.Matrix3d;
+import javax.vecmath.Matrix4d;
+import javax.vecmath.Vector3d;
+
+import org.biojava.nbio.structure.jama.Matrix;
+
+/**
+ * Matrices contains static methods to operate and transform matrices used in 3D
+ * geometry (transformation matrices and rotation matrices).
+ *
+ * This class complements and extends the functionallity of vecmath and JAMA.
+ *
+ * @author Aleix Lafita
+ * @since 5.0.0
+ *
+ */
+public class Matrices {
+
+ /** Prevent instantiation */
+ private Matrices(){}
+
+ /**
+ * Convert a transformation matrix into a JAMA rotation matrix. Because the
+ * JAMA matrix is a pre-multiplication matrix and the Vecmath matrix is a
+ * post-multiplication one, the rotation matrix is transposed to ensure that
+ * the transformation they produce is the same.
+ *
+ * @param transform
+ * Matrix4d with transposed rotation matrix
+ * @return rotation matrix as JAMA object
+ */
+ public static Matrix getRotationJAMA(Matrix4d transform) {
+
+ Matrix rot = new Matrix(3, 3);
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ rot.set(j, i, transform.getElement(i, j)); // transposed
+ }
+ }
+ return rot;
+ }
+
+ /**
+ * Convert a transformation matrix into a rotation matrix.
+ *
+ * @param transform
+ * Matrix4d
+ * @return rotation matrix
+ */
+ public static Matrix3d getRotationMatrix(Matrix4d transform) {
+
+ Matrix3d rot = new Matrix3d();
+ transform.setRotationScale(rot);
+ return rot;
+ }
+
+ /**
+ * Extract the translational vector of a transformation matrix.
+ *
+ * @param transform
+ * Matrix4d
+ * @return Vector3d translation vector
+ */
+ public static Vector3d getTranslationVector(Matrix4d transform) {
+ Vector3d transl = new Vector3d();
+ transform.get(transl);
+ return transl;
+ }
+
+ /**
+ * Convert JAMA rotation and translation to a Vecmath transformation matrix.
+ * Because the JAMA matrix is a pre-multiplication matrix and the Vecmath
+ * matrix is a post-multiplication one, the rotation matrix is transposed to
+ * ensure that the transformation they produce is the same.
+ *
+ * @param rot
+ * 3x3 Rotation matrix
+ * @param trans
+ * 3x1 Translation matrix
+ * @return 4x4 transformation matrix
+ */
+ public static Matrix4d getTransformation(Matrix rot, Matrix trans) {
+ return new Matrix4d(new Matrix3d(rot.getColumnPackedCopy()),
+ new Vector3d(trans.getColumnPackedCopy()), 1.0);
+ }
+
+}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPosition.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPosition.java
index f7abe7deb0..b5844da94d 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPosition.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPosition.java
@@ -1,305 +1,79 @@
-/*
- * BioJava development 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.
- *
- * For more information on the BioJava project and its aims,
- * or to join the biojava-l mailing list, visit the home page
- * at:
- *
- * http://www.biojava.org/
- *
- */
package org.biojava.nbio.structure.geometry;
-import javax.vecmath.*;
+import javax.vecmath.Matrix4d;
+import javax.vecmath.Point3d;
/**
- * The SuperPosition is a static Class that contains methods to superpose arrays
- * of Points in 3D.
- *
- * @author Peter Rose
+ * The SuperPosition interface defines and documents the required methods for
+ * any superpostion algorithm implementation, so that the input and expected
+ * output are uniform.
*
+ * @author Aleix Lafita
+ * @since 5.0.0
+ *
*/
-public final class SuperPosition {
-
- /** Prevent instantiation */
- private SuperPosition() {
- }
-
- /**
- * This method superposes x onto y, so it transforms the coordinates of x.
- * The translation component of the matrix is wrong.
- *
- * @param x
- * point array, point coordinates will be modified
- * @param y
- * point array, point coordinates will not be modified
- * @return the Matrix4d used for superposition, without translation
- * @deprecated use {@link #superposeWithTranslation(Point3d[], Point3d[])}
- */
- @Deprecated
- public static Matrix4d superpose(Point3d[] x, Point3d[] y) {
- // superpose x onto y
- Point3d[] ref = CalcPoint.clonePoint3dArray(y);
-
- Point3d ytrans = CalcPoint.centroid(ref);
- ytrans.negate();
- CalcPoint.translate(ytrans, ref);
-
- CalcPoint.center(x);
-
- // calculate quaternion from relative orientation
- Quat4d q = UnitQuaternions.relativeOrientation(x, ref);
-
- Matrix4d rotTrans = new Matrix4d();
- rotTrans.set(q);
-
- // set translational component of transformation matrix TODO wrong
- ytrans.negate();
- rotTrans.setTranslation(new Vector3d(ytrans));
-
- // tranform coordinates
- CalcPoint.transform(rotTrans, x);
-
- return rotTrans;
- }
+public interface SuperPosition {
/**
- * This method superposes x onto y, so it transforms the coordinates of x.
+ * Obtain the superposition matrix that minimizes the RMSD between two
+ * arrays of equivalent points.
+ *
+ * The two point arrays have to be of the same length and the order of
+ * points have to be the same, so that a specific position in the one array
+ * is equivalent to the same position in the other array.
*
- * @param x
- * point array, point coordinates will be modified
- * @param y
- * point array, point coordinates will not be modified
- * @return the Matrix4d used for superposition
+ * @param fixed
+ * point array as reference, onto which the other point array is
+ * superposed. Original coordinates will not be modified.
+ * @param moved
+ * point array to which the resulting transformation matrix is
+ * applied. Original coordinates will not be modified.
+ * @return transformation matrix as a Matrix4d to superpose moved onto fixed
+ * point arrays
*/
- public static Matrix4d superposeWithTranslation(Point3d[] x, Point3d[] y) {
-
- // translate to origin
- Point3d[] xref = CalcPoint.clonePoint3dArray(x);
- Point3d xtrans = CalcPoint.centroid(xref);
- xtrans.negate();
- CalcPoint.translate(xtrans, xref);
-
- Point3d[] yref = CalcPoint.clonePoint3dArray(y);
- Point3d ytrans = CalcPoint.centroid(yref);
- ytrans.negate();
- CalcPoint.translate(ytrans, yref);
-
- // calculate rotational component (rotation around origin)
- Quat4d q = UnitQuaternions.relativeOrientation(xref, yref);
- Matrix4d rotTrans = new Matrix4d();
- rotTrans.set(q);
-
- // combine with x -> origin translation
- Matrix4d trans = new Matrix4d();
- trans.setIdentity();
- trans.setTranslation(new Vector3d(xtrans));
- rotTrans.mul(rotTrans, trans);
-
- // combine with origin -> y translation
- ytrans.negate();
- Matrix4d transInverse = new Matrix4d();
- transInverse.setIdentity();
- transInverse.setTranslation(new Vector3d(ytrans));
- rotTrans.mul(transInverse, rotTrans);
-
- // transform x coordinates onto y coordinate frame
- CalcPoint.transform(rotTrans, x);
-
- return rotTrans;
- }
+ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved);
/**
- * This method superposes x onto y, so it transforms the coordinates of x.
- * It requires the points to be centered (use
- * {@link CalcPoint#center(Point3d[])}, because the transformation matrix
- * will not have translation component (faster).
+ * Transform an array of points so that the coordinates of its points
+ * minimize the RMSD to the other array of equivalent points, and return the
+ * transformation matrix applied.
+ *
+ * The two point arrays have to be of the same length and the order of
+ * points have to be the same, so that a specific position in the one array
+ * is equivalent to the same position in the other array.
*
- * @param x
- * point array, point coordinates will be modified
- * @param y
- * point array, point coordinates will not be modified
- * @return the Matrix4d used for superposition, without translation
+ * @param fixed
+ * point array as reference, onto which the other point array is
+ * superposed. Original coordinates will not be modified.
+ * @param moved
+ * point array to which the resulting transformation matrix is
+ * applied. Original coordinates will be transformed.
+ * @return transformation matrix as a Matrix4d to superpose moved onto fixed
+ * point arrays
*/
- public static Matrix4d superposeAtOrigin(Point3d[] x, Point3d[] y) {
- Quat4d q = UnitQuaternions.relativeOrientation(x, y);
-
- Matrix4d rotTrans = new Matrix4d();
- rotTrans.set(q);
-
- CalcPoint.transform(rotTrans, x);
-
- return rotTrans;
- }
+ public Matrix4d superposeAndTransform(Point3d[] fixed, Point3d[] moved);
/**
- * This method superposes x onto y, so it transforms the coordinates of x.
- * It requires the points to be centered (use
- * {@link CalcPoint#center(Point3d[])}, because the transformation matrix
- * will not have translation component (faster).
+ * Calculate the RMSD between two arrays of equivalent points that are not
+ * superposed.
+ *
+ * This is equivalent to first superposing the point arrays with
+ * {@link SuperPosition#superposeAndTransform(Point3d[], Point3d[])} and
+ * then calculating the RMSD of the superposed point arrays with
+ * {@link CalcPoint#rmsd(Point3d[], Point3d[])}, but it will be faster when
+ * the transformation matrix is not needed.
+ *
+ * The two point arrays have to be of the same length and the order of
+ * points have to be the same, so that a specific position in the one array
+ * is equivalent to the same position in the other array.
*
* @param x
- * point array, point coordinates will be modified
- * @param y
- * point array, point coordinates will not be modified
- * @param axisAngle
- * to be filled with rotation axis
- * @return the Matrix4d used for superposition, without translation
- */
- public static Matrix4d superposeAtOrigin(Point3d[] x, Point3d[] y,
- AxisAngle4d axisAngle) {
- Quat4d q = UnitQuaternions.relativeOrientation(x, y);
- Matrix4d rotTrans = new Matrix4d();
- rotTrans.setIdentity();
- rotTrans.set(q);
- axisAngle.set(q);
- Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
- if (axis.lengthSquared() < 1.0E-6) {
- // System.err.println("Error: SuperPosition.superposeAtOrigin: axis vector undefined!");
- axisAngle.x = 0;
- axisAngle.y = 0;
- axisAngle.z = 1;
- axisAngle.angle = 0;
- } else {
- axis.normalize();
- axisAngle.x = axis.x;
- axisAngle.y = axis.y;
- axisAngle.z = axis.z;
- }
- CalcPoint.transform(rotTrans, x);
-
- return rotTrans;
- }
-
- /**
- * Calculate the RMSD of two point arrays, already superposed.
- *
- * @param x
- * array of points superposed to y
- * @param y
- * array of points superposed to x
- * @return RMSD
- */
- public static double rmsd(Point3d[] x, Point3d[] y) {
- double sum = 0.0;
- for (int i = 0; i < x.length; i++) {
- sum += x[i].distanceSquared(y[i]);
- }
- return Math.sqrt(sum / x.length);
- }
-
- /*
- * Needs documentation!
- *
- * @param x
- *
- * @param y
- *
- * @return
- */
- public static double rmsdMin(Point3d[] x, Point3d[] y) {
- double sum = 0.0;
- for (int i = 0; i < x.length; i++) {
- double minDist = Double.MAX_VALUE;
- for (int j = 0; j < y.length; j++) {
- minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
- }
- sum += minDist;
- }
- return Math.sqrt(sum / x.length);
- }
-
- /**
- * Returns the TM-Score for two superimposed sets of coordinates Yang Zhang
- * and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics
- * 57:702–710 (2004)
- *
- * @param x
- * coordinate set 1
- * @param y
- * coordinate set 2
- * @param lengthNative
- * total length of native sequence
- * @return
- */
- public static double TMScore(Point3d[] x, Point3d[] y, int lengthNative) {
- double d0 = 1.24 * Math.cbrt(x.length - 15.0) - 1.8;
- double d0Sq = d0 * d0;
-
- double sum = 0;
- for (int i = 0; i < x.length; i++) {
- sum += 1.0 / (1.0 + x[i].distanceSquared(y[i]) / d0Sq);
- }
-
- return sum / lengthNative;
- }
-
- /*
- * Needs documentation!
- *
- * @param x
- * @param y
- * @return
- */
- public static double GTSlikeScore(Point3d[] x, Point3d[] y) {
- int contacts = 0;
-
- for (Point3d px : x) {
- double minDist = Double.MAX_VALUE;
-
- for (Point3d py : y) {
- minDist = Math.min(minDist, px.distanceSquared(py));
- }
-
- if (minDist > 64)
- continue;
- contacts++;
-
- if (minDist > 16)
- continue;
- contacts++;
-
- if (minDist > 4)
- continue;
- contacts++;
-
- if (minDist > 1)
- continue;
- contacts++;
- }
-
- return contacts * 25.0 / x.length;
- }
-
- /*
- * Needs documentation!
- * @param x
+ * an array of points. Original coordinates will not be modified.
* @param y
- * @param maxDistance
- * @return
+ * an array of points. Original coordinates will not be modified.
+ * @return the minimum RMSD between the equivalent point arrays (after
+ * superposition)
*/
- public static int contacts(Point3d[] x, Point3d[] y, double maxDistance) {
- int contacts = 0;
- for (int i = 0; i < x.length; i++) {
- double minDist = Double.MAX_VALUE;
- for (int j = 0; j < y.length; j++) {
- minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
- }
- if (minDist < maxDistance * maxDistance) {
- contacts++;
- }
- }
- return contacts;
- }
+ public double getRmsd(Point3d[] x, Point3d[] y);
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionAbstract.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionAbstract.java
new file mode 100644
index 0000000000..6107e6c266
--- /dev/null
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionAbstract.java
@@ -0,0 +1,50 @@
+package org.biojava.nbio.structure.geometry;
+
+import javax.vecmath.Matrix4d;
+import javax.vecmath.Point3d;
+
+/**
+ * The SuperPositionAbstract contains common code shared by all SuperPosition
+ * algorithm implementations.
+ *
+ * @author Aleix Lafita
+ * @since 5.0.0
+ *
+ */
+public abstract class SuperPositionAbstract implements SuperPosition {
+
+ protected boolean centered;
+
+ public SuperPositionAbstract(boolean centered) {
+ this.centered = centered;
+ }
+
+ @Override
+ public Matrix4d superposeAndTransform(Point3d[] fixed, Point3d[] moved) {
+ Matrix4d rotTrans = superpose(fixed, moved);
+ CalcPoint.transform(rotTrans, moved);
+ return rotTrans;
+ }
+
+ /**
+ * Check that the input to the superposition algorithms is valid.
+ *
+ * @param fixed
+ * @param moved
+ */
+ protected void checkInput(Point3d[] fixed, Point3d[] moved) {
+ if (fixed.length != moved.length)
+ throw new IllegalArgumentException(
+ "Point arrays to superpose are of different lengths.");
+ }
+
+ /**
+ * @param centered
+ * true if the point arrays are centered at the origin (faster),
+ * false otherwise
+ */
+ public void setCentered(boolean centered) {
+ this.centered = centered;
+ }
+
+}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQCP.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQCP.java
index 765762d5e4..a5d9a4aa1c 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQCP.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQCP.java
@@ -107,20 +107,13 @@
* @author Peter Rose (adopted to Java)
* @author Aleix Lafita (adopted to Java)
*/
-public final class SuperPositionQCP {
+public final class SuperPositionQCP extends SuperPositionAbstract {
private static final Logger logger = LoggerFactory
.getLogger(SuperPositionQCP.class);
- /**
- * The required eigenvector precission
- */
- private static final double EVEC_PREC = 1E-6;
-
- /**
- * The required eigenvalue precission
- */
- private static final double EVAL_PREC = 1E-11;
+ private double evec_prec = 1E-6;
+ private double eval_prec = 1E-11;
private Point3d[] x;
private Point3d[] y;
@@ -145,20 +138,30 @@ public final class SuperPositionQCP {
private boolean transformationCalculated = false;
private boolean centered = false;
- /**
- * Default constructor
- */
- public SuperPositionQCP() {
- this.centered = false;
+ /**
+ * Default constructor for the quaternion based superposition algorithm.
+ *
+ * @param centered
+ * true if the point arrays are centered at the origin (faster),
+ * false otherwise
+ */
+ public SuperPositionQCP(boolean centered) {
+ super(centered);
}
/**
- * Constructor with option to set centered flag. This constructor
- * should be used if both coordinate input set have been centered at the origin.
- * @param centered if set true, the input coordinates are already centered at the origin
+ * Constructor with option to set the precision values.
+ *
+ * @param centered
+ * true if the point arrays are centered at the origin (faster),
+ * false otherwise
+ * @param evec_prec required eigenvector precision
+ * @param eval_prec required eigenvalue precision
*/
- public SuperPositionQCP(boolean centered) {
- this.centered = centered;
+ public SuperPositionQCP(boolean centered, double evec_prec, double eval_prec) {
+ super(centered);
+ this.evec_prec = evec_prec;
+ this.eval_prec = eval_prec;
}
/**
@@ -167,7 +170,7 @@ public SuperPositionQCP(boolean centered) {
* @param x 3d points of reference coordinate set
* @param y 3d points of coordinate set for superposition
*/
- public void set(Point3d[] x, Point3d[] y) {
+ private void set(Point3d[] x, Point3d[] y) {
this.x = x;
this.y = y;
rmsdCalculated = false;
@@ -182,7 +185,7 @@ public void set(Point3d[] x, Point3d[] y) {
* @param y 3d points of coordinate set for superposition
* @param weight a weight in the inclusive range [0,1] for each point
*/
- public void set(Point3d[] x, Point3d[] y, double[] weight) {
+ private void set(Point3d[] x, Point3d[] y, double[] weight) {
this.x = x;
this.y = y;
this.weight = weight;
@@ -190,10 +193,6 @@ public void set(Point3d[] x, Point3d[] y, double[] weight) {
transformationCalculated = false;
}
- public void setCentered(boolean centered) {
- this.centered = centered;
- }
-
/**
* Return the RMSD of the superposition of input coordinate set y onto x.
* Note, this is the fasted way to calculate an RMSD without actually
@@ -201,19 +200,16 @@ public void setCentered(boolean centered) {
* calculations are only performed if necessary.
* @return root mean square deviation for superposition of y onto x
*/
- public double getRmsd() {
+ private double getRmsd() {
if (!rmsdCalculated) {
calcRmsd(x, y);
}
return rmsd;
}
- /**
- * Returns a 4x4 transformation matrix that transforms the y coordinates onto the x coordinates.
- * The calculation is performed "lazy", meaning calculations are only performed if necessary.
- * @return 4x4 transformation matrix to transform y coordinates onto x
- */
- public Matrix4d getTransformationMatrix() {
+ @Override
+ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
+ set(moved, fixed);
getRotationMatrix();
if (!centered) {
calcTransformation();
@@ -222,8 +218,27 @@ public Matrix4d getTransformationMatrix() {
}
return transformation;
}
+
+ /**
+ * Weighted superposition.
+ *
+ * @param fixed
+ * @param moved
+ * @param weight array of weigths for each equivalent point position
+ * @return
+ */
+ public Matrix4d weightedSuperpose(Point3d[] fixed, Point3d[] moved, double[] weight) {
+ set(moved, fixed, weight);
+ getRotationMatrix();
+ if (!centered) {
+ calcTransformation();
+ } else {
+ transformation.set(rotmat);
+ }
+ return transformation;
+ }
- public Matrix3d getRotationMatrix() {
+ private Matrix3d getRotationMatrix() {
getRmsd();
if (!transformationCalculated) {
calcRotationMatrix();
@@ -231,16 +246,6 @@ public Matrix3d getRotationMatrix() {
return rotmat;
}
- /**
- * Returns the transformed (superposed) y coordinates
- * TODO this is actually returning a transformation of the x coordinates, we should review if that's correct and change it if not
- * @return transformed y coordinates
- */
- public Point3d[] getTransformedCoordinates() {
- CalcPoint.transform(transformation, x);
- return x;
- }
-
/**
* Calculates the RMSD value for superposition of y onto x.
* This requires the coordinates to be precentered.
@@ -433,7 +438,7 @@ private int calcRmsd(double len) {
double delta = ((a * mxEigenV + c0) / (2.0 * x2 * mxEigenV + b + a));
mxEigenV -= delta;
- if (Math.abs(mxEigenV - oldg) < Math.abs(EVAL_PREC*mxEigenV))
+ if (Math.abs(mxEigenV - oldg) < Math.abs(eval_prec*mxEigenV))
break;
}
@@ -488,14 +493,14 @@ private int calcRotationMatrix() {
* commented block will never be activated. To be absolutely safe this
* should be uncommented, but it is most likely unnecessary.
*/
- if (qsqr < EVEC_PREC) {
+ if (qsqr < evec_prec) {
q1 = a12 * a3344_4334 - a13 * a3244_4234 + a14 * a3243_4233;
q2 = -a11 * a3344_4334 + a13 * a3144_4134 - a14 * a3143_4133;
q3 = a11 * a3244_4234 - a12 * a3144_4134 + a14 * a3142_4132;
q4 = -a11 * a3243_4233 + a12 * a3143_4133 - a13 * a3142_4132;
qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
- if (qsqr < EVEC_PREC) {
+ if (qsqr < evec_prec) {
double a1324_1423 = a13 * a24 - a14 * a23, a1224_1422 = a12
* a24 - a14 * a22;
double a1223_1322 = a12 * a23 - a13 * a22, a1124_1421 = a11
@@ -509,7 +514,7 @@ private int calcRotationMatrix() {
q4 = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221;
qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
- if (qsqr < EVEC_PREC) {
+ if (qsqr < evec_prec) {
q1 = a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322;
q2 = -a31 * a1324_1423 + a33 * a1124_1421 - a34
* a1123_1321;
@@ -518,7 +523,7 @@ private int calcRotationMatrix() {
* a1122_1221;
qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
- if (qsqr < EVEC_PREC) {
+ if (qsqr < evec_prec) {
/*
* if qsqr is still too small, return the identity
* matrix.
@@ -566,4 +571,21 @@ private int calcRotationMatrix() {
return 1;
}
+ @Override
+ public double getRmsd(Point3d[] fixed, Point3d[] moved) {
+ set(moved, fixed);
+ return getRmsd();
+ }
+
+ /**
+ * @param fixed
+ * @param moved
+ * @param weight array of weigths for each equivalent point position
+ * @return weighted RMSD.
+ */
+ public double getWeightedRmsd(Point3d[] fixed, Point3d[] moved, double[] weight) {
+ set(moved, fixed, weight);
+ return getRmsd();
+ }
+
}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQuat.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQuat.java
new file mode 100644
index 0000000000..1888c0d6b9
--- /dev/null
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQuat.java
@@ -0,0 +1,98 @@
+/*
+ * BioJava development 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.
+ *
+ * For more information on the BioJava project and its aims,
+ * or to join the biojava-l mailing list, visit the home page
+ * at:
+ *
+ * http://www.biojava.org/
+ *
+ */
+package org.biojava.nbio.structure.geometry;
+
+import javax.vecmath.*;
+
+/**
+ * The SuperPositionQuat implements a quaternion based algorithm to superpose
+ * arrays of Points in 3D.
+ *
+ * @author Peter Rose
+ * @author Aleix Lafita
+ *
+ */
+public final class SuperPositionQuat extends SuperPositionAbstract {
+
+ /**
+ * Constructor for the quaternion based superposition algorithm.
+ *
+ * @param centered
+ * true if the point arrays are centered at the origin (faster),
+ * false otherwise
+ */
+ public SuperPositionQuat(boolean centered) {
+ super(centered);
+ }
+
+ @Override
+ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
+
+ checkInput(fixed, moved);
+
+ if (centered) {
+ Quat4d q = UnitQuaternions.relativeOrientation(fixed, moved);
+ Matrix4d rotTrans = new Matrix4d();
+ rotTrans.set(q);
+ return rotTrans;
+ }
+
+ // translate to origin
+ Point3d[] xref = CalcPoint.clonePoint3dArray(fixed);
+ Point3d xtrans = CalcPoint.centroid(xref);
+ xtrans.negate();
+ CalcPoint.translate(xtrans, xref);
+
+ Point3d[] yref = CalcPoint.clonePoint3dArray(moved);
+ Point3d ytrans = CalcPoint.centroid(yref);
+ ytrans.negate();
+ CalcPoint.translate(ytrans, yref);
+
+ // calculate rotational component (rotation around origin)
+ Quat4d q = UnitQuaternions.relativeOrientation(xref, yref);
+ Matrix4d rotTrans = new Matrix4d();
+ rotTrans.set(q);
+
+ // combine with moved -> origin translation
+ Matrix4d trans = new Matrix4d();
+ trans.setIdentity();
+ trans.setTranslation(new Vector3d(ytrans));
+ rotTrans.mul(rotTrans, trans);
+
+ // combine with origin -> fixed translation
+ xtrans.negate();
+ Matrix4d transInverse = new Matrix4d();
+ transInverse.setIdentity();
+ transInverse.setTranslation(new Vector3d(xtrans));
+ rotTrans.mul(transInverse, rotTrans);
+
+ return rotTrans;
+
+ }
+
+ @Override
+ public double getRmsd(Point3d[] x, Point3d[] y) {
+ Point3d[] yclone = CalcPoint.clonePoint3dArray(y);
+ superposeAndTransform(x, yclone);
+ return CalcPoint.rmsd(x, yclone);
+ }
+
+}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionSVD.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionSVD.java
new file mode 100644
index 0000000000..aa91d254bd
--- /dev/null
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionSVD.java
@@ -0,0 +1,149 @@
+/*
+ * BioJava development 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.
+ *
+ * For more information on the BioJava project and its aims,
+ * or to join the biojava-l mailing list, visit the home page
+ * at:
+ *
+ * http://www.biojava.org/
+ *
+ * Created on Dec 4, 2005
+ *
+ */
+package org.biojava.nbio.structure.geometry;
+
+import javax.vecmath.Matrix4d;
+import javax.vecmath.Point3d;
+
+import org.biojava.nbio.structure.Atom;
+import org.biojava.nbio.structure.Calc;
+import org.biojava.nbio.structure.jama.Matrix;
+import org.biojava.nbio.structure.jama.SingularValueDecomposition;
+
+/**
+ * A class that calculates the superposition between two sets of points using an
+ * SVD Matrix Decomposition. It was introduced by Wolfgang Kabsch, hence the
+ * alternative name Kabsh algorithm. Inspired by the biopython SVDSuperimposer
+ * class.
+ *
+ * @author Andreas Prlic
+ * @author Aleix Lafita
+ * @since 1.5
+ * @version %I% %G%
+ *
+ */
+public class SuperPositionSVD extends SuperPositionAbstract {
+
+ /**
+ * Constructor for the SVD superposition algorithm.
+ *
+ * @param centered
+ * true if the point arrays are centered at the origin (faster),
+ * false otherwise
+ */
+ public SuperPositionSVD(boolean centered) {
+ super(centered);
+ }
+
+ public Matrix4d superpose(Atom[] fixed, Atom[] moved) {
+ return superpose(Calc.atomsToPoints(fixed), Calc.atomsToPoints(moved));
+ }
+
+ @Override
+ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
+
+ checkInput(fixed, moved);
+
+ Point3d cena = CalcPoint.centroid(fixed);
+ Point3d cenb = CalcPoint.centroid(moved);
+
+ double[][] centAcoords = new double[][] { { cena.x, cena.y, cena.z } };
+ Matrix centroidA = new Matrix(centAcoords);
+
+ double[][] centBcoords = new double[][] { { cenb.x, cenb.x, cenb.x } };
+ Matrix centroidB = new Matrix(centBcoords);
+
+ // center at centroid
+
+ cena.negate();
+ cenb.negate();
+
+ Point3d[] ats1 = CalcPoint.clonePoint3dArray(fixed);
+ CalcPoint.translate(cena, ats1);
+
+ Point3d[] ats2 = CalcPoint.clonePoint3dArray(moved);
+ CalcPoint.translate(cenb, ats2);
+
+ double[][] coordSet1 = new double[ats1.length][3];
+ double[][] coordSet2 = new double[ats2.length][3];
+
+ // copy the atoms into the internal coords;
+ for (int i = 0; i < ats1.length; i++) {
+ coordSet1[i] = new double[3];
+ ats1[i].get(coordSet1[i]);
+ coordSet2[i] = new double[3];
+ ats2[i].get(coordSet2[i]);
+ }
+
+ // now this is the bridge to the Jama package:
+ Matrix a = new Matrix(coordSet1);
+ Matrix b = new Matrix(coordSet2);
+
+ // # correlation matrix
+
+ Matrix b_trans = b.transpose();
+ Matrix corr = b_trans.times(a);
+
+ SingularValueDecomposition svd = corr.svd();
+
+ Matrix u = svd.getU();
+ // v is alreaady transposed ! difference to numermic python ...
+ Matrix vt = svd.getV();
+
+ Matrix vt_orig = (Matrix) vt.clone();
+ Matrix u_transp = u.transpose();
+
+ Matrix rot_nottrans = vt.times(u_transp);
+ Matrix rot = rot_nottrans.transpose();
+
+ // check if we have found a reflection
+
+ double det = rot.det();
+
+ if (det < 0) {
+ vt = vt_orig.transpose();
+ vt.set(2, 0, (0 - vt.get(2, 0)));
+ vt.set(2, 1, (0 - vt.get(2, 1)));
+ vt.set(2, 2, (0 - vt.get(2, 2)));
+
+ Matrix nv_transp = vt.transpose();
+ rot_nottrans = nv_transp.times(u_transp);
+ rot = rot_nottrans.transpose();
+
+ }
+
+ Matrix cb_tmp = centroidB.times(rot);
+ Matrix tran = centroidA.minus(cb_tmp);
+
+ return Matrices.getTransformation(rot, tran);
+
+ }
+
+ @Override
+ public double getRmsd(Point3d[] x, Point3d[] y) {
+ Point3d[] yclone = CalcPoint.clonePoint3dArray(y);
+ superposeAndTransform(x, yclone);
+ return CalcPoint.rmsd(x, yclone);
+ }
+
+}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositions.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositions.java
new file mode 100644
index 0000000000..56a788a669
--- /dev/null
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositions.java
@@ -0,0 +1,87 @@
+package org.biojava.nbio.structure.geometry;
+
+import javax.vecmath.Matrix4d;
+import javax.vecmath.Point3d;
+
+/**
+ * SuperPositions is a Class that provides static helper methods and an easy
+ * access to the whole family of {@link SuperPosition} algorithms.
+ *
+ * It defines a static SuperPosition object and uses it for calculation.
+ *
+ * @author Aleix Lafita
+ * @since 5.0.0
+ *
+ */
+public class SuperPositions {
+
+ private static SuperPositionAbstract superposer = new SuperPositionQuat(
+ false);
+
+ /** Prevent instantiation */
+ private SuperPositions(){}
+
+ /**
+ * Use the {@link SuperPosition#superpose(Point3d[], Point3d[])} method of
+ * the default static SuperPosition algorithm contained in this Class.
+ */
+ public static Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
+ superposer.setCentered(false);
+ return superposer.superpose(fixed, moved);
+ }
+
+ /**
+ * Use the {@link SuperPosition#superpose(Point3d[], Point3d[])} method of
+ * the default static SuperPosition algorithm contained in this Class,
+ * assuming that the point arrays are centered at the origin.
+ */
+ public static Matrix4d superposeAtOrigin(Point3d[] fixed, Point3d[] moved) {
+ superposer.setCentered(true);
+ return superposer.superpose(fixed, moved);
+ }
+
+ /**
+ * Use the {@link SuperPosition#superposeAndTransform(Point3d[], Point3d[])}
+ * method of the default static SuperPosition algorithm contained in this
+ * Class.
+ */
+ public static Matrix4d superposeAndTransform(Point3d[] fixed,
+ Point3d[] moved) {
+ superposer.setCentered(false);
+ return superposer.superposeAndTransform(fixed, moved);
+ }
+
+ /**
+ * Use the {@link SuperPosition#superposeAndTransform(Point3d[], Point3d[])}
+ * method of the default static SuperPosition algorithm contained in this
+ * Class, assuming that the point arrays are centered at the origin.
+ */
+ public static Matrix4d superposeAndTransformAtOrigin(Point3d[] fixed,
+ Point3d[] moved) {
+ superposer.setCentered(true);
+ return superposer.superposeAndTransform(fixed, moved);
+ }
+
+ /**
+ * Use the {@link SuperPosition#getRmsd(Point3d[], Point3d[])} method of the
+ * default static SuperPosition algorithm contained in this Class.
+ */
+ public static double getRmsd(Point3d[] fixed, Point3d[] moved) {
+ superposer.setCentered(false);
+ return superposer.getRmsd(fixed, moved);
+ }
+
+ /**
+ * Use the {@link SuperPosition#getRmsd(Point3d[], Point3d[])} method of the
+ * default static SuperPosition algorithm contained in this Class, assuming
+ * that the point arrays are centered at the origin.
+ */
+ public static double getRmsdAtOrigin(Point3d[] fixed, Point3d[] moved) {
+ superposer.setCentered(true);
+ return superposer.getRmsd(fixed, moved);
+ }
+
+ public static void setDefaultSuperPosition(SuperPositionAbstract defaultAlgorithm) {
+ superposer = defaultAlgorithm;
+ }
+}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/UnitQuaternions.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/UnitQuaternions.java
index bacb845eca..9b780bb942 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/UnitQuaternions.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/UnitQuaternions.java
@@ -125,14 +125,18 @@ public static double angle(Quat4d q) {
* The arrays of points need to be centered at the origin. To center the
* points use {@link CalcPoint#center(Point3d[])}.
*
- * @param a
- * array of Point3d, centered at origin
- * @param b
- * array of Point3d, centered at origin
- * @return the angle in radians of the relative orientation of the points
+ * @param fixed
+ * array of Point3d, centered at origin. Original coordinates
+ * will not be modified.
+ * @param moved
+ * array of Point3d, centered at origin. Original coordinates
+ * will not be modified.
+ * @return the angle in radians of the relative orientation of the points,
+ * angle to rotate moved to bring it to the same orientation as
+ * fixed.
*/
- public static double orientationAngle(Point3d[] a, Point3d[] b) {
- Quat4d q = relativeOrientation(a, b);
+ public static double orientationAngle(Point3d[] fixed, Point3d[] moved) {
+ Quat4d q = relativeOrientation(fixed, moved);
return angle(q);
}
@@ -141,36 +145,39 @@ public static double orientationAngle(Point3d[] a, Point3d[] b) {
* Equivalent to {@link #angle(Quat4d)} of the unit quaternion obtained by
* {@link #relativeOrientation(Point3d[], Point3d[])}.
*
- * @param a
- * array of Point3d
- * @param b
- * array of Point3d
+ * @param fixed
+ * array of Point3d. Original coordinates will not be modified.
+ * @param moved
+ * array of Point3d. Original coordinates will not be modified.
* @param centered
* true if the points are already centered at the origin
- * @return the angle in radians of the relative orientation of the points
+ * @return the angle in radians of the relative orientation of the points,
+ * angle to rotate moved to bring it to the same orientation as
+ * fixed.
*/
- public static double orientationAngle(Point3d[] a, Point3d[] b,
+ public static double orientationAngle(Point3d[] fixed, Point3d[] moved,
boolean centered) {
if (!centered) {
- a = CalcPoint.clonePoint3dArray(a);
- b = CalcPoint.clonePoint3dArray(b);
- CalcPoint.center(a);
- CalcPoint.center(b);
+ fixed = CalcPoint.clonePoint3dArray(fixed);
+ moved = CalcPoint.clonePoint3dArray(moved);
+ CalcPoint.center(fixed);
+ CalcPoint.center(moved);
}
- return orientationAngle(a, b);
+ return orientationAngle(fixed, moved);
}
/**
* Calculate the relative quaternion orientation of two arrays of points.
*
- * @param a
- * point array
- * @param b
- * point array
- * @return a unit quaternion representing the relative orientation
+ * @param fixed
+ * point array, coordinates will not be modified
+ * @param moved
+ * point array, coordinates will not be modified
+ * @return a unit quaternion representing the relative orientation, to
+ * rotate moved to bring it to the same orientation as fixed.
*/
- public static Quat4d relativeOrientation(Point3d[] a, Point3d[] b) {
- Matrix m = CalcPoint.formMatrix(a, b);
+ public static Quat4d relativeOrientation(Point3d[] fixed, Point3d[] moved) {
+ Matrix m = CalcPoint.formMatrix(moved, fixed); // inverse
EigenvalueDecomposition eig = m.eig();
double[][] v = eig.getV().getArray();
Quat4d q = new Quat4d(v[1][3], v[2][3], v[3][3], v[0][3]);
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/HelixAxisAligner.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/HelixAxisAligner.java
index 037dd8e277..e90e657a22 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/HelixAxisAligner.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/HelixAxisAligner.java
@@ -20,17 +20,25 @@
*/
package org.biojava.nbio.structure.symmetry.axis;
+import org.biojava.nbio.structure.geometry.CalcPoint;
import org.biojava.nbio.structure.geometry.MomentsOfInertia;
-import org.biojava.nbio.structure.geometry.SuperPosition;
import org.biojava.nbio.structure.symmetry.core.Helix;
import org.biojava.nbio.structure.symmetry.core.HelixLayers;
import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults;
import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
import javax.vecmath.*;
import java.util.*;
public class HelixAxisAligner extends AxisAligner {
+
+ private static final Logger logger = LoggerFactory
+ .getLogger(HelixAxisAligner.class);
+
+
private static final Vector3d Y_AXIS = new Vector3d(0,1,0);
private static final Vector3d Z_AXIS = new Vector3d(0,0,1);
@@ -483,8 +491,9 @@ private Matrix4d alignAxes(Vector3d[] axisVectors, Vector3d[] referenceVectors)
Point3d[] ref = new Point3d[2];
ref[0] = new Point3d(referenceVectors[0]);
ref[1] = new Point3d(referenceVectors[1]);
- if (SuperPosition.rmsd(axes, ref) > 0.1) {
- System.out.println("Warning: AxisTransformation: axes alignment is off. RMSD: " + SuperPosition.rmsd(axes, ref));
+ if (CalcPoint.rmsd(axes, ref) > 0.1) {
+ logger.warn("AxisTransformation: axes alignment is off. RMSD: "
+ + CalcPoint.rmsd(axes, ref));
}
return m2;
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/RotationAxisAligner.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/RotationAxisAligner.java
index 8e1e73fa53..9dbb87f3cd 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/RotationAxisAligner.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/axis/RotationAxisAligner.java
@@ -20,18 +20,24 @@
*/
package org.biojava.nbio.structure.symmetry.axis;
+import org.biojava.nbio.structure.geometry.CalcPoint;
import org.biojava.nbio.structure.geometry.MomentsOfInertia;
-import org.biojava.nbio.structure.geometry.SuperPosition;
import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults;
import org.biojava.nbio.structure.symmetry.core.Rotation;
import org.biojava.nbio.structure.symmetry.core.RotationGroup;
import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
import javax.vecmath.*;
import java.util.*;
public class RotationAxisAligner extends AxisAligner{
+
+ private static final Logger logger = LoggerFactory
+ .getLogger(RotationAxisAligner.class);
+
private static final Vector3d X_AXIS = new Vector3d(1,0,0);
private static final Vector3d Y_AXIS = new Vector3d(0,1,0);
private static final Vector3d Z_AXIS = new Vector3d(0,0,1);
@@ -463,8 +469,9 @@ private Matrix4d alignAxes(Vector3d[] axisVectors, Vector3d[] referenceVectors)
Point3d[] ref = new Point3d[2];
ref[0] = new Point3d(referenceVectors[0]);
ref[1] = new Point3d(referenceVectors[1]);
- if (SuperPosition.rmsd(axes, ref) > 0.1) {
- System.out.println("Warning: AxisTransformation: axes alignment is off. RMSD: " + SuperPosition.rmsd(axes, ref));
+ if (CalcPoint.rmsd(axes, ref) > 0.1) {
+ logger.warn("AxisTransformation: axes alignment is off. RMSD: "
+ + CalcPoint.rmsd(axes, ref));
}
return m2;
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/C2RotationSolver.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/C2RotationSolver.java
index b691dda050..a001d6462c 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/C2RotationSolver.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/C2RotationSolver.java
@@ -22,11 +22,12 @@
package org.biojava.nbio.structure.symmetry.core;
import org.biojava.nbio.structure.geometry.CalcPoint;
-import org.biojava.nbio.structure.geometry.SuperPosition;
+import org.biojava.nbio.structure.geometry.UnitQuaternions;
import javax.vecmath.AxisAngle4d;
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3d;
+import javax.vecmath.Quat4d;
import javax.vecmath.Vector3d;
import java.util.ArrayList;
@@ -78,9 +79,29 @@ private void solve() {
Point3d[] y = CalcPoint.clonePoint3dArray(traces.get(1));
CalcPoint.translate(new Point3d(trans), y);
+ // TODO implement this piece of code using at origin superposition
+ Quat4d quat = UnitQuaternions.relativeOrientation(
+ x, y);
AxisAngle4d axisAngle = new AxisAngle4d();
-
- Matrix4d transformation = SuperPosition.superposeAtOrigin(x, y, axisAngle);
+ Matrix4d transformation = new Matrix4d();
+
+ transformation.set(quat);
+ axisAngle.set(quat);
+
+ Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
+ if (axis.lengthSquared() < 1.0E-6) {
+ axisAngle.x = 0;
+ axisAngle.y = 0;
+ axisAngle.z = 1;
+ axisAngle.angle = 0;
+ } else {
+ axis.normalize();
+ axisAngle.x = axis.x;
+ axisAngle.y = axis.y;
+ axisAngle.z = axis.z;
+ }
+
+ CalcPoint.transform(transformation, y);
// if rmsd or angle deviation is above threshold, stop
double angleThresholdRadians = Math.toRadians(parameters.getAngleThreshold());
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/HelixSolver.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/HelixSolver.java
index 9ff28c64dd..9d60e12485 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/HelixSolver.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/HelixSolver.java
@@ -21,7 +21,7 @@
package org.biojava.nbio.structure.symmetry.core;
import org.biojava.nbio.structure.geometry.CalcPoint;
-import org.biojava.nbio.structure.geometry.SuperPosition;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -90,17 +90,14 @@ private void solve() {
logger.debug("HelixSolver: pair: " + Arrays.toString(pair));
int contacts = entry.getValue();
- Point3d[] h1 = null;
- Point3d[] h2 = null;
+ Point3d[] h1 = CalcPoint.clonePoint3dArray(repeatUnits.get(pair[0]));
+ Point3d[] h2 = CalcPoint.clonePoint3dArray(repeatUnits.get(pair[1]));
// trial superposition of repeat unit pairs to get a seed
// permutation
- h1 = CalcPoint.clonePoint3dArray(repeatUnits.get(pair[0]));
- h2 = CalcPoint.clonePoint3dArray(repeatUnits.get(pair[1]));
- Matrix4d transformation = SuperPosition.superposeWithTranslation(
- h1, h2);
+ Matrix4d transformation = SuperPositions.superposeAndTransform(h2, h1);
- double rmsd = SuperPosition.rmsd(h1, h2);
+ double rmsd = CalcPoint.rmsd(h1, h2);
double rise = getRise(transformation,
repeatUnitCenters.get(pair[0]),
repeatUnitCenters.get(pair[1]));
@@ -185,9 +182,9 @@ private void solve() {
// calculate subunit rmsd if at least 3 subunits are available
double subunitRmsd = 0;
if (point1.size() > 2) {
- transformation = SuperPosition.superposeWithTranslation(h1, h2);
+ transformation = SuperPositions.superposeAndTransform(h2, h1);
- subunitRmsd = SuperPosition.rmsd(h1, h2);
+ subunitRmsd = CalcPoint.rmsd(h1, h2);
rise = getRise(transformation, repeatUnitCenters.get(pair[0]),
repeatUnitCenters.get(pair[1]));
angle = getAngle(transformation);
@@ -229,13 +226,13 @@ private void solve() {
point1.toArray(h1);
point2.toArray(h2);
Point3d[] h3 = CalcPoint.clonePoint3dArray(h1);
- transformation = SuperPosition.superposeWithTranslation(h1, h2);
+ transformation = SuperPositions.superposeAndTransform(h2, h1);
Point3d xtrans = CalcPoint.centroid(h3);
xtrans.negate();
- double traceRmsd = SuperPosition.rmsd(h1, h2);
+ double traceRmsd = CalcPoint.rmsd(h1, h2);
rise = getRise(transformation, repeatUnitCenters.get(pair[0]),
repeatUnitCenters.get(pair[1]));
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/QuatSuperpositionScorer.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/QuatSuperpositionScorer.java
index e1e87d0c05..5df1278563 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/QuatSuperpositionScorer.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/QuatSuperpositionScorer.java
@@ -21,7 +21,7 @@
package org.biojava.nbio.structure.symmetry.core;
-import org.biojava.nbio.structure.geometry.SuperPosition;
+import org.biojava.nbio.structure.geometry.SuperPositions;
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3d;
@@ -150,7 +150,7 @@ private static void calcIntrasubunitScores(QuatSymmetrySubunits subunits, Matrix
}
// superpose individual subunits
- SuperPosition.superposeWithTranslation(trans, orig);
+ SuperPositions.superposeAndTransform(orig, trans);
for (int j = 0; j < orig.length; j++) {
double dSq = orig[j].distanceSquared(trans[j]);
sumTm += 1.0/(1.0 + dSq/d0Sq);
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java
index e2645c2201..2a879df3d1 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java
@@ -21,14 +21,16 @@
package org.biojava.nbio.structure.symmetry.core;
+import org.biojava.nbio.structure.geometry.CalcPoint;
import org.biojava.nbio.structure.geometry.MomentsOfInertia;
-import org.biojava.nbio.structure.geometry.SuperPosition;
+import org.biojava.nbio.structure.geometry.UnitQuaternions;
import org.biojava.nbio.structure.symmetry.geometry.DistanceBox;
import org.biojava.nbio.structure.symmetry.geometry.SphereSampler;
import javax.vecmath.AxisAngle4d;
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3d;
+import javax.vecmath.Quat4d;
import javax.vecmath.Vector3d;
import java.util.ArrayList;
@@ -166,10 +168,31 @@ private boolean evaluatePermutation(List permutation) {
int fold = PermutationGroup.getOrder(permutation);
- // get optimal transformation and axisangle by superimposing subunits
+ // get optimal transformation and axisangle by subunit superposition
+ // TODO implement this piece of code using at origin superposition
+ Quat4d quat = UnitQuaternions.relativeOrientation(
+ originalCoords, transformedCoords);
AxisAngle4d axisAngle = new AxisAngle4d();
- Matrix4d transformation = SuperPosition.superposeAtOrigin(transformedCoords, originalCoords, axisAngle);
- double subunitRmsd = SuperPosition.rmsd(transformedCoords, originalCoords);
+ Matrix4d transformation = new Matrix4d();
+
+ transformation.set(quat);
+ axisAngle.set(quat);
+
+ Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
+ if (axis.lengthSquared() < 1.0E-6) {
+ axisAngle.x = 0;
+ axisAngle.y = 0;
+ axisAngle.z = 1;
+ axisAngle.angle = 0;
+ } else {
+ axis.normalize();
+ axisAngle.x = axis.x;
+ axisAngle.y = axis.y;
+ axisAngle.z = axis.z;
+ }
+
+ CalcPoint.transform(transformation, transformedCoords);
+ double subunitRmsd = CalcPoint.rmsd(transformedCoords, originalCoords);
if (subunitRmsd < parameters.getRmsdThreshold()) {
combineWithTranslation(transformation);
@@ -252,7 +275,6 @@ private boolean isAllowedPermutation(List permutation) {
}
/**
* Adds translational component to rotation matrix
- * @param rotTrans
* @param rotation
* @return
*/
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/SystematicSolver.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/SystematicSolver.java
index fbd351c4c8..64fae2e92b 100644
--- a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/SystematicSolver.java
+++ b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/SystematicSolver.java
@@ -21,12 +21,14 @@
package org.biojava.nbio.structure.symmetry.core;
-import org.biojava.nbio.structure.geometry.SuperPosition;
+import org.biojava.nbio.structure.geometry.CalcPoint;
+import org.biojava.nbio.structure.geometry.UnitQuaternions;
import org.biojava.nbio.structure.symmetry.utils.PermutationGenerator;
import javax.vecmath.AxisAngle4d;
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3d;
+import javax.vecmath.Quat4d;
import javax.vecmath.Vector3d;
import java.util.ArrayList;
@@ -188,10 +190,32 @@ private boolean evaluatePermutation(List permutation) {
}
int fold = PermutationGroup.getOrder(permutation);
- // get optimal transformation and axisangle by superimposing subunits
+
+ // TODO implement this piece of code using at origin superposition
+ Quat4d quat = UnitQuaternions.relativeOrientation(
+ originalCoords, transformedCoords);
AxisAngle4d axisAngle = new AxisAngle4d();
- Matrix4d transformation = SuperPosition.superposeAtOrigin(transformedCoords, originalCoords, axisAngle);
- double subunitRmsd = SuperPosition.rmsd(transformedCoords, originalCoords);
+ Matrix4d transformation = new Matrix4d();
+
+ transformation.set(quat);
+ axisAngle.set(quat);
+
+ Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
+ if (axis.lengthSquared() < 1.0E-6) {
+ axisAngle.x = 0;
+ axisAngle.y = 0;
+ axisAngle.z = 1;
+ axisAngle.angle = 0;
+ } else {
+ axis.normalize();
+ axisAngle.x = axis.x;
+ axisAngle.y = axis.y;
+ axisAngle.z = axis.z;
+ }
+
+ CalcPoint.transform(transformation, transformedCoords);
+
+ double subunitRmsd = CalcPoint.rmsd(transformedCoords, originalCoords);
if (subunitRmsd 0 && arr2.length > 0) {
- SVDSuperimposer svd = new SVDSuperimposer(arr1, arr2);
- Matrix4d axis = svd.getTransformation();
+ Matrix4d axis = SuperPositions.superpose(
+ Calc.atomsToPoints(arr1),
+ Calc.atomsToPoints(arr2));
axes.updateAxis(level, axis);
}
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/TestBioAssemblyIdentifier.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/TestBioAssemblyIdentifier.java
index e5f47ee399..845ea0f51c 100644
--- a/biojava-structure/src/test/java/org/biojava/nbio/structure/TestBioAssemblyIdentifier.java
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/TestBioAssemblyIdentifier.java
@@ -40,13 +40,13 @@ public void test() throws IOException, StructureException {
// first assembly
id = new BioAssemblyIdentifier("BIO:2ehz:1");
s = cache.getStructure(id);
- assertEquals("Number of models",8, s.nrModels());
- assertEquals("Number of chains per model",11, s.getChains(0).size());
+ assertEquals("Number of models",1, s.nrModels());
+ assertEquals("Number of chains",88, s.getChains().size());
// equivalent
id = new BioAssemblyIdentifier("BIO:2ehz");
s = cache.getStructure(id);
- assertEquals("Number of models",8, s.nrModels());
- assertEquals("Number of chains per model",1,s.getPolyChains(0).size());
+ assertEquals("Number of models",1, s.nrModels());
+ assertEquals("Number of chains",8,s.getPolyChains().size());
// No second
id = new BioAssemblyIdentifier("BIO:2ehz:2");
try {
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/TestCalc.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/TestCalc.java
index c84c1e30e4..ccea9fa2c5 100644
--- a/biojava-structure/src/test/java/org/biojava/nbio/structure/TestCalc.java
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/TestCalc.java
@@ -25,6 +25,7 @@
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3d;
+import org.biojava.nbio.structure.geometry.Matrices;
import org.biojava.nbio.structure.jama.Matrix;
import org.junit.Test;
@@ -121,7 +122,7 @@ public void testJamaTransformation() {
//Sample transform: calc transposes automatically the matrix
//because it is a pre-multiplication rotation matrix
- Matrix sampleR = Calc.getRotationMatrix(getSampleTransform());
+ Matrix sampleR = Matrices.getRotationJAMA(getSampleTransform());
Atom sampleT = Calc.getTranslationVector(getSampleTransform());
Calc.rotate(atom, sampleR);
Calc.shift(atom, sampleT);
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/align/multiple/util/TestMultipleAlignmentWriter.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/align/multiple/util/TestMultipleAlignmentWriter.java
index 27e8db5732..0f54c654fe 100644
--- a/biojava-structure/src/test/java/org/biojava/nbio/structure/align/multiple/util/TestMultipleAlignmentWriter.java
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/align/multiple/util/TestMultipleAlignmentWriter.java
@@ -26,7 +26,6 @@
import org.biojava.nbio.structure.StructureException;
import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
-import org.biojava.nbio.structure.align.multiple.MultipleAlignmentEnsemble;
import org.biojava.nbio.structure.align.multiple.TestSampleGenerator;
import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentWriter;
import org.junit.Test;
@@ -34,14 +33,13 @@
import static org.junit.Assert.*;
/**
- * Test the correctness of various Text outputs
- * for {@link MultipleAlignment}s.
+ * Test the correctness of various Text outputs for {@link MultipleAlignment}s.
+ *
* Currently tested:
- *
- FASTA
+ *
+ * - FASTA
*
- FatCat format
*
- Aligned Residues
- *
- Transformation Matrices
- *
- XML format
*
*
* @author Aleix Lafita
@@ -53,27 +51,24 @@ public class TestMultipleAlignmentWriter {
private MultipleAlignment alignment2;
/**
- * Build the alignments in common for every writer output,
- * so that they do not have to be created each time.
+ * Build the alignments in common for every writer output, so that they do
+ * not have to be created each time.
+ *
* @throws IOException
* @throws StructureException
*/
- public TestMultipleAlignmentWriter()
- throws StructureException, IOException{
-
+ public TestMultipleAlignmentWriter() throws StructureException, IOException {
+
alignment1 = TestSampleGenerator.testAlignment1();
alignment2 = TestSampleGenerator.testAlignment2();
}
-
-
@Test
- public void testFASTA1() throws IOException{
+ public void testFASTA1() throws IOException {
String result = MultipleAlignmentWriter.toFASTA(alignment1);
- FileReader file = new FileReader(
- "src/test/resources/testMSTA1.fasta");
+ FileReader file = new FileReader("src/test/resources/testMSTA1.fasta");
BufferedReader reader = new BufferedReader(file);
String line = null;
StringBuilder stringBuilder = new StringBuilder();
@@ -85,7 +80,7 @@ public void testFASTA1() throws IOException{
reader.close();
String expected = stringBuilder.toString();
- assertEquals(expected,result);
+ assertEquals(expected, result);
}
@Test
@@ -93,8 +88,7 @@ public void testFASTA2() throws IOException {
String result = MultipleAlignmentWriter.toFASTA(alignment2);
- FileReader file = new FileReader(
- "src/test/resources/testMSTA2.fasta");
+ FileReader file = new FileReader("src/test/resources/testMSTA2.fasta");
BufferedReader reader = new BufferedReader(file);
String line = null;
StringBuilder stringBuilder = new StringBuilder();
@@ -106,16 +100,15 @@ public void testFASTA2() throws IOException {
reader.close();
String expected = stringBuilder.toString();
- assertEquals(expected,result);
+ assertEquals(expected, result);
}
@Test
- public void testFatCat1() throws IOException{
+ public void testFatCat1() throws IOException {
String result = MultipleAlignmentWriter.toFatCat(alignment1);
- FileReader file = new FileReader(
- "src/test/resources/testMSTA1.fatcat");
+ FileReader file = new FileReader("src/test/resources/testMSTA1.fatcat");
BufferedReader reader = new BufferedReader(file);
String line = null;
StringBuilder stringBuilder = new StringBuilder();
@@ -127,16 +120,15 @@ public void testFatCat1() throws IOException{
reader.close();
String expected = stringBuilder.toString();
- assertEquals(expected,result);
+ assertEquals(expected, result);
}
@Test
- public void testFatCat2() throws IOException{
+ public void testFatCat2() throws IOException {
String result = MultipleAlignmentWriter.toFatCat(alignment2);
- FileReader file = new FileReader(
- "src/test/resources/testMSTA2.fatcat");
+ FileReader file = new FileReader("src/test/resources/testMSTA2.fatcat");
BufferedReader reader = new BufferedReader(file);
String line = null;
StringBuilder stringBuilder = new StringBuilder();
@@ -148,11 +140,11 @@ public void testFatCat2() throws IOException{
reader.close();
String expected = stringBuilder.toString();
- assertEquals(expected,result);
+ assertEquals(expected, result);
}
@Test
- public void testAlignedResidues1() throws IOException{
+ public void testAlignedResidues1() throws IOException {
String result = MultipleAlignmentWriter.toAlignedResidues(alignment1);
@@ -169,11 +161,11 @@ public void testAlignedResidues1() throws IOException{
reader.close();
String expected = stringBuilder.toString();
- assertEquals(expected,result);
+ assertEquals(expected, result);
}
@Test
- public void testAlignedResidues2() throws IOException{
+ public void testAlignedResidues2() throws IOException {
String result = MultipleAlignmentWriter.toAlignedResidues(alignment2);
@@ -190,98 +182,7 @@ public void testAlignedResidues2() throws IOException{
reader.close();
String expected = stringBuilder.toString();
- assertEquals(expected,result);
- }
-
- @Test
- public void testTransformMatrices1() throws IOException{
-
- String result = MultipleAlignmentWriter.
- toTransformMatrices(alignment1);
-
- FileReader file = new FileReader(
- "src/test/resources/testMSTA1.transforms");
- BufferedReader reader = new BufferedReader(file);
- String line = null;
- StringBuilder stringBuilder = new StringBuilder();
-
- while ((line = reader.readLine()) != null) {
- stringBuilder.append(line);
- stringBuilder.append("\n");
- }
- reader.close();
-
- String expected = stringBuilder.toString();
- assertEquals(expected,result);
- }
-
- @Test
- public void testTransformMatrices2() throws IOException{
-
- String result = MultipleAlignmentWriter.
- toTransformMatrices(alignment2);
-
-
- FileReader file = new FileReader(
- "src/test/resources/testMSTA2.transforms");
- BufferedReader reader = new BufferedReader(file);
- String line = null;
- StringBuilder stringBuilder = new StringBuilder();
-
- while ((line = reader.readLine()) != null) {
- stringBuilder.append(line);
- stringBuilder.append("\n");
- }
- reader.close();
-
- String expected = stringBuilder.toString();
- // TODO THIS SHOULD BE REWRITTEN TO USE FLOATS
- //assertEquals(expected,result);
-
+ assertEquals(expected, result);
}
- @Test
- public void testXMLformat1() throws IOException{
-
- MultipleAlignmentEnsemble ensemble = alignment1.getEnsemble();
- String result = MultipleAlignmentWriter.toXML(ensemble);
-
- FileReader file = new FileReader(
- "src/test/resources/testMSTA1.xml");
- BufferedReader reader = new BufferedReader(file);
- String line = null;
- StringBuilder stringBuilder = new StringBuilder();
-
- while ((line = reader.readLine()) != null) {
- stringBuilder.append(line);
- stringBuilder.append("\n");
- }
- reader.close();
-
- String expected = stringBuilder.toString();
- assertEquals(expected,result);
- }
-
- @Test
- public void testXMLformat2() throws IOException{
-
- MultipleAlignmentEnsemble ensemble = alignment2.getEnsemble();
- String result = MultipleAlignmentWriter.toXML(ensemble);
-
- FileReader file = new FileReader(
- "src/test/resources/testMSTA2.xml");
- BufferedReader reader = new BufferedReader(file);
- String line = null;
- StringBuilder stringBuilder = new StringBuilder();
-
- while ((line = reader.readLine()) != null) {
- stringBuilder.append(line);
- stringBuilder.append("\n");
- }
- reader.close();
-
- String expected = stringBuilder.toString();
- // TODO THIS SHOULD BE REWRITTEN TO DO FLOAT AND NOT TEXT COMPARISON
-// assertEquals(expected,result);
- }
}
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPosition.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPosition.java
new file mode 100644
index 0000000000..b7f693e311
--- /dev/null
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPosition.java
@@ -0,0 +1,167 @@
+package org.biojava.nbio.structure.geometry;
+
+import static org.junit.Assert.*;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+import java.util.Random;
+
+import javax.vecmath.AxisAngle4d;
+import javax.vecmath.Matrix4d;
+import javax.vecmath.Point3d;
+import javax.vecmath.Vector3d;
+
+import org.biojava.nbio.structure.StructureException;
+import org.biojava.nbio.structure.geometry.SuperPositionQuat;
+import org.biojava.nbio.structure.geometry.SuperPositionQCP;
+import org.junit.Before;
+import org.junit.Test;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+/**
+ * Test the Quaternion-Based Characteristic Polynomial {@link SuperPositionQCP}
+ * algorithm for RMSD and Superposition calculations.
+ *
+ * @author Aleix Lafita
+ * @since 5.0.0
+ *
+ */
+public class TestSuperPosition {
+
+ private static final Logger logger = LoggerFactory
+ .getLogger(TestSuperPosition.class);
+
+ private List cloud1;
+ private List cloud2;
+
+ private AxisAngle4d rotAxis;
+ private Vector3d translation;
+ private Matrix4d transform;
+
+ /**
+ * Generate two clouds of random points of different sizes to test
+ * correctness and performance of superposition algorithms.
+ *
+ * @throws StructureException
+ */
+ @Before
+ public void setUp() throws StructureException {
+
+ cloud1 = new ArrayList(5);
+ cloud2 = new ArrayList(5);
+
+ Random rnd = new Random(0);
+
+ rotAxis = new AxisAngle4d(0.440, 0.302, 0.845, 1.570);
+ translation = new Vector3d(0.345, 2.453, 5.324);
+ transform = new Matrix4d();
+ transform.set(rotAxis);
+ transform.setTranslation(translation);
+
+ List sizes = Arrays.asList(5, 50, 500, 5000, 50000, 500000);
+
+ for (Integer size : sizes) {
+
+ Point3d[] c1 = new Point3d[size];
+ Point3d[] c2 = new Point3d[size];
+
+ for (int p = 0; p < size; p++) {
+
+ Point3d a = new Point3d(rnd.nextInt(100), rnd.nextInt(50),
+ rnd.nextInt(150));
+ c1[p] = a;
+
+ // Add some noise
+ Point3d b = new Point3d(a.x + rnd.nextDouble(), a.y
+ + rnd.nextDouble(), a.z + rnd.nextDouble());
+ c2[p] = b;
+ }
+
+ CalcPoint.center(c1);
+ CalcPoint.center(c2);
+
+ CalcPoint.transform(transform, c1);
+
+ cloud1.add(c1);
+ cloud2.add(c2);
+ }
+
+ }
+
+ /**
+ * Test method to obtain the transformation matrix from superpositions.
+ */
+ @Test
+ public void testSuperposition() {
+
+ for (int c = 0; c < cloud1.size(); c++) {
+
+ // Use SVD superposition
+ SuperPosition svd = new SuperPositionSVD(false);
+ long svdStart = System.nanoTime();
+ Matrix4d svdTransform = svd.superpose(cloud1.get(c), cloud2.get(c));
+ long svdTime = (System.nanoTime() - svdStart) / 1000;
+
+ // Use quaternion superposition
+ SuperPosition quat = new SuperPositionQuat(false);
+ long quatStart = System.nanoTime();
+ Matrix4d quatTransform = quat.superpose(cloud1.get(c), cloud2.get(c));
+ long quatTime = (System.nanoTime() - quatStart) / 1000;
+
+ // Use QCP algorithm
+ SuperPosition qcp = new SuperPositionQCP(false);
+ long qcpStart = System.nanoTime();
+ Matrix4d qcpTransform = qcp.superpose(cloud1.get(c), cloud2.get(c));
+ long qcpTime = (System.nanoTime() - qcpStart) / 1000;
+
+ logger.error(String.format("Transformation Matrix %d points: "
+ + "SVD time %d us, SP time: %d us, QCP time: %d us",
+ cloud1.get(c).length, svdTime, quatTime, qcpTime));
+
+ // Check that the transformation matrix was recovered
+ assertTrue(transform.epsilonEquals(svdTransform, 0.01));
+ assertTrue(transform.epsilonEquals(quatTransform, 0.01));
+ assertTrue(transform.epsilonEquals(qcpTransform, 0.01));
+ }
+
+ }
+
+ /**
+ * Test method to obtain the RMSD of a superposition.
+ */
+ @Test
+ public void testRMSD() {
+
+ for (int c = 0; c < cloud1.size(); c++) {
+
+ // Use SVD superposition
+ SuperPosition svd = new SuperPositionSVD(false);
+ long svdStart = System.nanoTime();
+ double svdrmsd = svd.getRmsd(cloud1.get(c), cloud2.get(c));
+ long svdTime = (System.nanoTime() - svdStart) / 1000;
+
+ // Use quaternion superposition
+ SuperPosition quat = new SuperPositionQuat(false);
+ long quatStart = System.nanoTime();
+ double quatrmsd = quat.getRmsd(cloud1.get(c), cloud2.get(c));
+ long quatTime = (System.nanoTime() - quatStart) / 1000;
+
+ // Use QCP algorithm
+ SuperPosition qcp = new SuperPositionQCP(false);
+ long qcpStart = System.nanoTime();
+ double qcprmsd = qcp.getRmsd(cloud1.get(c), cloud2.get(c));
+ long qcpTime = (System.nanoTime() - qcpStart) / 1000;
+
+ logger.error(String.format("RMSD %d points: SVD time %d us, "
+ + "Quat time: %d us, QCP time: %d us", cloud1.get(c).length,
+ svdTime, quatTime, qcpTime));
+
+ // Check that the returned RMSDs are equal
+ assertEquals(svdrmsd, quatrmsd, 0.001);
+ assertEquals(svdrmsd, qcprmsd, 0.001);
+ }
+ }
+
+}
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPositionQCP.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPositionQCP.java
index e0d134a9a5..d815276919 100644
--- a/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPositionQCP.java
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPositionQCP.java
@@ -2,21 +2,13 @@
import static org.junit.Assert.*;
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.List;
-import java.util.Random;
-
import javax.vecmath.AxisAngle4d;
import javax.vecmath.Matrix4d;
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
-import org.biojava.nbio.structure.SVDSuperimposer;
-import org.biojava.nbio.structure.StructureException;
-import org.biojava.nbio.structure.geometry.SuperPosition;
+import org.biojava.nbio.structure.geometry.SuperPositionQuat;
import org.biojava.nbio.structure.geometry.SuperPositionQCP;
-import org.junit.Before;
import org.junit.Test;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -34,157 +26,14 @@ public class TestSuperPositionQCP {
private static final Logger logger = LoggerFactory
.getLogger(TestSuperPositionQCP.class);
- private List cloud1;
- private List cloud2;
-
- private AxisAngle4d rotAxis;
- private Vector3d translation;
- private Matrix4d transform;
-
- /**
- * Generate two clouds of random points of different sizes to test
- * correctness and performance of superposition algorithms.
- *
- * @throws StructureException
- */
- @Before
- public void setUp() throws StructureException {
-
- cloud1 = new ArrayList(5);
- cloud2 = new ArrayList(5);
-
- Random rnd = new Random(0);
-
- rotAxis = new AxisAngle4d(0.440, 0.302, 0.845, 1.570);
- translation = new Vector3d(0.345, 2.453, 5.324);
- transform = new Matrix4d();
- transform.set(rotAxis);
- transform.setTranslation(translation);
-
- List sizes = Arrays.asList(5, 50, 500, 5000, 50000, 500000);
-
- for (Integer size : sizes) {
-
- Point3d[] c1 = new Point3d[size];
- Point3d[] c2 = new Point3d[size];
-
- for (int p = 0; p < size; p++) {
-
- Point3d a = new Point3d(rnd.nextInt(100), rnd.nextInt(50),
- rnd.nextInt(150));
- c1[p] = a;
-
- // Add some noise
- Point3d b = new Point3d(a.x + rnd.nextDouble(), a.y
- + rnd.nextDouble(), a.z + rnd.nextDouble());
- c2[p] = b;
- }
-
- CalcPoint.center(c1);
- CalcPoint.center(c2);
-
- CalcPoint.transform(transform, c1);
-
- cloud1.add(c1);
- cloud2.add(c2);
- }
-
- }
-
- /**
- * Test method to obtain the transformation matrix from superposition
- * {@link SuperPositionQCP#getTransformationMatrix()},
- * {@link SuperPositionQCP#getRmsd()}.
- *
- * @throws StructureException
- */
- @Test
- public void testTransformationMatrix() throws StructureException {
-
- for (int c = 0; c < cloud1.size(); c++) {
- // Use SVD superposition to obtain the optimal transformation matrix
- long svdStart = System.nanoTime();
- SVDSuperimposer svd = new SVDSuperimposer(cloud1.get(c),
- cloud2.get(c));
- Matrix4d svdTransform = svd.getTransformation();
- long svdTime = (System.nanoTime() - svdStart) / 1000;
-
- Point3d[] c2c = CalcPoint.clonePoint3dArray(cloud2.get(c));
-
- // Use SuperPosition to obtain the optimal transformation matrix
- long spStart = System.nanoTime();
- Matrix4d spTransform = SuperPosition.superposeWithTranslation(c2c,
- cloud1.get(c));
- long spTime = (System.nanoTime() - spStart) / 1000;
-
- // Use QCP algorithm to get the optimal transformation matrix
- SuperPositionQCP qcp = new SuperPositionQCP();
- qcp.set(cloud2.get(c), cloud1.get(c));
- long qcpStart = System.nanoTime();
- Matrix4d qcpTransform = qcp.getTransformationMatrix();
- long qcpTime = (System.nanoTime() - qcpStart) / 1000;
-
- logger.info(String.format("Transformation Matrix %d points: "
- + "SVD time %d us, SP time: %d us, QCP time: %d us",
- cloud1.get(c).length, svdTime, spTime, qcpTime));
-
- // Check that the transformation matrix was recovered
- assertTrue(transform.epsilonEquals(svdTransform, 0.01));
- assertTrue(transform.epsilonEquals(spTransform, 0.01));
- assertTrue(transform.epsilonEquals(qcpTransform, 0.01));
- }
-
- }
-
/**
- * Test method to obtain the RMSD of a superposition
- * {@link SuperPositionQCP#getRmsd()}.
- *
- * @throws StructureException
+ * Test case proposed by Peter Rose from his observations about quaternary
+ * symmetry artifacts with the QCP algorithm.
*/
- @Test
- public void testRMSD() throws StructureException {
-
- for (int c = 0; c < cloud1.size(); c++) {
- Point3d[] c2c = CalcPoint.clonePoint3dArray(cloud2.get(c));
-
- // Use SVD superposition to obtain the RMSD
- long svdStart = System.nanoTime();
- SVDSuperimposer svd = new SVDSuperimposer(cloud1.get(c),
- cloud2.get(c));
- Matrix4d svdTransform = svd.getTransformation();
- CalcPoint.transform(svdTransform, c2c);
- double svdrmsd = SuperPosition.rmsd(cloud1.get(c), c2c);
- long svdTime = (System.nanoTime() - svdStart) / 1000;
-
- c2c = CalcPoint.clonePoint3dArray(cloud2.get(c));
-
- // Use SVD superposition to obtain the RMSD
- long spStart = System.nanoTime();
- SuperPosition.superposeWithTranslation(c2c, cloud1.get(c));
- double sprmsd = SuperPosition.rmsd(c2c, cloud1.get(c));
- long spTime = (System.nanoTime() - spStart) / 1000;
-
- // Use QCP algorithm to obtain the RMSD
- SuperPositionQCP qcp = new SuperPositionQCP();
- qcp.set(cloud2.get(c), cloud1.get(c));
- long qcpStart = System.nanoTime();
- double qcprmsd = qcp.getRmsd();
- long qcpTime = (System.nanoTime() - qcpStart) / 1000;
-
- logger.info(String.format("RMSD %d points: SVD time %d us, "
- + "SP time: %d us, QCP time: %d us", cloud1.get(c).length,
- svdTime, spTime, qcpTime));
-
- // Check that the returned RMSDs are equal
- assertEquals(svdrmsd, sprmsd, 0.001);
- assertEquals(svdrmsd, qcprmsd, 0.001);
- }
- }
-
@Test
public void testSymmetryQCP() {
+ // Generate an array of points with symmetry
Point3d[] set1 = new Point3d[16];
set1[0] = new Point3d(14.065934, 47.068832, -32.895836);
set1[1] = new Point3d(-14.065934, -47.068832, -32.895836);
@@ -204,28 +53,32 @@ public void testSymmetryQCP() {
set1[15] = new Point3d(-22.748293, -43.813946, 32.14434);
Point3d[] set2 = CalcPoint.clonePoint3dArray(set1);
- CalcPoint.transform(transform, set2);
- // Use SP superposition to obtain the RMSD
- long spStart = System.nanoTime();
- SuperPosition.superposeWithTranslation(set1, set2);
- double sprmsd = SuperPosition.rmsd(set1, set2);
- long spTime = (System.nanoTime() - spStart) / 1000;
+ // Use a random transformation to set2
+ AxisAngle4d rotAxis = new AxisAngle4d(0.440, 0.302, 0.845, 1.570);
+ Vector3d translation = new Vector3d(0.345, 2.453, 5.324);
+ Matrix4d transform = new Matrix4d();
+ transform.set(rotAxis);
+ transform.setTranslation(translation);
+ CalcPoint.transform(transform, set2);
- set2 = CalcPoint.clonePoint3dArray(set1);
+ // Use Quaternion superposition to obtain the RMSD
+ SuperPosition algorithm = new SuperPositionQuat(false);
+ long quatStart = System.nanoTime();
+ double quatrmsd = algorithm.getRmsd(set1, set2);
+ long quatTime = (System.nanoTime() - quatStart) / 1000;
// Use QCP algorithm to get the RMSD
- SuperPositionQCP qcp = new SuperPositionQCP();
- qcp.set(set1, set2);
+ algorithm = new SuperPositionQCP(false);
long qcpStart = System.nanoTime();
- double qcprmsd = qcp.getRmsd();
+ double qcprmsd = algorithm.getRmsd(set1, set2);
long qcpTime = (System.nanoTime() - qcpStart) / 1000;
- logger.info(String.format("RMSD Symmetry: SP time: %d us"
- + ", QCP time: %d us", spTime, qcpTime));
+ logger.info(String.format("RMSD Symmetry: Quat time: %d us"
+ + ", QCP time: %d us", quatTime, qcpTime));
// Check that the returned RMSDs are equal
- assertEquals(sprmsd, qcprmsd, 0.001);
+ assertEquals(quatrmsd, qcprmsd, 0.001);
}
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestUnitQuaternions.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestUnitQuaternions.java
index 1bbc19d629..5dfc35e6ba 100644
--- a/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestUnitQuaternions.java
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestUnitQuaternions.java
@@ -222,12 +222,12 @@ public void testRelativeOrientation() throws IOException,
double angle = UnitQuaternions.orientationAngle(cloud, cloud2, false);
assertEquals(angle, 0, 0.001);
- // Apply a 30 degree rotation to cloud 2
+ // Apply a 30 degree rotation to cloud
AxisAngle4d axis = new AxisAngle4d(new Vector3d(1,1,1), Math.PI / 6);
Matrix4d transform = new Matrix4d();
transform.set(axis);
- CalcPoint.transform(transform, cloud2);
+ CalcPoint.transform(transform, cloud);
angle = UnitQuaternions.orientationAngle(cloud, cloud2, false);
angle = Math.min(Math.abs(2 * Math.PI - angle), angle);
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/io/TestHardBioUnits.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/io/TestHardBioUnits.java
index e513906c8e..81f62f6561 100644
--- a/biojava-structure/src/test/java/org/biojava/nbio/structure/io/TestHardBioUnits.java
+++ b/biojava-structure/src/test/java/org/biojava/nbio/structure/io/TestHardBioUnits.java
@@ -63,18 +63,12 @@ public void test4A1I() throws IOException, StructureException {
//System.out.println(bioAssembly.toPDB());
- assertEquals(bioAssembly.nrModels(), 2);
+ assertEquals(1, bioAssembly.nrModels());
- assertTrue(bioAssembly.getChains().size() > 0);
+ assertEquals(2, bioAssembly.getPolyChains().size());
- // Check that one model has Chain B and the other has chain G
- Chain g = bioAssembly.getPolyChainByPDB("G", 0);
- Chain b = bioAssembly.getPolyChainByPDB("B", 1);
- // Anthony Bradley - I think the model numbering here is arbitrary???
- if(g==null){
- g = bioAssembly.getPolyChainByPDB("G", 1);
- b = bioAssembly.getPolyChainByPDB("B", 0);
- }
+ Chain g = bioAssembly.getPolyChainByPDB("G_1");
+ Chain b = bioAssembly.getPolyChainByPDB("B_2");
assertNotNull(g);
@@ -84,8 +78,6 @@ public void test4A1I() throws IOException, StructureException {
assertFalse(bioAssembly.hasChain("H"));
- assertEquals(1,bioAssembly.getPolyChains(0).size());
- assertEquals(1,bioAssembly.getPolyChains(1).size());
}
diff --git a/biojava-structure/src/test/resources/testMSTA1.transforms b/biojava-structure/src/test/resources/testMSTA1.transforms
deleted file mode 100644
index c084de4591..0000000000
--- a/biojava-structure/src/test/resources/testMSTA1.transforms
+++ /dev/null
@@ -1,26 +0,0 @@
-Operations for block 1
- X1 = ( 1.000000)*Xref + ( 0.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Y1 = ( 0.000000)*Xref + ( 1.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Z1 = ( 0.000000)*Xref + ( 0.000000)*Yref + ( 1.000000)*Zref + ( 0.000000)
-
- X2 = ( 1.000000)*Xref + (-0.000000)*Yref + (-0.000000)*Zref + ( -0.000000)
- Y2 = (-0.000000)*Xref + ( 1.000000)*Yref + (-0.000000)*Zref + ( 0.000000)
- Z2 = ( 0.000000)*Xref + ( 0.000000)*Yref + ( 1.000000)*Zref + ( 0.000000)
-
- X3 = ( 0.994348)*Xref + ( 0.098958)*Yref + ( 0.038458)*Zref + ( 3.279491)
- Y3 = (-0.102907)*Xref + ( 0.987440)*Yref + ( 0.119883)*Zref + ( -1.341072)
- Z3 = (-0.026112)*Xref + (-0.123163)*Yref + ( 0.992043)*Zref + ( -4.769127)
-
-Operations for block 2
- X1 = ( 1.000000)*Xref + ( 0.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Y1 = ( 0.000000)*Xref + ( 1.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Z1 = ( 0.000000)*Xref + ( 0.000000)*Yref + ( 1.000000)*Zref + ( 0.000000)
-
- X2 = ( 0.980507)*Xref + (-0.035883)*Yref + ( 0.193180)*Zref + ( -2.815110)
- Y2 = ( 0.022123)*Xref + ( 0.997092)*Yref + ( 0.072921)*Zref + ( -2.750880)
- Z2 = (-0.195235)*Xref + (-0.067225)*Yref + ( 0.978450)*Zref + ( -1.900286)
-
- X3 = ( 0.418403)*Xref + ( 0.836888)*Yref + (-0.352926)*Zref + ( 46.918691)
- Y3 = (-0.834517)*Xref + ( 0.507594)*Yref + ( 0.214310)*Zref + ( -15.835306)
- Z3 = ( 0.358497)*Xref + ( 0.204855)*Yref + ( 0.910777)*Zref + ( 2.596739)
-
diff --git a/biojava-structure/src/test/resources/testMSTA1.xml b/biojava-structure/src/test/resources/testMSTA1.xml
deleted file mode 100644
index c493e13aea..0000000000
--- a/biojava-structure/src/test/resources/testMSTA1.xml
+++ /dev/null
@@ -1,74 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
\ No newline at end of file
diff --git a/biojava-structure/src/test/resources/testMSTA2.transforms b/biojava-structure/src/test/resources/testMSTA2.transforms
deleted file mode 100644
index 8ee3cdac59..0000000000
--- a/biojava-structure/src/test/resources/testMSTA2.transforms
+++ /dev/null
@@ -1,34 +0,0 @@
-Operations for block 1
- X1 = ( 1.000000)*Xref + ( 0.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Y1 = ( 0.000000)*Xref + ( 1.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Z1 = ( 0.000000)*Xref + ( 0.000000)*Yref + ( 1.000000)*Zref + ( 0.000000)
-
- X2 = ( 0.814927)*Xref + (-0.239779)*Yref + ( 0.527636)*Zref + ( -0.198745)
- Y2 = ( 0.467768)*Xref + ( 0.809629)*Yref + (-0.354534)*Zref + ( 2.791076)
- Z2 = (-0.342180)*Xref + ( 0.535730)*Yref + ( 0.771950)*Zref + ( -8.230708)
-
- X3 = (-0.330341)*Xref + ( 0.933003)*Yref + (-0.142761)*Zref + ( 6.552893)
- Y3 = ( 0.678439)*Xref + ( 0.339868)*Yref + ( 0.651314)*Zref + ( 1.762178)
- Z3 = ( 0.656198)*Xref + ( 0.118301)*Yref + (-0.745258)*Zref + ( 3.834594)
-
- X4 = ( 0.518511)*Xref + (-0.228254)*Yref + (-0.824043)*Zref + ( 25.789098)
- Y4 = ( 0.499493)*Xref + ( 0.863045)*Yref + ( 0.075237)*Zref + ( -19.128131)
- Z4 = ( 0.694013)*Xref + (-0.450614)*Yref + ( 0.561509)*Zref + ( 12.348022)
-
-Operations for block 2
- X1 = ( 1.000000)*Xref + ( 0.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Y1 = ( 0.000000)*Xref + ( 1.000000)*Yref + ( 0.000000)*Zref + ( 0.000000)
- Z1 = ( 0.000000)*Xref + ( 0.000000)*Yref + ( 1.000000)*Zref + ( 0.000000)
-
- X2 = ( 0.842373)*Xref + (-0.345647)*Yref + ( 0.413445)*Zref + ( 4.318617)
- Y2 = ( 0.521783)*Xref + ( 0.714942)*Yref + (-0.465403)*Zref + ( 7.141563)
- Z2 = (-0.134724)*Xref + ( 0.607772)*Yref + ( 0.782600)*Zref + ( -10.200162)
-
- X3 = (-0.359346)*Xref + ( 0.905779)*Yref + (-0.224576)*Zref + ( 6.158376)
- Y3 = ( 0.649531)*Xref + ( 0.415555)*Yref + ( 0.636729)*Zref + ( 1.430105)
- Z3 = ( 0.670059)*Xref + ( 0.082937)*Yref + (-0.737660)*Zref + ( 3.658757)
-
- X4 = ( 0.529785)*Xref + (-0.108909)*Yref + (-0.841110)*Zref + ( 21.801521)
- Y4 = ( 0.305436)*Xref + ( 0.949679)*Yref + ( 0.069417)*Zref + ( -20.677779)
- Z4 = ( 0.791225)*Xref + (-0.293682)*Yref + ( 0.536390)*Zref + ( 7.940602)
-
diff --git a/biojava-structure/src/test/resources/testMSTA2.xml b/biojava-structure/src/test/resources/testMSTA2.xml
deleted file mode 100644
index 93cf8b5d99..0000000000
--- a/biojava-structure/src/test/resources/testMSTA2.xml
+++ /dev/null
@@ -1,76 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
\ No newline at end of file