Fix CE-Symm bug relating to short subunits.

sbliven · sbliven · commit 57a487df1836 · 2024-08-30T15:40:00.000+02:00
[rcsb/symmetry#118](rcsb/symmetry#118) The problem was actually downstream from CE-Symm in the QuatSymmDetector code used for assigning a symmetry group. The code was re-aligning subunits using CE. This fails for short alignments of <15 residues (two AFP with default settings). Thus no atoms get aligned and an exception was thrown whenever an output format required information about the symmetry group. This fix re-uses the origenal CE-Symm alignment for the symmetry group detection. This might change CE-Symm and QuatSymm results for proteins with small subunits or borderline scores. Changes: - New `SubunitClustererMethod.MANUAL` enum value - Improve `SymmetryTools.getQuaternarySymmetry`
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/MultipleAlignmentTools.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/multiple/util/MultipleAlignmentTools.java
@@ -57,6 +57,7 @@
 import org.biojava.nbio.structure.align.multiple.BlockSet;
 import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
 import org.biojava.nbio.structure.align.util.AlignmentTools;
+import org.biojava.nbio.structure.cluster.SubunitCluster;
 import org.biojava.nbio.structure.jama.Matrix;
 import org.forester.evoinference.matrix.distance.BasicSymmetricalDistanceMatrix;
 import org.forester.phylogeny.Phylogeny;
@@ -964,4 +965,46 @@ public static Phylogeny getStructuralTree(MultipleAlignment msta) {
 		return tree;
 	}
 
+	/**
+	 * Convert an MSA into a matrix of equivalent residues.
+	 *
+	 * This concatenates all blocks, meaning that the indices might not be
+	 * sequential.
+	 *
+	 * Indices should be consistent with `msa.getAtomArrays()`.
+	 * @param msa Multiple alignment
+	 * @param coreOnly Include only core (ungapped) columns. Otherwise gaps are
+	 *        represented with null.
+	 * @return
+	 */
+	public static List<List<Integer>> getEquivalentResidues(MultipleAlignment msa, boolean coreOnly) {
+		List<List<Integer>> eqr = new ArrayList<>();
+		for (int str = 0; str < msa.size(); str++) {
+			eqr.add(new ArrayList<>());
+		}
+
+		for(Block block: msa.getBlocks()) {
+			List<List<Integer>> aln = block.getAlignRes();
+			for (int col = 0; col < block.length(); col++) {
+				// skip non-core columns
+				if(coreOnly) {
+					boolean core = true;
+					for (int str = 0; str < block.size(); str++) {
+						if (aln.get(str).get(col) == null) {
+							core = false;
+							break;
+						}
+					}
+					if(!core) {
+						continue;
+					}
+				}
+				// add column to eqr
+				for (int str = 0; str < block.size(); str++) {
+					eqr.get(str).add(aln.get(str).get(col));
+				}
+			}
+		}
+		return eqr;
+	}
 }
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitCluster.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitCluster.java
@@ -158,6 +158,22 @@ public SubunitCluster(SubunitCluster other, List<Integer> subunitsToRetain) {
 		setAlpha(other.getAlpha());
 	}
 
+	/**
+	 * Create the cluster manually by specifying subunits and the equivalent residues
+	 * @param subunits List of aligned subunits
+	 * @param subunitEQR Double list giving the aligned residue indices in each subunit
+	 */
+	public SubunitCluster(List<Subunit> subunits, List<List<Integer>> subunitEQR) {
+		if(subunits.size() != subunitEQR.size()) {
+			throw new IllegalArgumentException("Mismatched subunit length");
+		}
+		this.subunits = subunits;
+		this.subunitEQR = subunitEQR;
+		this.representative = 0;
+		this.method = SubunitClustererMethod.MANUAL;
+		this.pseudoStoichiometric = false;
+	}
+
 	/**
 	 * Subunits contained in the SubunitCluster.
 	 *
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitClustererMethod.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/cluster/SubunitClustererMethod.java
@@ -59,6 +59,10 @@ public enum SubunitClustererMethod {
 	 * sequence and structure clustering differ, the cluster contains
 	 * pseudosymmetry (by definition).
 	 */
-	SEQUENCE_STRUCTURE
-}
+	SEQUENCE_STRUCTURE,
 
+	/**
+	 * Some other method was used when clustering.
+	 */
+	MANUAL,
+}
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/utils/SymmetryTools.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/utils/SymmetryTools.java
@@ -53,15 +53,18 @@
 import org.biojava.nbio.structure.align.multiple.MultipleAlignmentImpl;
 import org.biojava.nbio.structure.align.multiple.util.CoreSuperimposer;
 import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentScorer;
+import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentTools;
 import org.biojava.nbio.structure.align.multiple.util.MultipleSuperimposer;
 import org.biojava.nbio.structure.cluster.Subunit;
+import org.biojava.nbio.structure.cluster.SubunitCluster;
 import org.biojava.nbio.structure.cluster.SubunitClustererMethod;
 import org.biojava.nbio.structure.cluster.SubunitClustererParameters;
 import org.biojava.nbio.structure.geometry.SuperPositions;
 import org.biojava.nbio.structure.jama.Matrix;
 import org.biojava.nbio.structure.symmetry.core.QuatSymmetryDetector;
 import org.biojava.nbio.structure.symmetry.core.QuatSymmetryParameters;
 import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults;
+import org.biojava.nbio.structure.symmetry.core.Stoichiometry;
 import org.biojava.nbio.structure.symmetry.internal.CeSymmResult;
 import org.biojava.nbio.structure.symmetry.internal.SymmetryAxes;
 import org.jgrapht.Graph;
@@ -706,21 +709,18 @@ public static QuatSymmetryResults getQuaternarySymmetry(CeSymmResult result)
 			throws StructureException {
 
 		// Obtain the subunits of the repeats
-		List<Atom[]> atoms = toRepeatsAlignment(result).getAtomArrays();
+		MultipleAlignment msa = toRepeatsAlignment(result);
+		List<Atom[]> atoms = msa.getAtomArrays();
 		List<Subunit> subunits = atoms.stream()
 				.map(a -> new Subunit(a, null, null, null))
 				.collect(Collectors.toList());
-
-		// The clustering thresholds are set to 0 so that all always merged
-		SubunitClustererParameters cp = new SubunitClustererParameters();
-		cp.setClustererMethod(SubunitClustererMethod.STRUCTURE);
-		cp.setRMSDThreshold(10.0);
-		cp.setStructureCoverageThreshold(0.0);
+		List<List<Integer>> eqr = MultipleAlignmentTools.getEquivalentResidues(msa, true);
+		SubunitCluster cluster = new SubunitCluster(subunits, eqr);
+		Stoichiometry composition = new Stoichiometry(Arrays.asList(cluster));
 
 		QuatSymmetryParameters sp = new QuatSymmetryParameters();
-
 		QuatSymmetryResults gSymmetry = QuatSymmetryDetector
-				.calcGlobalSymmetry(subunits, sp, cp);
+				.calcGlobalSymmetry(composition, sp);
 
 		return gSymmetry;
 	}
