Avoid ASA calculation for redundant NCS interfaces (#823)

josemduarte · web-flow · commit ebabdf7b6439 · 2019-01-24T16:47:50.000-08:00
* WIP: towards calculating interface ASA non-redundantly for NCS structures

* Working implementation

* More assertions

* Docs

* More efficient

* Not storing the ASA values in all redundant interfaces, it can take a lot of memory

* Bugfix: removing also ncs-clustered interfaces in removeInterfacesBelowArea

* Bugfix: now removing interfaces below area properly. Test for it

* Bugfix: clustering of interfaces depended on the ids. It should work whether ids are there or not

* Better solution: propagate the interface area values

* Sorting with lambda and some docs

* Set groupAsas as a reference to the same objects, thus using no more mem

* No wildcard imports

* Bugfix: order of molecules in interfaces can be inverted with respect to reference

* Extracting constant

* A much better solution without a regex
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/GroupAsa.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/asa/GroupAsa.java
@@ -97,8 +97,8 @@ public GroupAsa(Group g) {
 		this.g = g;
 
 		int groupNoHSize = getGroupNoHSize();
-		atomAsaUs = new ArrayList<Double>(groupNoHSize);
-		atomAsaCs = new ArrayList<Double>(groupNoHSize);
+		atomAsaUs = new ArrayList<>(groupNoHSize);
+		atomAsaCs = new ArrayList<>(groupNoHSize);
 	}
 
 	private int getGroupNoHSize() {
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterface.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterface.java
@@ -104,18 +104,18 @@ public StructureInterface(
 			AtomContactSet contacts,
 			CrystalTransform firstTransf, CrystalTransform secondTransf) {
 
-		this.molecules = new Pair<Atom[]>(firstMolecule, secondMolecule);
-		this.moleculeIds = new Pair<String>(firstMoleculeId,secondMoleculeId);
+		this.molecules = new Pair<>(firstMolecule, secondMolecule);
+		this.moleculeIds = new Pair<>(firstMoleculeId,secondMoleculeId);
 		this.contacts = contacts;
-		this.transforms = new Pair<CrystalTransform>(firstTransf, secondTransf);
+		this.transforms = new Pair<>(firstTransf, secondTransf);
 	}
 
 	/**
 	 * Constructs an empty StructureInterface
 	 */
 	public StructureInterface() {
-		this.groupAsas1 = new TreeMap<ResidueNumber, GroupAsa>();
-		this.groupAsas2 = new TreeMap<ResidueNumber, GroupAsa>();
+		this.groupAsas1 = new TreeMap<>();
+		this.groupAsas2 = new TreeMap<>();
 	}
 
 	public int getId() {
@@ -133,7 +133,7 @@ public void setId(int id) {
 	 * @return
 	 */
 	public Pair<String> getCrystalIds() {
-		return new Pair<String>(
+		return new Pair<>(
 			moleculeIds.getFirst()+transforms.getFirst().getTransformId()+transforms.getFirst().getCrystalTranslation(),
 			moleculeIds.getSecond()+transforms.getSecond().getTransformId()+transforms.getSecond().getCrystalTranslation());
 	}
@@ -197,7 +197,16 @@ public void setTransforms(Pair<CrystalTransform> transforms) {
 		this.transforms = transforms;
 	}
 
-	protected void setAsas(double[] asas1, double[] asas2, int nSpherePoints, int nThreads, int cofactorSizeToUse) {
+	/**
+	 * Set ASA annotations by passing the uncomplexed ASA values of the 2 partners.
+	 * This will calculate complexed ASA and set the ASA values in the member variables.
+	 * @param asas1 ASA values for atoms of partner 1
+	 * @param asas2 ASA values for atoms of partner 2
+	 * @param nSpherePoints the number of sphere points to be used for complexed ASA calculation
+	 * @param nThreads the number of threads to be used for complexed ASA calculation
+	 * @param cofactorSizeToUse the minimum size of cofactor molecule (non-chain HET atoms) that will be used in ASA calculation
+	 */
+	void setAsas(double[] asas1, double[] asas2, int nSpherePoints, int nThreads, int cofactorSizeToUse) {
 
 		Atom[] atoms = getAtomsForAsa(cofactorSizeToUse);
 		AsaCalculator asaCalc = new AsaCalculator(atoms,
@@ -296,7 +305,7 @@ protected Atom[] getAtomsForAsa(int cofactorSizeToUse) {
 	 * @return
 	 */
 	private static final Atom[] getAllNonHAtomArray(Atom[] m, int minSizeHetAtomToInclude) {
-		List<Atom> atoms = new ArrayList<Atom>();
+		List<Atom> atoms = new ArrayList<>();
 
 		for (Atom a:m){
 
@@ -411,6 +420,14 @@ public void setFirstGroupAsa(GroupAsa groupAsa) {
 		groupAsas1.put(groupAsa.getGroup().getResidueNumber(), groupAsa);
 	}
 
+	void setFirstGroupAsas(Map<ResidueNumber, GroupAsa> firstGroupAsas) {
+		this.groupAsas1 = firstGroupAsas;
+	}
+
+	void setSecondGroupAsas(Map<ResidueNumber, GroupAsa> secondGroupAsas) {
+		this.groupAsas2 = secondGroupAsas;
+	}
+
 	/**
 	 * Gets a map of ResidueNumbers to GroupAsas for all groups of second chain.
 	 * @return
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java
@@ -21,7 +21,14 @@
 package org.biojava.nbio.structure.contact;
 
 import java.io.Serializable;
-import java.util.*;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Iterator;
+import java.util.List;
+import java.util.Map;
+import java.util.Optional;
+import java.util.Set;
+import java.util.TreeMap;
 
 import org.biojava.nbio.core.util.SingleLinkageClusterer;
 import org.biojava.nbio.structure.Atom;
@@ -72,6 +79,8 @@ public class StructureInterfaceList implements Serializable, Iterable<StructureI
 	private List<StructureInterfaceCluster> clusters = null;
 	private List<StructureInterfaceCluster> clustersNcs = null;
 
+	private Map<String, String> chainOrigNamesMap;
+
 	public StructureInterfaceList() {
 		this.list = new ArrayList<>();
 	}
@@ -129,8 +138,20 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) {
 		Map<String, Atom[]> uniqAsaChains = new TreeMap<>();
 		Map<String, double[]> chainAsas = new TreeMap<>();
 
+		List<StructureInterface> redundancyReducedList;
+		if (clustersNcs != null) {
+			redundancyReducedList = new ArrayList<>();
+			for (StructureInterfaceCluster ncsCluster : clustersNcs) {
+				// we use the first one in list as the only one for which we calculate ASAs
+				redundancyReducedList.add(ncsCluster.getMembers().get(0));
+			}
+
+		} else {
+			redundancyReducedList = list;
+		}
+
 		// first we gather rotation-unique chains (in terms of AU id and transform id)
-		for (StructureInterface interf:list) {
+		for (StructureInterface interf:redundancyReducedList) {
 			String molecId1 = interf.getMoleculeIds().getFirst()+interf.getTransforms().getFirst().getTransformId();
 			String molecId2 = interf.getMoleculeIds().getSecond()+interf.getTransforms().getSecond().getTransformId();
 
@@ -159,12 +180,12 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) {
 
 		logger.debug("Calculated uncomplexed ASA for {} orientation-unique chains. Time: {} s", uniqAsaChains.size(), ((end-start)/1000.0));
 
-		logger.debug ("Will calculate complexed ASA for {} pairwise complexes.", list.size());
+		logger.debug ("Will calculate complexed ASA for {} pairwise complexes.", redundancyReducedList.size());
 
 		start = System.currentTimeMillis();
 
 		// now we calculate the ASAs for each of the complexes
-		for (StructureInterface interf:list) {
+		for (StructureInterface interf:redundancyReducedList) {
 
 			String molecId1 = interf.getMoleculeIds().getFirst()+interf.getTransforms().getFirst().getTransformId();
 			String molecId2 = interf.getMoleculeIds().getSecond()+interf.getTransforms().getSecond().getTransformId();
@@ -176,15 +197,53 @@ public void calcAsas(int nSpherePoints, int nThreads, int cofactorSizeToUse) {
 		}
 		end = System.currentTimeMillis();
 
-		logger.debug("Calculated complexes ASA for {} pairwise complexes. Time: {} s", list.size(), ((end-start)/1000.0));
+		logger.debug("Calculated complexes ASA for {} pairwise complexes. Time: {} s", redundancyReducedList.size(), ((end-start)/1000.0));
 
+		// now let's populate the interface area value for the NCS-redundant ones from the reference interface (first one in list)
+		if (clustersNcs!=null) {
+			if (chainOrigNamesMap==null)  {
+				logger.warn("No chainOrigNamesMap is set. Considering NCS interfaces in same order as reference. This is likely a bug.");
+			}
+			for (StructureInterfaceCluster ncsCluster : clustersNcs) {
+				StructureInterface refInterf = ncsCluster.getMembers().get(0);
+				String refMolecId1 = refInterf.getMoleculeIds().getFirst();
+				for (int i=1;i<ncsCluster.getMembers().size();i++) {
+					StructureInterface member = ncsCluster.getMembers().get(i);
+					member.setTotalArea(refInterf.getTotalArea());
+					String molecId1 = member.getMoleculeIds().getFirst();
+					if (areMolecIdsSameOrder(refMolecId1, molecId1)) {
+						// we add the reference interface GroupAsas as the GroupAsas for all other members, like that
+						// ResidueNumbers won't match in their chain ids, but otherwise all info is there without using a lot of memory
+						member.setFirstGroupAsas(refInterf.getFirstGroupAsas());
+						member.setSecondGroupAsas(refInterf.getSecondGroupAsas());
+					} else {
+						member.setFirstGroupAsas(refInterf.getSecondGroupAsas());
+						member.setSecondGroupAsas(refInterf.getFirstGroupAsas());
+					}
+				}
+			}
+		}
 
 		// finally we sort based on the ChainInterface.comparable() (based in interfaceArea)
 		sort();
 	}
 
+	private boolean areMolecIdsSameOrder(String refMolecId, String molecId) {
+
+		if (chainOrigNamesMap==null) {
+			// we've already warned above
+			return true;
+		}
+
+		String refMolecIdOrig = chainOrigNamesMap.get(refMolecId);
+		String molecIdOrig = chainOrigNamesMap.get(molecId);
+
+		return (refMolecIdOrig.equals(molecIdOrig));
+	}
+
 	/**
 	 * Sorts the interface list and reassigns ids based on new sorting
+	 *
 	 */
 	public void sort() {
 		Collections.sort(list);
@@ -251,9 +310,21 @@ public void addNcsEquivalent(StructureInterface interfaceNew, StructureInterface
 		clustersNcs.add(newCluster);
 	}
 
+	/**
+	 * Sets a map with mapping from NCS chain names to original chain names.
+	 * Necessary when {@link #addNcsEquivalent(StructureInterface, StructureInterface)} is used and NCS equivalent
+	 * interfaces exist in this list and their names need mapping when setting ASAs.
+	 * @param chainOrigNamesMap a map of NCS chain name to original chain name
+	 */
+	public void setChainOrigNamesMap(Map<String, String> chainOrigNamesMap) {
+		this.chainOrigNamesMap = chainOrigNamesMap;
+	}
+
 	/**
 	 * Removes from this interface list all interfaces with areas
-	 * below the default cutoff area
+	 * below the default cutoff area.
+     * Note that this must be called after {@link #calcAsas(int, int, int)}, otherwise all areas would
+     * be 0 and thus all removed.
 	 * @see #DEFAULT_MINIMUM_INTERFACE_AREA
 	 */
 	public void removeInterfacesBelowArea() {
@@ -262,17 +333,18 @@ public void removeInterfacesBelowArea() {
 
 	/**
 	 * Removes from this interface list all interfaces with areas
-	 * below the given cutoff area
-	 * @param area
+	 * below the given cutoff area.
+     * Note that this must be called after {@link #calcAsas(int, int, int)}, otherwise all areas would
+     * be 0 and thus all removed.
+	 * @param area the minimum interface buried surface area to keep. Interfaces below this value will be removed.
 	 */
 	public void removeInterfacesBelowArea(double area) {
-		Iterator<StructureInterface> it = iterator();
-		while (it.hasNext()) {
-			StructureInterface interf = it.next();
-			if (interf.getTotalArea()<area) {
-				it.remove();
-			}
-		}
+
+		list.removeIf(interf -> interf.getTotalArea() < area);
+
+	    if (clustersNcs != null) {
+	    	clustersNcs.removeIf(ncsCluster -> ncsCluster.getMembers().get(0).getTotalArea() < area);
+        }
 	}
 
 	/**
@@ -291,6 +363,7 @@ public List<StructureInterfaceCluster> getClusters() {
 	 * Calculate the interface clusters for this StructureInterfaceList
 	 * using a contact overlap score to measure the similarity of interfaces.
 	 * Subsequent calls will use the cached value without recomputing the clusters.
+	 * The clusters will be assigned ids by sorting descending by {@link StructureInterfaceCluster#getTotalArea()}
 	 * @param contactOverlapScoreClusterCutoff the contact overlap score above which a pair will be
 	 * clustered
 	 * @return
@@ -300,7 +373,7 @@ public List<StructureInterfaceCluster> getClusters(double contactOverlapScoreClu
 			return clusters;
 		}
 
-		clusters = new ArrayList<StructureInterfaceCluster>();
+		clusters = new ArrayList<>();
 
 		// nothing to do if we have no interfaces
 		if (list.size()==0) return clusters;
@@ -324,9 +397,9 @@ public List<StructureInterfaceCluster> getClusters(double contactOverlapScoreClu
 
 		SingleLinkageClusterer slc = new SingleLinkageClusterer(matrix, true);
 
-		Map<Integer,Set<Integer>> clusteredIndices = slc.getClusters(contactOverlapScoreClusterCutoff);
+		Map<Integer, Set<Integer>> clusteredIndices = slc.getClusters(contactOverlapScoreClusterCutoff);
 		for (int clusterIdx:clusteredIndices.keySet()) {
-			List<StructureInterface> members = new ArrayList<StructureInterface>();
+			List<StructureInterface> members = new ArrayList<>();
 			for (int idx:clusteredIndices.get(clusterIdx)) {
 				members.add(list.get(idx));
 			}
@@ -335,8 +408,9 @@ public List<StructureInterfaceCluster> getClusters(double contactOverlapScoreClu
 			double averageScore = 0.0;
 			int countPairs = 0;
 			for (int i=0;i<members.size();i++) {
+                int iIdx = list.indexOf(members.get(i));
 				for (int j=i+1;j<members.size();j++) {
-					averageScore += matrix[members.get(i).getId()-1][members.get(j).getId()-1];
+					averageScore += matrix[iIdx][list.indexOf(members.get(j))];
 					countPairs++;
 				}
 			}
@@ -358,19 +432,15 @@ public List<StructureInterfaceCluster> getClusters(double contactOverlapScoreClu
 		}
 
 		// now we sort by areas (descending) and assign ids based on that sorting
-		Collections.sort(clusters, new Comparator<StructureInterfaceCluster>() {
-			@Override
-			public int compare(StructureInterfaceCluster o1, StructureInterfaceCluster o2) {
-				return Double.compare(o2.getTotalArea(), o1.getTotalArea()); //note we invert so that sorting is descending
-			}
+		clusters.sort((o1, o2) -> {
+			return Double.compare(o2.getTotalArea(), o1.getTotalArea()); //note we invert so that sorting is descending
 		});
 		int id = 1;
 		for (StructureInterfaceCluster cluster:clusters) {
 			cluster.setId(id);
 			id++;
 		}
 
-
 		return clusters;
 	}
 
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/xtal/CrystalBuilder.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/xtal/CrystalBuilder.java
@@ -44,6 +44,8 @@
 
 public class CrystalBuilder {
 
+	public static final String NCS_CHAINID_SUFFIX_CHAR = "n";
+
 	// Default number of cell neighbors to try in interface search (in 3 directions of space).
 	// In the search, only bounding box overlaps are tried, thus there's not so much overhead in adding
 	// more cells. We actually tested it and using numCells from 1 to 10 didn't change runtimes at all.
@@ -219,7 +221,10 @@ public StructureInterfaceList getUniqueInterfaces(double cutoff) {
 			return set;
 		}
 
-
+		// pass the chainOrigNames map in NCS case so that StructureInterfaceList can deal with original to NCS chain names conversion
+		if (chainOrigNames!=null) {
+			set.setChainOrigNamesMap(chainOrigNames);
+		}
 
 		// initialising the visited ArrayList for keeping track of symmetry redundancy
 		initialiseVisited();
@@ -319,7 +324,7 @@ private void calcInterfacesCrystal(StructureInterfaceList set, double cutoff) {
 
 						// 3) an operator can be "self redundant" if it is the inverse of itself (involutory, e.g. all pure 2-folds with no translation)
 						if (tt.isEquivalent(tt)) {
-							logger.debug("Transform "+tt+" is equivalent to itself, will skip half of i-chains to j-chains comparisons");
+							logger.debug("Transform {} is equivalent to itself, will skip half of i-chains to j-chains comparisons", tt.toString());
 							// in this case we can't skip the operator, but we can skip half of the matrix comparisons e.g. j>i
 							// we set a flag and do that within the loop below
 							selfEquivalent = true;
@@ -552,7 +557,7 @@ public void translate(Matrix4d m, Vector3d translation) {
 
 	/**
 	 * Apply the NCS operators in the given Structure adding new chains as needed.
-	 * All chains are (re)assigned ids of the form: original_chain_id+ncs_operator_index+"n".
+	 * All chains are (re)assigned ids of the form: original_chain_id+ncs_operator_index+{@value #NCS_CHAINID_SUFFIX_CHAR}.
 	 * @param structure
 	 *          the structure to expand
 	 * @param chainOrigNames
@@ -583,8 +588,8 @@ public static void expandNcsOps(Structure structure, Map<String,String> chainOri
 				Matrix4d m = ncsOps[iOperator];
 
 				Chain clonedChain = (Chain)c.clone();
-				String newChainId = cOrigId+(iOperator+1)+"n";
-				String newChainName = cOrigName+(iOperator+1)+"n";
+				String newChainId = cOrigId+(iOperator+1)+NCS_CHAINID_SUFFIX_CHAR;
+				String newChainName = cOrigName+(iOperator+1)+NCS_CHAINID_SUFFIX_CHAR;
 				clonedChain.setId(newChainId);
 				clonedChain.setName(newChainName);
 
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/xtal/TestInterfaceClustering.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/xtal/TestInterfaceClustering.java

Original file line number	Diff line number	Diff line change
`@@ -97,8 +97,8 @@ public GroupAsa(Group g) {`
`97`	`97`	`this.g = g;`
`98`	`98`
`99`	`99`	`int groupNoHSize = getGroupNoHSize();`
`100`		`- atomAsaUs = new ArrayList<Double>(groupNoHSize);`
`101`		`- atomAsaCs = new ArrayList<Double>(groupNoHSize);`
	`100`	`+ atomAsaUs = new ArrayList<>(groupNoHSize);`
	`101`	`+ atomAsaCs = new ArrayList<>(groupNoHSize);`
`102`	`102`	`}`
`103`	`103`
`104`	`104`	`private int getGroupNoHSize() {`