Skip to content

Commit 8c94928

Browse files
lafitalafita
committed
Complete the QsAlign algorithm
Introduce a Demo and a skeleton for a Test.
1 parent 9a27e7f commit 8c94928

13 files changed

Lines changed: 589 additions & 119 deletions

File tree

‎biojava-structure/src/main/java/demo/DemoQsAlign.java‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
package demo;
2+
3+
import java.io.IOException;
4+
5+
import org.biojava.nbio.structure.Structure;
6+
import org.biojava.nbio.structure.StructureException;
7+
import org.biojava.nbio.structure.StructureIO;
8+
import org.biojava.nbio.structure.align.quaternary.QsAlign;
9+
import org.biojava.nbio.structure.align.quaternary.QsAlignParameters;
10+
import org.biojava.nbio.structure.align.quaternary.QsAlignResult;
11+
import org.biojava.nbio.structure.align.util.AtomCache;
12+
import org.biojava.nbio.structure.cluster.SubunitClustererParameters;
13+
14+
/**
15+
* Demo on how to use programatically {@link QsAlign} for the alignment of
16+
* quaternary structures.
17+
* <p>
18+
* Small oligomers: proliferating cell nuclear antigens (1PLR, 3HI8),
19+
* photosynthetic reaction centers (2JIY, 1DXR) ]
20+
* <p>
21+
* Big oligomers: cytochrome bc1 complexes (1bcc, 1kb9, 1qcr), phycocyanin
22+
* (2VML, 2BV8), bacterial ribosome (1FJG, 4V54).
23+
*
24+
*
25+
* @author Aleix Lafita
26+
* @since 5.0.0
27+
*
28+
*/
29+
public class DemoQsAlign {
30+
31+
public static void main(String[] args) throws IOException,
32+
StructureException {
33+
34+
// Remove lines when bug with ChemComp in MMTF is fixed TODO
35+
AtomCache cache = new AtomCache();
36+
cache.setUseMmCif(true);
37+
StructureIO.setAtomCache(cache);
38+
39+
// Align two trimeric DNA clamps
40+
Structure s1 = StructureIO.getStructure("1bcc");
41+
Structure s2 = StructureIO.getStructure("1kb9");
42+
43+
// Select the parameters for clustering and alignment
44+
SubunitClustererParameters clusterParams = new SubunitClustererParameters();
45+
QsAlignParameters alignParams = new QsAlignParameters();
46+
47+
QsAlignResult result = QsAlign
48+
.align(s1, s2, clusterParams, alignParams);
49+
50+
System.out.println(result);
51+
52+
}
53+
}

‎biojava-structure/src/main/java/org/biojava/nbio/structure/align/quaternary/QsAlign.java‎

Lines changed: 213 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,35 @@
11
package org.biojava.nbio.structure.align.quaternary;
22

3+
import java.util.ArrayList;
4+
import java.util.Arrays;
5+
import java.util.HashMap;
36
import java.util.List;
7+
import java.util.Map;
8+
import java.util.Map.Entry;
9+
import javax.vecmath.Matrix4d;
410

11+
import org.biojava.nbio.core.exceptions.CompoundNotFoundException;
12+
import org.biojava.nbio.structure.Atom;
13+
import org.biojava.nbio.structure.Calc;
14+
import org.biojava.nbio.structure.SVDSuperimposer;
515
import org.biojava.nbio.structure.Structure;
16+
import org.biojava.nbio.structure.StructureException;
17+
import org.biojava.nbio.structure.align.multiple.Block;
18+
import org.biojava.nbio.structure.align.multiple.BlockImpl;
19+
import org.biojava.nbio.structure.align.multiple.BlockSet;
20+
import org.biojava.nbio.structure.align.multiple.BlockSetImpl;
21+
import org.biojava.nbio.structure.align.multiple.MultipleAlignment;
22+
import org.biojava.nbio.structure.align.multiple.MultipleAlignmentEnsembleImpl;
23+
import org.biojava.nbio.structure.align.multiple.MultipleAlignmentImpl;
24+
import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentScorer;
25+
import org.biojava.nbio.structure.align.multiple.util.ReferenceSuperimposer;
626
import org.biojava.nbio.structure.cluster.Subunit;
27+
import org.biojava.nbio.structure.cluster.SubunitCluster;
28+
import org.biojava.nbio.structure.cluster.SubunitClusterer;
729
import org.biojava.nbio.structure.cluster.SubunitClustererParameters;
830
import org.biojava.nbio.structure.cluster.SubunitExtractor;
31+
import org.slf4j.Logger;
32+
import org.slf4j.LoggerFactory;
933

1034
/**
1135
* Quaternary Structure Alignment (QS-Align). The algorithm takes as input two
@@ -19,16 +43,200 @@
1943
*/
2044
public class QsAlign {
2145

46+
private static final Logger logger = LoggerFactory.getLogger(QsAlign.class);
47+
2248
public static QsAlignResult align(Structure s1, Structure s2,
23-
SubunitClustererParameters cParams, QsAlignParameters aParams) {
49+
SubunitClustererParameters cParams, QsAlignParameters aParams)
50+
throws StructureException {
2451
return align(SubunitExtractor.extractSubunits(s1, cParams),
2552
SubunitExtractor.extractSubunits(s2, cParams), cParams, aParams);
2653
}
2754

2855
public static QsAlignResult align(List<Subunit> s1, List<Subunit> s2,
29-
SubunitClustererParameters cParams, QsAlignParameters aParams) {
30-
31-
return null;
32-
}
56+
SubunitClustererParameters cParams, QsAlignParameters aParams)
57+
throws StructureException {
58+
59+
QsAlignResult result = new QsAlignResult(s1, s2);
60+
61+
// SETP 1: cluster each group of subunits O(N^2*L^2) - intra
62+
List<SubunitCluster> c1 = SubunitClusterer.cluster(s1, cParams);
63+
List<SubunitCluster> c2 = SubunitClusterer.cluster(s2, cParams);
64+
65+
// STEP 2: match each subunit cluster between groups O(N^2*L^2) - inter
66+
Map<Integer, Integer> clusterMap = new HashMap<Integer, Integer>();
67+
for (int i = 0; i < c1.size(); i++) {
68+
for (int j = 0; j < c2.size(); j++) {
69+
70+
if (clusterMap.keySet().contains(i))
71+
break;
72+
if (clusterMap.values().contains(j))
73+
continue;
74+
75+
switch (cParams.getClustererMethod()) {
76+
77+
case IDENTITY:
78+
if (c1.get(i).mergeIdentical(c2.get(j)))
79+
clusterMap.put(i, j);
80+
break;
81+
82+
case SEQUENCE:
83+
try {
84+
if (c1.get(i).mergeSequence(c2.get(j),
85+
cParams.getSequenceIdentityThreshold(),
86+
cParams.getCoverageThreshold()))
87+
clusterMap.put(i, j);
88+
} catch (CompoundNotFoundException e) {
89+
logger.warn("Could compare by Sequence. {}",
90+
e.getMessage());
91+
}
92+
break;
93+
94+
default: // case STRUCTURE:
95+
if (c1.get(i).mergeStructure(c2.get(j),
96+
cParams.getRmsdThreshold(),
97+
cParams.getCoverageThreshold()))
98+
clusterMap.put(i, j);
99+
break;
100+
}
101+
}
102+
}
103+
104+
// STEP 3: Align the assemblies for each cluster match O(L^2+N^2)
105+
for (int globalKey : clusterMap.keySet()) {
106+
107+
// Obtain the clusters
108+
SubunitCluster clust1 = c1.get(globalKey);
109+
SubunitCluster clust2 = c2.get(clusterMap.get(globalKey));
110+
111+
// Take the cluster match as reference and obtain transformation
112+
int index1 = 0;
113+
int index2 = clust1.size() - clust2.size();
114+
115+
Atom[] atoms1 = clust1.getAlignedAtomsSubunit(index1);
116+
Atom[] atoms2 = clust1.getAlignedAtomsSubunit(index2);
117+
118+
SVDSuperimposer svd = new SVDSuperimposer(atoms1, atoms2);
119+
Matrix4d trans = svd.getTransformation();
120+
121+
// Map the subunits of each cluster to their spatial equivalents
122+
Map<Integer, Map<Integer, Integer>> clustSubunitMap = new HashMap<Integer, Map<Integer, Integer>>();
123+
124+
for (int key : clusterMap.keySet()) {
125+
126+
// Obtain the clusters
127+
clust1 = c1.get(key);
128+
clust2 = c2.get(clusterMap.get(key));
129+
130+
// Take the cluster match as reference and obtain transformation
131+
index1 = 0;
132+
index2 = clust1.size() - clust2.size();
133+
134+
// Map the subunits of the cluster to their spatial equivalents
135+
Map<Integer, Integer> subunitMap = new HashMap<Integer, Integer>();
136+
137+
for (int i = 0; i < index2; i++) {
138+
for (int j = index2; j < clust1.size(); j++) {
33139

140+
if (subunitMap.keySet().contains(i))
141+
break;
142+
if (subunitMap.values().contains(j))
143+
continue;
144+
145+
// Obtain centroids and transform the second
146+
Atom centr1 = Calc.getCentroid(clust1
147+
.getAlignedAtomsSubunit(i));
148+
Atom centr2 = Calc.getCentroid(clust1
149+
.getAlignedAtomsSubunit(j));
150+
Calc.transform(centr2, trans);
151+
152+
if (Calc.getDistance(centr1, centr2) < aParams
153+
.getdCutoff())
154+
subunitMap.put(i, j);
155+
}
156+
}
157+
158+
clustSubunitMap.put(key, subunitMap);
159+
}
160+
161+
// Unfold the nested map into subunit map and alignment
162+
Map<Integer, Integer> subunitMap = new HashMap<Integer, Integer>();
163+
List<Integer> alignRes1 = new ArrayList<Integer>();
164+
List<Integer> alignRes2 = new ArrayList<Integer>();
165+
List<Atom> atomArray1 = new ArrayList<Atom>();
166+
List<Atom> atomArray2 = new ArrayList<Atom>();
167+
168+
for (int key : clustSubunitMap.keySet()) {
169+
170+
// Obtain the cluster and the alignment in it
171+
SubunitCluster cluster = c1.get(key);
172+
List<List<Integer>> clusterEqrs = cluster
173+
.getMultipleAlignment().getBlock(0).getAlignRes();
174+
175+
for (Entry<Integer, Integer> pair : clustSubunitMap.get(key)
176+
.entrySet()) {
177+
178+
int i = pair.getKey();
179+
int j = pair.getValue();
180+
181+
// Obtain the indices of the original Subunit Lists
182+
int orig1 = s1.indexOf(cluster.getSubunits().get(i));
183+
int orig2 = s2.indexOf(cluster.getSubunits().get(j));
184+
185+
// Append rescaled aligned residue indices
186+
for (Integer eqr : clusterEqrs.get(i))
187+
alignRes1.add(eqr + atomArray1.size());
188+
for (Integer eqr : clusterEqrs.get(j))
189+
alignRes2.add(eqr + atomArray2.size());
190+
191+
// Apend atoms to the arrays
192+
atomArray1.addAll(Arrays.asList(s1.get(orig1)
193+
.getRepresentativeAtoms()));
194+
atomArray2.addAll(Arrays.asList(s2.get(orig2)
195+
.getRepresentativeAtoms()));
196+
197+
subunitMap.put(orig1, orig2);
198+
}
199+
}
200+
201+
// Evaluate the goodness of the match with an alignment object
202+
MultipleAlignment msa = new MultipleAlignmentImpl();
203+
msa.setEnsemble(new MultipleAlignmentEnsembleImpl());
204+
msa.getEnsemble().setAtomArrays(
205+
Arrays.asList(new Atom[][] {
206+
atomArray1.toArray(new Atom[atomArray1.size()]),
207+
atomArray2.toArray(new Atom[atomArray2.size()]) }));
208+
209+
// Fill in the alignment information
210+
BlockSet bs = new BlockSetImpl(msa);
211+
Block b = new BlockImpl(bs);
212+
List<List<Integer>> alignRes = new ArrayList<List<Integer>>(2);
213+
alignRes.add(alignRes1);
214+
alignRes.add(alignRes2);
215+
b.setAlignRes(alignRes);
216+
217+
// Fill in the transformation matrices
218+
new ReferenceSuperimposer().superimpose(msa);
219+
220+
// Calculate some scores
221+
MultipleAlignmentScorer.calculateScores(msa);
222+
223+
// If it is the best match found so far store it
224+
if (subunitMap.size() > result.getSubunitMap().size()) {
225+
result.setSubunitMap(subunitMap);
226+
result.setAlignment(msa);
227+
} else if (subunitMap.size() == result.getSubunitMap().size()) {
228+
if (result.getAlignment() == null) {
229+
result.setSubunitMap(subunitMap);
230+
result.setAlignment(msa);
231+
} else if (msa.getScore(MultipleAlignmentScorer.RMSD) < result
232+
.getRmsd()) {
233+
result.setSubunitMap(subunitMap);
234+
result.setAlignment(msa);
235+
}
236+
}
237+
238+
}
239+
240+
return result;
241+
}
34242
}

‎biojava-structure/src/main/java/org/biojava/nbio/structure/align/quaternary/QsAlignParameters.java‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,15 +12,21 @@ public class QsAlignParameters {
1212
private double dCutoff = 10.0;
1313

1414
/**
15-
* The maximum allowed distance between two aligned residues of equivalent
16-
* Subunits.
15+
* The maximum allowed distance between the centroids of two equivalent
16+
* Subunits, in A.
1717
*
18-
* @return
18+
* @return dCutoff
1919
*/
2020
public double getdCutoff() {
2121
return dCutoff;
2222
}
2323

24+
/**
25+
* The maximum allowed distance between the centroids of two equivalent
26+
* Subunits, in A.
27+
*
28+
* @param dCutoff
29+
*/
2430
public void setdCutoff(double dCutoff) {
2531
this.dCutoff = dCutoff;
2632
}

0 commit comments

Comments
 (0)