Skip to content

Commit ac19532

Browse files
sblivensbliven
committed
Fix Quaternary symmetry symmetry detection bug
In special cases it was possible to have pairs of axes that pass our quality thresholds but form invalid symmetry groups. Now alignments that superimpose a subunit onto itself are forbidden, which should restrict it to the chiral groups. Also add numerous notes and documentation.
1 parent 89ac11a commit ac19532

File tree

5 files changed

+192
-60
lines changed

5 files changed

+192
-60
lines changed

‎biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/PermutationGroup.java‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,14 +23,15 @@
2323

2424
import java.util.ArrayList;
2525
import java.util.HashSet;
26+
import java.util.Iterator;
2627
import java.util.List;
2728
import java.util.Set;
2829

2930
/**
3031
*
3132
* @author Peter
3233
*/
33-
public class PermutationGroup {
34+
public class PermutationGroup implements Iterable<List<Integer>> {
3435
List<List<Integer>> permutations = new ArrayList<List<Integer>>();
3536

3637
public void addPermutation(List<Integer> permutation) {
@@ -139,4 +140,9 @@ public String getGroupTable() {
139140
public int hashCode() {
140141
return getGroupTable().hashCode();
141142
}
143+
144+
@Override
145+
public Iterator<List<Integer>> iterator() {
146+
return permutations.iterator();
147+
}
142148
}

‎biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationGroup.java‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,13 +27,14 @@
2727
import java.util.ArrayList;
2828
import java.util.Collections;
2929
import java.util.Comparator;
30+
import java.util.Iterator;
3031
import java.util.List;
3132

3233
/**
3334
*
3435
* @author Peter
3536
*/
36-
public class RotationGroup {
37+
public class RotationGroup implements Iterable<Rotation> {
3738
private List<Rotation> rotations = new ArrayList<Rotation>();
3839
private int principalAxisIndex = 0;
3940
private int higherOrderRotationAxis = 0;
@@ -369,4 +370,9 @@ public int compare(Rotation o1, Rotation o2) {
369370
}
370371
});
371372
}
373+
374+
@Override
375+
public Iterator<Rotation> iterator() {
376+
return rotations.iterator();
377+
}
372378
}

‎biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/RotationSolver.java‎

Lines changed: 126 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -21,19 +21,22 @@
2121

2222
package org.biojava.nbio.structure.symmetry.core;
2323

24-
import org.biojava.nbio.structure.symmetry.geometry.DistanceBox;
25-
import org.biojava.nbio.structure.symmetry.geometry.MomentsOfInertia;
26-
import org.biojava.nbio.structure.symmetry.geometry.SphereSampler;
27-
import org.biojava.nbio.structure.symmetry.geometry.SuperPosition;
24+
import java.util.ArrayList;
25+
import java.util.HashMap;
26+
import java.util.HashSet;
27+
import java.util.List;
28+
import java.util.Map;
29+
import java.util.Set;
2830

2931
import javax.vecmath.AxisAngle4d;
3032
import javax.vecmath.Matrix4d;
3133
import javax.vecmath.Point3d;
3234
import javax.vecmath.Vector3d;
33-
import java.util.ArrayList;
34-
import java.util.HashSet;
35-
import java.util.List;
36-
import java.util.Set;
35+
36+
import org.biojava.nbio.structure.symmetry.geometry.DistanceBox;
37+
import org.biojava.nbio.structure.symmetry.geometry.MomentsOfInertia;
38+
import org.biojava.nbio.structure.symmetry.geometry.SphereSampler;
39+
import org.biojava.nbio.structure.symmetry.geometry.SuperPosition;
3740

3841

3942
/**
@@ -50,7 +53,8 @@ public class RotationSolver implements QuatSymmetrySolver {
5053
private Matrix4d centroidInverse = new Matrix4d();
5154
private Point3d[] originalCoords = null;
5255
private Point3d[] transformedCoords = null;
53-
private Set<List<Integer>> hashCodes = new HashSet<List<Integer>>();
56+
// Cache whether a permutation is invalid (null) vs has been added to rotations
57+
private Map<List<Integer>,Rotation> evaluatedPermutations = new HashMap<>();
5458

5559
private RotationGroup rotations = new RotationGroup();
5660

@@ -96,8 +100,13 @@ private void solve() {
96100
List<Double> angles = getAngles();
97101

98102
for (int i = 0; i < sphereCount; i++) {
103+
// Sampled orientation
104+
//TODO The SphereSampler samples 4D orientation space. We really
105+
// only need to sample 3D unit vectors, since we use limited
106+
// angles. -SB
99107
SphereSampler.getAxisAngle(i, sphereAngle);
100108

109+
// Each valid rotation angle
101110
for (double angle : angles) {
102111
// apply rotation
103112
sphereAngle.angle = angle;
@@ -113,14 +122,14 @@ private void solve() {
113122
List<Integer> permutation = getPermutation();
114123
// System.out.println("Rotation Solver: permutation: " + i + ": " + permutation);
115124

116-
boolean isValidPermuation = isValidPermutation(permutation);
117-
if (! isValidPermuation) {
118-
continue;
125+
// check if novel
126+
if ( evaluatedPermutations.containsKey(permutation)) {
127+
continue; //either invalid or already added
119128
}
120129

121-
boolean newPermutation = evaluatePermutation(permutation);
122-
if (newPermutation) {
123-
completeRotationGroup();
130+
Rotation newPermutation = isValidPermutation(permutation);
131+
if (newPermutation != null) {
132+
completeRotationGroup(newPermutation);
124133
}
125134

126135
// check if all symmetry operations have been found.
@@ -131,33 +140,84 @@ private void solve() {
131140
}
132141
}
133142

134-
private void completeRotationGroup() {
143+
/**
144+
* Combine current rotations to make all possible permutations.
145+
* If these are all valid, add them to the rotations
146+
* @param additionalRots Additional rotations we are considering adding to this.rotations
147+
* @return whether the rotations were valid and added
148+
*/
149+
private boolean completeRotationGroup(Rotation... additionalRots) {
135150
PermutationGroup g = new PermutationGroup();
136-
for (int i = 0; i < rotations.getOrder(); i++) {
137-
Rotation s = rotations.getRotation(i);
151+
for (Rotation s : rotations) {
138152
g.addPermutation(s.getPermutation());
139153
}
154+
for( Rotation s : additionalRots) {
155+
g.addPermutation(s.getPermutation());
156+
// inputs should not have been added already
157+
assert evaluatedPermutations.get(s.getPermutation()) == null;
158+
}
140159
g.completeGroup();
141160

142161
// the group is complete, nothing to do
143-
if (g.getOrder() == rotations.getOrder()) {
144-
return;
162+
if (g.getOrder() == rotations.getOrder()+additionalRots.length) {
163+
for( Rotation s : additionalRots) {
164+
addRotation(s);
165+
}
166+
return true;
145167
}
146168

147169
// try to complete the group
148-
for (int i = 0; i < g.getOrder(); i++) {
149-
List<Integer> permutation = g.getPermutation(i);
150-
151-
boolean isValidPermutation = isValidPermutation(permutation);
152-
if (isValidPermutation) {
170+
List<Rotation> newRots = new ArrayList<>(g.getOrder());
171+
// First, quick check for whether they're allowed
172+
for (List<Integer> permutation : g) {
173+
if( evaluatedPermutations.containsKey(permutation)) {
174+
Rotation rot = evaluatedPermutations.get(permutation);
175+
if( rot == null ) {
176+
return false;
177+
}
178+
} else {
179+
if( ! isAllowedPermutation(permutation)) {
180+
return false;
181+
}
182+
}
183+
}
184+
// Slower check including the superpositions
185+
for (List<Integer> permutation : g) {
186+
Rotation rot;
187+
if( evaluatedPermutations.containsKey(permutation)) {
188+
rot = evaluatedPermutations.get(permutation);
189+
} else {
190+
rot = isValidPermutation(permutation);
191+
}
153192

154-
// perform permutation of subunits
155-
evaluatePermutation(permutation);
193+
if( rot == null ) {
194+
// if any induced rotation is invalid, abort
195+
return false;
196+
}
197+
if(!evaluatedPermutations.containsKey(permutation)){ //novel
198+
newRots.add(rot);
156199
}
157200
}
201+
// Add rotations
202+
for( Rotation rot : newRots) {
203+
addRotation(rot);
204+
}
205+
return true;
206+
}
207+
208+
private void addRotation(Rotation rot) {
209+
evaluatedPermutations.put(rot.getPermutation(),rot);
210+
rotations.addRotation(rot);
158211
}
159212

160-
private boolean evaluatePermutation(List<Integer> permutation) {
213+
/**
214+
* Superimpose subunits based on the given permutation. Then check whether
215+
* the superposition passes RMSD thresholds and create a Rotation to
216+
* represent it if so.
217+
* @param permutation A list specifying which subunits should be aligned by the current transformation
218+
* @return A Rotation representing the permutation, or null if the superposition did not meet thresholds.
219+
*/
220+
private Rotation superimposePermutation(List<Integer> permutation) {
161221
// permutate subunits
162222
for (int j = 0, n = subunits.getSubunitCount(); j < n; j++) {
163223
transformedCoords[j].set(originalCoords[permutation.get(j)]);
@@ -176,17 +236,20 @@ private boolean evaluatePermutation(List<Integer> permutation) {
176236
// evaluate superposition of CA traces
177237
QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
178238
if (scores.getRmsd() < 0.0 || scores.getRmsd() > parameters.getRmsdThreshold()) {
179-
return false;
239+
return null;
180240
}
181241

182242
scores.setRmsdCenters(subunitRmsd);
183243
Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
184-
rotations.addRotation(symmetryOperation);
185-
return true;
244+
return symmetryOperation;
186245
}
187-
return false;
246+
return null;
188247
}
189248

249+
/**
250+
* Get valid rotation angles given the number of subunits
251+
* @return The rotation angle corresponding to each fold of {@link Subunits#getFolds()}
252+
*/
190253
private List<Double> getAngles() {
191254
int n = subunits.getSubunitCount();
192255
// for spherical symmetric cases, n cannot be higher than 60
@@ -210,23 +273,31 @@ private boolean isSpherical() {
210273
return m.getSymmetryClass(0.05) == MomentsOfInertia.SymmetryClass.SYMMETRIC;
211274
}
212275

213-
private boolean isValidPermutation(List<Integer> permutation) {
214-
if (permutation.size() == 0) {
215-
return false;
216-
}
276+
/**
277+
* Checks if a particular permutation is allowed and superimposes well.
278+
* Caches results.
279+
* @param permutation
280+
* @return null if invalid, or a rotation if valid
281+
*/
282+
private Rotation isValidPermutation(List<Integer> permutation) {
283+
if (permutation.size() == 0) {
284+
return null;
285+
}
217286

218-
// if this permutation is a duplicate, return false
219-
if (hashCodes.contains(permutation)) {
220-
return false;
287+
// cached value
288+
if (evaluatedPermutations.containsKey(permutation)) {
289+
return evaluatedPermutations.get(permutation);
221290
}
222291

223292
// check if permutation is allowed
224293
if (! isAllowedPermutation(permutation)) {
225-
return false;
294+
evaluatedPermutations.put(permutation, null);
295+
return null;
226296
}
227297

228-
// if this permutation is a duplicate, returns false
229-
return hashCodes.add(permutation);
298+
// check if superimposes
299+
Rotation rot = superimposePermutation(permutation);
300+
return rot;
230301
}
231302

232303
/**
@@ -237,13 +308,18 @@ private boolean isValidPermutation(List<Integer> permutation) {
237308
*/
238309
private boolean isAllowedPermutation(List<Integer> permutation) {
239310
List<Integer> seqClusterId = subunits.getSequenceClusterIds();
311+
int selfaligned = 0;
240312
for (int i = 0; i < permutation.size(); i++) {
241313
int j = permutation.get(i);
242-
if (i == j || seqClusterId.get(i) != seqClusterId.get(j)) {
314+
if ( seqClusterId.get(i) != seqClusterId.get(j)) {
243315
return false;
244316
}
317+
if(i == j ) {
318+
selfaligned++;
319+
}
245320
}
246-
return true;
321+
// either identity (all self aligned) or all unique
322+
return selfaligned == 0 || selfaligned == permutation.size();
247323
}
248324
/**
249325
* Adds translational component to rotation matrix
@@ -256,7 +332,7 @@ private void combineWithTranslation(Matrix4d rotation) {
256332
rotation.mul(rotation, centroidInverse);
257333
}
258334

259-
private Rotation createSymmetryOperation(List<Integer> permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) {
335+
private static Rotation createSymmetryOperation(List<Integer> permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) {
260336
Rotation s = new Rotation();
261337
s.setPermutation(new ArrayList<Integer>(permutation));
262338
s.setTransformation(new Matrix4d(transformation));
@@ -295,6 +371,11 @@ private double calcDistanceThreshold() {
295371
return distanceThreshold;
296372
}
297373

374+
/**
375+
* Compare this.transformedCoords with the original coords. For each
376+
* subunit, return the transformed subunit with the closest position.
377+
* @return A list mapping each subunit to the closest transformed subunit
378+
*/
298379
private List<Integer> getPermutation() {
299380
List<Integer> permutation = new ArrayList<Integer>(transformedCoords.length);
300381
double sum = 0.0f;

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

Lines changed: 11 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -24,9 +24,12 @@
2424
import javax.vecmath.Quat4d;
2525

2626
/**
27-
*
27+
* Represents an even coverage of quaternion space by 60 points. This grid is
28+
* defined by the vertices of one half of a hexacosichoron (600-cell). It
29+
* approximates all possible orientations to within 44.48 degrees.
2830
* @author Peter
2931
*/
32+
// This would be better named HexacosichoronSampler, since it's sampling 4D space. -SB
3033
public final class IcosahedralSampler {
3134
private static Quat4d quat = new Quat4d();
3235

@@ -39,7 +42,7 @@ public static int getSphereCount() {
3942
}
4043

4144
public static Quat4d getQuat4d(int index) {
42-
Quat4d q = new Quat4d(orientations[index]);
45+
Quat4d q = new Quat4d(orientations[index]); //ignores 5th element
4346
return q;
4447
}
4548

@@ -48,10 +51,12 @@ public static void getAxisAngle(int index, AxisAngle4d axisAngle) {
4851
axisAngle.set(quat);
4952
}
5053

51-
// # Orientation set c600v, number = 60, radius = 44.48 degrees
52-
// # $Id: c600v.quat 6102 2006-02-21 19:45:40Z ckarney $
53-
// # For more information, eee http://charles.karney.info/orientation/
54-
// format quaternion
54+
// # Orientation set c600v, number = 60, radius = 44.48 degrees
55+
// # $Id: c600v.quat 6102 2006-02-21 19:45:40Z ckarney $
56+
// # For more information, eee http://charles.karney.info/orientation/
57+
// format quaternion
58+
// The fifth column gives a weighting factor. Since the 600-cell is regular, all
59+
// orientations cover an equal fraction of orientation space and have equal weight.
5560
private static double[][] orientations = {
5661
{1.000000000f, 0.000000000f, 0.000000000f, 0.000000000f, 1.000000f},
5762
{0.000000000f, 1.000000000f, 0.000000000f, 0.000000000f, 1.000000f},

0 commit comments

Comments
 (0)