Skip to content

Commit 7c521ed

Browse files
josemduartejosemduarte
committed
Use a greedy approach (leader algorithm) for a much more scalable interface clustering
1 parent a9dd08f commit 7c521ed

2 files changed

Lines changed: 171 additions & 57 deletions

File tree

‎biojava-structure/src/main/java/org/biojava/nbio/structure/contact/StructureInterfaceList.java‎

Lines changed: 110 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -27,11 +27,12 @@
2727
import java.util.List;
2828
import java.util.Map;
2929
import java.util.Optional;
30-
import java.util.Set;
3130
import java.util.TreeMap;
31+
import java.util.function.Function;
3232

33-
import org.biojava.nbio.core.util.SingleLinkageClusterer;
3433
import org.biojava.nbio.structure.Atom;
34+
import org.biojava.nbio.structure.Chain;
35+
import org.biojava.nbio.structure.EntityInfo;
3536
import org.biojava.nbio.structure.Structure;
3637
import org.biojava.nbio.structure.asa.AsaCalculator;
3738
import org.biojava.nbio.structure.xtal.CrystalBuilder;
@@ -368,6 +369,10 @@ public List<StructureInterfaceCluster> getClusters() {
368369
* using Jaccard contact set scores to measure the similarity of interfaces.
369370
* Subsequent calls will use the cached value without recomputing the clusters.
370371
* The clusters will be assigned ids by sorting descending by {@link StructureInterfaceCluster#getTotalArea()}
372+
* <p>
373+
* Interfaces are compared by entity pairs (see {@link #ENTITY_ID_PAIR}) and grouped with the leader
374+
* algorithm described in {@link #clusterInterfaces(List, Function, double)}, thus the result depends on
375+
* the order of the list (by default descending by area, see {@link #sort()}).
371376
* @param contactOverlapScoreClusterCutoff the contact overlap score above which a pair will be
372377
* clustered
373378
* @return
@@ -377,77 +382,125 @@ public List<StructureInterfaceCluster> getClusters(double contactOverlapScoreClu
377382
return clusters;
378383
}
379384

380-
clusters = new ArrayList<>();
381-
382-
// nothing to do if we have no interfaces
383-
if (list.isEmpty()) return clusters;
384-
385-
logger.debug("Calculating all-vs-all Jaccard scores for {} interfaces", list.size());
386-
double[][] matrix = new double[list.size()][list.size()];
385+
List<StructureInterfaceCluster> singletons = new ArrayList<>(list.size());
386+
for (StructureInterface interf : list) {
387+
StructureInterfaceCluster cluster = new StructureInterfaceCluster();
388+
cluster.addMember(interf);
389+
singletons.add(cluster);
390+
}
387391

388-
for (int i=0;i<list.size();i++) {
389-
for (int j=i+1;j<list.size();j++) {
390-
StructureInterface iInterf = list.get(i);
391-
StructureInterface jInterf = list.get(j);
392+
clusters = clusterInterfaces(singletons, ENTITY_ID_PAIR, contactOverlapScoreClusterCutoff);
392393

393-
double scoreDirect = iInterf.getContactOverlapScore(jInterf, false);
394-
double scoreInvert = iInterf.getContactOverlapScore(jInterf, true);
394+
// now we sort by areas (descending) and assign ids based on that sorting
395+
clusters.sort((o1, o2) -> Double.compare(o2.getTotalArea(), o1.getTotalArea())); //note we invert so that sorting is descending
395396

396-
double maxScore = Math.max(scoreDirect, scoreInvert);
397+
int id = 1;
398+
for (StructureInterfaceCluster cluster:clusters) {
399+
cluster.setId(id);
400+
id++;
401+
}
397402

398-
matrix[i][j] = maxScore;
399-
}
403+
return clusters;
404+
}
400405

406+
/**
407+
* Identifies each side of an interface by the entity id ({@link EntityInfo#getMolId()}) of its parent chain.
408+
* Returns null if any of the parent chains or entities can't be found.
409+
* @see #clusterInterfaces(List, Function, double)
410+
*/
411+
public static final Function<StructureInterface, Pair<Integer>> ENTITY_ID_PAIR = interf -> {
412+
Pair<Chain> chains = interf.getParentChains();
413+
if (chains == null || chains.getFirst().getEntityInfo() == null || chains.getSecond().getEntityInfo() == null) {
414+
return null;
401415
}
416+
return new Pair<>(chains.getFirst().getEntityInfo().getMolId(), chains.getSecond().getEntityInfo().getMolId());
417+
};
402418

403-
logger.debug("Will now cluster {} interfaces based on full all-vs-all Jaccard scores matrix", list.size());
404-
SingleLinkageClusterer slc = new SingleLinkageClusterer(matrix, true);
419+
/**
420+
* Groups interface clusters further, merging any two whose representatives (first members) have a
421+
* contact overlap score above the given cutoff.
422+
* <p>
423+
* This is the leader algorithm: the representative of the first cluster is compared to the representatives of all
424+
* subsequent clusters, merging those above cutoff into it; then the same is repeated for the next remaining cluster.
425+
* It performs at most n(n-1)/2 score calculations, and about n*k for n input and k output clusters.
426+
* Unlike single linkage clustering, similarity is not transitive through non-representative members
427+
* and the result depends on the input order: placing larger interfaces first makes them the representatives.
428+
* <p>
429+
* Two interfaces are only compared if their identifier pairs match, in the same or inverted order.
430+
* The score is then {@link StructureInterface#getContactOverlapScore(StructureInterface, boolean)} non-inverted
431+
* or inverted respectively, or the maximum of both if all 4 identifiers are equal.
432+
* <p>
433+
* The input clusters are not modified. The output clusters are new objects whose members are the
434+
* representative's cluster members followed by the members of the clusters merged into it, and whose average
435+
* score is the average of the scores that caused merges (1.0 if no merges). The cluster back-references of all
436+
* members are set to the output clusters. Ids are not assigned.
437+
* @param initialClusters the starting clusters, e.g. one singleton cluster per interface. The first member of
438+
* each is its representative
439+
* @param idPairFunction gives the pair of identifiers for the 2 sides of an interface, e.g. {@link #ENTITY_ID_PAIR}.
440+
* If it returns null the interface is not merged with any other
441+
* @param contactOverlapScoreClusterCutoff the contact overlap score above which a pair of clusters is merged
442+
* @return the merged clusters, in the order of their representatives in the input
443+
* @since 7.3.0
444+
*/
445+
public static <T> List<StructureInterfaceCluster> clusterInterfaces(List<StructureInterfaceCluster> initialClusters,
446+
Function<StructureInterface, Pair<T>> idPairFunction,
447+
double contactOverlapScoreClusterCutoff) {
448+
int n = initialClusters.size();
449+
logger.debug("Clustering {} interface clusters by contact overlap score", n);
450+
451+
List<StructureInterface> reps = new ArrayList<>(n);
452+
List<Pair<T>> idPairs = new ArrayList<>(n);
453+
for (StructureInterfaceCluster cluster : initialClusters) {
454+
StructureInterface rep = cluster.getMembers().get(0);
455+
reps.add(rep);
456+
idPairs.add(idPairFunction.apply(rep));
457+
}
405458

406-
Map<Integer, Set<Integer>> clusteredIndices = slc.getClusters(contactOverlapScoreClusterCutoff);
407-
for (int clusterIdx:clusteredIndices.keySet()) {
408-
List<StructureInterface> members = new ArrayList<>();
409-
for (int idx:clusteredIndices.get(clusterIdx)) {
410-
members.add(list.get(idx));
411-
}
459+
List<StructureInterfaceCluster> result = new ArrayList<>();
460+
boolean[] merged = new boolean[n];
461+
for (int i = 0; i < n; i++) {
462+
if (merged[i]) continue;
412463
StructureInterfaceCluster cluster = new StructureInterfaceCluster();
413-
cluster.setMembers(members);
414-
double averageScore = 0.0;
415-
int countPairs = 0;
416-
for (int i=0;i<members.size();i++) {
417-
int iIdx = list.indexOf(members.get(i));
418-
for (int j=i+1;j<members.size();j++) {
419-
averageScore += matrix[iIdx][list.indexOf(members.get(j))];
420-
countPairs++;
464+
cluster.getMembers().addAll(initialClusters.get(i).getMembers());
465+
double sumScores = 0;
466+
int countMerges = 0;
467+
// descending, so that merged members are appended in the same order as in previous implementations
468+
for (int j = n - 1; j > i; j--) {
469+
if (merged[j]) continue;
470+
double score = getContactOverlapScore(reps.get(i), idPairs.get(i), reps.get(j), idPairs.get(j));
471+
if (score > contactOverlapScoreClusterCutoff) {
472+
cluster.getMembers().addAll(initialClusters.get(j).getMembers());
473+
merged[j] = true;
474+
sumScores += score;
475+
countMerges++;
421476
}
422477
}
423-
if (countPairs>0) {
424-
averageScore = averageScore/countPairs;
425-
} else {
426-
// if only one interface in cluster we set the score to the maximum
427-
averageScore = 1.0;
478+
cluster.setAverageScore(countMerges > 0 ? sumScores / countMerges : 1.0);
479+
for (StructureInterface member : cluster.getMembers()) {
480+
member.setCluster(cluster);
428481
}
429-
cluster.setAverageScore(averageScore);
430-
clusters.add(cluster);
482+
result.add(cluster);
431483
}
432484

433-
// finally we have to set the back-references in each StructureInterface
434-
for (StructureInterfaceCluster cluster:clusters) {
435-
for (StructureInterface interf:cluster.getMembers()) {
436-
interf.setCluster(cluster);
437-
}
438-
}
439-
logger.debug("Done clustering {} interfaces based on full all-vs-all Jaccard scores matrix. Found a total of {} clusters", list.size(), clusters.size());
440-
441-
// now we sort by areas (descending) and assign ids based on that sorting
442-
clusters.sort((o1, o2) -> Double.compare(o2.getTotalArea(), o1.getTotalArea())); //note we invert so that sorting is descending
485+
logger.debug("Found {} clusters from {} input clusters at contact overlap score cutoff {}", result.size(), n, contactOverlapScoreClusterCutoff);
486+
return result;
487+
}
443488

444-
int id = 1;
445-
for (StructureInterfaceCluster cluster:clusters) {
446-
cluster.setId(id);
447-
id++;
489+
private static <T> double getContactOverlapScore(StructureInterface iInterf, Pair<T> iIds, StructureInterface jInterf, Pair<T> jIds) {
490+
if (iIds == null || jIds == null) {
491+
return 0;
448492
}
449-
450-
return clusters;
493+
boolean direct = iIds.getFirst().equals(jIds.getFirst()) && iIds.getSecond().equals(jIds.getSecond());
494+
boolean inverted = iIds.getFirst().equals(jIds.getSecond()) && iIds.getSecond().equals(jIds.getFirst());
495+
if (direct && inverted) {
496+
// all 4 ids are the same: the order of the sides is not known, so we try both
497+
return Math.max(iInterf.getContactOverlapScore(jInterf, false), iInterf.getContactOverlapScore(jInterf, true));
498+
} else if (direct) {
499+
return iInterf.getContactOverlapScore(jInterf, false);
500+
} else if (inverted) {
501+
return iInterf.getContactOverlapScore(jInterf, true);
502+
}
503+
return 0;
451504
}
452505

453506
@Override

‎biojava-structure/src/test/java/org/biojava/nbio/structure/xtal/TestInterfaceClustering.java‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,14 +36,18 @@
3636
import org.biojava.nbio.structure.StructureIO;
3737
import org.biojava.nbio.structure.align.util.AtomCache;
3838
import org.biojava.nbio.structure.asa.GroupAsa;
39+
import org.biojava.nbio.structure.contact.InterfaceFinder;
40+
import org.biojava.nbio.structure.contact.Pair;
3941
import org.biojava.nbio.structure.contact.StructureInterface;
4042
import org.biojava.nbio.structure.contact.StructureInterfaceCluster;
4143
import org.biojava.nbio.structure.contact.StructureInterfaceList;
4244
import org.biojava.nbio.structure.io.FileParsingParameters;
4345
import org.biojava.nbio.structure.io.PDBFileParser;
46+
import org.biojava.nbio.structure.quaternary.BiologicalAssemblyBuilder;
4447
import org.junit.Test;
4548

4649
import javax.vecmath.Matrix4d;
50+
import java.util.function.Function;
4751

4852
public class TestInterfaceClustering {
4953

@@ -272,4 +276,61 @@ public void test3DDONoSeqresPdb() throws IOException, StructureException {
272276

273277

274278
}
279+
280+
/**
281+
* Two-step clustering of the interfaces of an assembly: first geometrically identical interfaces (same asym ids),
282+
* then similar interfaces between same entities, starting from the clusters of the first step.
283+
*/
284+
@Test
285+
public void testClusterInterfacesTwoSteps() throws IOException, StructureException {
286+
287+
AtomCache cache = new AtomCache();
288+
FileParsingParameters params = new FileParsingParameters();
289+
params.setAlignSeqRes(true);
290+
cache.setFileParsingParams(params);
291+
cache.setFiletype(StructureFiletype.CIF);
292+
StructureIO.setAtomCache(cache);
293+
294+
// 6cco assembly 1: D3 with 2 entities, 12 interfaces
295+
Structure assembly = StructureIO.getBiologicalAssembly("6cco", 1, false);
296+
List<StructureInterface> interfaces = new InterfaceFinder(assembly).getAllInterfaces().getList();
297+
assertEquals(12, interfaces.size());
298+
299+
Function<StructureInterface, Pair<String>> asymIdPair = interf -> new Pair<>(
300+
stripOperators(interf.getParentChains().getFirst().getId()),
301+
stripOperators(interf.getParentChains().getSecond().getId()));
302+
303+
List<StructureInterfaceCluster> singletons = new ArrayList<>();
304+
for (StructureInterface interf : interfaces) {
305+
StructureInterfaceCluster cluster = new StructureInterfaceCluster();
306+
cluster.addMember(interf);
307+
singletons.add(cluster);
308+
}
309+
310+
List<StructureInterfaceCluster> identicalClusters = StructureInterfaceList.clusterInterfaces(singletons, asymIdPair, 0.99);
311+
assertEquals(4, identicalClusters.size());
312+
for (StructureInterfaceCluster cluster : identicalClusters) {
313+
assertEquals(3, cluster.getMembers().size());
314+
cluster.getMembers().forEach(m -> assertSame(cluster, m.getCluster()));
315+
}
316+
// input is not modified
317+
singletons.forEach(c -> assertEquals(1, c.getMembers().size()));
318+
319+
List<StructureInterfaceCluster> entityClusters = StructureInterfaceList.clusterInterfaces(identicalClusters, StructureInterfaceList.ENTITY_ID_PAIR, StructureInterfaceList.DEFAULT_CONTACT_OVERLAP_SCORE_CLUSTER_CUTOFF);
320+
assertEquals(3, entityClusters.size());
321+
int total = 0;
322+
for (StructureInterfaceCluster cluster : entityClusters) {
323+
// representatives are preserved from the first step
324+
assertTrue(identicalClusters.stream().anyMatch(c -> c.getMembers().get(0) == cluster.getMembers().get(0)));
325+
cluster.getMembers().forEach(m -> assertSame(cluster, m.getCluster()));
326+
total += cluster.getMembers().size();
327+
}
328+
assertEquals(12, total);
329+
identicalClusters.forEach(c -> assertEquals(3, c.getMembers().size()));
330+
}
331+
332+
private static String stripOperators(String chainId) {
333+
return chainId.substring(0, chainId.indexOf(BiologicalAssemblyBuilder.SYM_CHAIN_ID_SEPARATOR));
334+
}
335+
275336
}

0 commit comments

Comments
 (0)