001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021package org.biojava.nbio.structure.align.quaternary; 022 023import java.util.ArrayList; 024import java.util.Arrays; 025import java.util.HashMap; 026import java.util.List; 027import java.util.Map; 028import java.util.Map.Entry; 029 030import javax.vecmath.Matrix4d; 031 032import org.biojava.nbio.structure.Atom; 033import org.biojava.nbio.structure.Calc; 034import org.biojava.nbio.structure.Structure; 035import org.biojava.nbio.structure.StructureException; 036import org.biojava.nbio.structure.StructureTools; 037import org.biojava.nbio.structure.align.multiple.Block; 038import org.biojava.nbio.structure.align.multiple.BlockImpl; 039import org.biojava.nbio.structure.align.multiple.BlockSet; 040import org.biojava.nbio.structure.align.multiple.BlockSetImpl; 041import org.biojava.nbio.structure.align.multiple.MultipleAlignment; 042import org.biojava.nbio.structure.align.multiple.MultipleAlignmentEnsembleImpl; 043import org.biojava.nbio.structure.align.multiple.MultipleAlignmentImpl; 044import org.biojava.nbio.structure.align.multiple.util.MultipleAlignmentScorer; 045import org.biojava.nbio.structure.align.multiple.util.ReferenceSuperimposer; 046import org.biojava.nbio.structure.cluster.Subunit; 047import org.biojava.nbio.structure.cluster.SubunitCluster; 048import org.biojava.nbio.structure.cluster.SubunitClusterer; 049import org.biojava.nbio.structure.cluster.SubunitClustererParameters; 050import org.biojava.nbio.structure.cluster.SubunitExtractor; 051import org.biojava.nbio.structure.contact.Pair; 052import org.biojava.nbio.structure.geometry.SuperPositions; 053import org.biojava.nbio.structure.geometry.UnitQuaternions; 054import org.slf4j.Logger; 055import org.slf4j.LoggerFactory; 056 057/** 058 * Quaternary Structure Alignment (QS-Align). The algorithm takes as input two 059 * protein structures at the quaternary structure level (multiple chains or 060 * subunits) and calculates the equivalent subunit matching and a residue-based 061 * alignment, together with usual alignment quality scores. 062 * 063 * @author Aleix Lafita 064 * @since 5.0.0 065 * 066 */ 067public class QsAlign { 068 069 private static final Logger logger = LoggerFactory.getLogger(QsAlign.class); 070 071 public static QsAlignResult align(Structure s1, Structure s2, 072 SubunitClustererParameters cParams, QsAlignParameters aParams) 073 throws StructureException { 074 return align( 075 SubunitExtractor.extractSubunits(s1, 076 cParams.getAbsoluteMinimumSequenceLength(), 077 cParams.getMinimumSequenceLengthFraction(), 078 cParams.getMinimumSequenceLength()), 079 SubunitExtractor.extractSubunits(s2, 080 cParams.getAbsoluteMinimumSequenceLength(), 081 cParams.getMinimumSequenceLengthFraction(), 082 cParams.getMinimumSequenceLength()), cParams, aParams); 083 } 084 085 public static QsAlignResult align(List<Subunit> s1, List<Subunit> s2, 086 SubunitClustererParameters cParams, QsAlignParameters aParams) 087 throws StructureException { 088 089 QsAlignResult result = new QsAlignResult(s1, s2); 090 091 // SETP 1: cluster each group of subunits O(N^2*L^2) - intra 092 093 List<SubunitCluster> c1 = SubunitClusterer.cluster(s1, cParams).getClusters(); 094 List<SubunitCluster> c2 = SubunitClusterer.cluster(s2, cParams).getClusters(); 095 096 // STEP 2: match each subunit cluster between groups O(N^2*L^2) - inter 097 Map<Integer, Integer> clusterMap = new HashMap<Integer, Integer>(); 098 for (int i = 0; i < c1.size(); i++) { 099 for (int j = 0; j < c2.size(); j++) { 100 101 if (clusterMap.keySet().contains(i)) 102 break; 103 if (clusterMap.values().contains(j)) 104 continue; 105 106 // Use structural alignment to match the subunit clusters 107 if (c1.get(i).mergeStructure(c2.get(j),cParams)) { 108 clusterMap.put(i, j); 109 } 110 } 111 } 112 113 logger.info("Cluster Map: " + clusterMap.toString()); 114 result.setClusters(c1); 115 116 // STEP 3: Align the assemblies for each cluster match O(N^2*L) 117 for (int globalKey : clusterMap.keySet()) { 118 119 // Obtain the clusters 120 SubunitCluster clust1 = c1.get(globalKey); 121 SubunitCluster clust2 = c2.get(clusterMap.get(globalKey)); 122 123 // Take a cluster match as reference 124 int index1 = 0; 125 int index2 = clust1.size() - clust2.size(); 126 Map<Integer, Integer> subunitMap = new HashMap<Integer, Integer>(); 127 subunitMap.put(index1, index2); 128 129 // Map cluster id to their subunit matching 130 Map<Integer, Map<Integer, Integer>> clustSubunitMap = new HashMap<Integer, Map<Integer, Integer>>(); 131 clustSubunitMap.put(globalKey, subunitMap); 132 133 // Change order of key set so that globalKey is first 134 List<Integer> keySet = new ArrayList<Integer>(clusterMap.keySet()); 135 keySet.remove((Integer) globalKey); 136 keySet.add(0, globalKey); 137 138 for (int key : clusterMap.keySet()) { 139 140 // Recover subunitMap if it is the reference, new one otherwise 141 if (key == globalKey) 142 subunitMap = clustSubunitMap.get(key); 143 else 144 subunitMap = new HashMap<Integer, Integer>(); 145 146 // Obtain the clusters of each subunit group 147 clust1 = c1.get(key); 148 clust2 = c2.get(clusterMap.get(key)); 149 150 // Get the initial subunit indices of each group 151 index1 = 0; 152 index2 = clust1.size() - clust2.size(); 153 154 for (int i = 0; i < index2; i++) { 155 for (int j = index2; j < clust1.size(); j++) { 156 157 if (subunitMap.keySet().contains(i)) 158 break; 159 if (subunitMap.values().contains(j)) 160 continue; 161 162 // Obtain cumulative transformation matrix 163 Matrix4d transform = getTransformForClusterSubunitMap( 164 c1, clustSubunitMap); 165 166 // Obtain Atom arrays of the subunit pair to match 167 Atom[] atoms1 = clust1.getAlignedAtomsSubunit(i); 168 Atom[] atoms2 = clust1.getAlignedAtomsSubunit(j); 169 170 // Obtain centroids and transform the second 171 Atom centr1 = Calc.getCentroid(atoms1); 172 Atom centr2 = Calc.getCentroid(atoms2); 173 Calc.transform(centr2, transform); 174 175 // 1- Check that the distance fulfills maximum 176 double dCentroid = Calc.getDistance(centr1, centr2); 177 if (dCentroid > aParams.getdCutoff()) { 178 logger.debug(String.format("Subunit matching %d " 179 + "vs %d of cluster %d could not be " 180 + "matched, because centroid distance is " 181 + "%.2f", index1, index2, key, dCentroid)); 182 continue; 183 } 184 185 // Transform coordinates of second 186 Atom[] atoms2c = StructureTools.cloneAtomArray(atoms2); 187 Calc.transform(atoms2c, transform); 188 189 // 2- Check the orientation metric condition 190 double qOrient = UnitQuaternions.orientationAngle( 191 Calc.atomsToPoints(atoms1), 192 Calc.atomsToPoints(atoms2c), false); 193 qOrient = Math.min(Math.abs(2*Math.PI - qOrient), qOrient); 194 if (qOrient > aParams.getMaxOrientationAngle()) { 195 logger.debug(String.format("Subunit matching %d " 196 + "vs %d of cluster %d could not be " 197 + "matched, because orientation metric is " 198 + "%.2f", i, j, key, qOrient)); 199 continue; 200 } 201 202 // 3- Check the RMSD condition 203 double rmsd = Calc.rmsd(atoms1, atoms2c); 204 if (rmsd > aParams.getMaxRmsd()) { 205 logger.debug(String.format("Subunit matching %d " 206 + "vs %d of cluster %d could not be " 207 + "matched, because RMSD is %.2f", i, 208 j, key, rmsd)); 209 continue; 210 } 211 212 logger.info(String.format("Subunit matching %d vs %d" 213 + " of cluster %d with centroid distance %.2f" 214 + ", orientation metric %.2f and RMSD %.2f", 215 i, j, key, dCentroid, qOrient, rmsd)); 216 217 subunitMap.put(i, j); 218 } 219 } 220 221 clustSubunitMap.put(key, subunitMap); 222 } 223 224 logger.info("Cluster Subunit Map: " + clustSubunitMap.toString()); 225 226 // Unfold the nested map into subunit map and alignment 227 subunitMap = new HashMap<Integer, Integer>(); 228 List<Integer> alignRes1 = new ArrayList<Integer>(); 229 List<Integer> alignRes2 = new ArrayList<Integer>(); 230 List<Atom> atomArray1 = new ArrayList<Atom>(); 231 List<Atom> atomArray2 = new ArrayList<Atom>(); 232 233 for (int key : clustSubunitMap.keySet()) { 234 235 // Obtain the cluster and the alignment in it 236 SubunitCluster cluster = c1.get(key); 237 List<List<Integer>> clusterEqrs = cluster 238 .getMultipleAlignment().getBlock(0).getAlignRes(); 239 240 for (Entry<Integer, Integer> pair : clustSubunitMap.get(key) 241 .entrySet()) { 242 243 int i = pair.getKey(); 244 int j = pair.getValue(); 245 246 // Obtain the indices of the original Subunit Lists 247 int orig1 = s1.indexOf(cluster.getSubunits().get(i)); 248 int orig2 = s2.indexOf(cluster.getSubunits().get(j)); 249 250 // Append rescaled aligned residue indices 251 for (Integer eqr : clusterEqrs.get(i)) 252 alignRes1.add(eqr + atomArray1.size()); 253 for (Integer eqr : clusterEqrs.get(j)) 254 alignRes2.add(eqr + atomArray2.size()); 255 256 // Apend atoms to the arrays 257 atomArray1.addAll(Arrays.asList(s1.get(orig1) 258 .getRepresentativeAtoms())); 259 atomArray2.addAll(Arrays.asList(s2.get(orig2) 260 .getRepresentativeAtoms())); 261 262 subunitMap.put(orig1, orig2); 263 } 264 } 265 266 // Evaluate the goodness of the match with an alignment object 267 MultipleAlignment msa = new MultipleAlignmentImpl(); 268 msa.setEnsemble(new MultipleAlignmentEnsembleImpl()); 269 msa.getEnsemble().setAtomArrays( 270 Arrays.asList(new Atom[][] { 271 atomArray1.toArray(new Atom[atomArray1.size()]), 272 atomArray2.toArray(new Atom[atomArray2.size()]) })); 273 274 // Fill in the alignment information 275 BlockSet bs = new BlockSetImpl(msa); 276 Block b = new BlockImpl(bs); 277 List<List<Integer>> alignRes = new ArrayList<List<Integer>>(2); 278 alignRes.add(alignRes1); 279 alignRes.add(alignRes2); 280 b.setAlignRes(alignRes); 281 282 // Fill in the transformation matrices 283 new ReferenceSuperimposer().superimpose(msa); 284 285 // Calculate some scores 286 MultipleAlignmentScorer.calculateScores(msa); 287 288 // If it is the best match found so far store it 289 if (subunitMap.size() > result.getSubunitMap().size()) { 290 result.setSubunitMap(subunitMap); 291 result.setAlignment(msa); 292 logger.info("Better result found: " + result.toString()); 293 } else if (subunitMap.size() == result.getSubunitMap().size()) { 294 if (result.getAlignment() == null) { 295 result.setSubunitMap(subunitMap); 296 result.setAlignment(msa); 297 } else if (msa.getScore(MultipleAlignmentScorer.RMSD) < result 298 .getRmsd()) { 299 result.setSubunitMap(subunitMap); 300 result.setAlignment(msa); 301 logger.info("Better result found: " + result.toString()); 302 } 303 } 304 305 } 306 307 return result; 308 } 309 310 /** 311 * Returns a pair of Atom arrays corresponding to the alignment of subunit 312 * matchings, in order of appearance. Superposition of the two Atom sets 313 * gives the transformation of the complex. 314 * <p> 315 * Utility method to cumulative calculate the alignment Atoms. 316 * 317 * @param clusters 318 * List of SubunitClusters 319 * @param clusterSubunitMap 320 * map from cluster id to subunit matching 321 * @return pair of atom arrays to be superposed 322 */ 323 private static Pair<Atom[]> getAlignedAtomsForClusterSubunitMap( 324 List<SubunitCluster> clusters, 325 Map<Integer, Map<Integer, Integer>> clusterSubunitMap) { 326 327 List<Atom> atomArray1 = new ArrayList<Atom>(); 328 List<Atom> atomArray2 = new ArrayList<Atom>(); 329 330 // For each cluster of subunits 331 for (int key : clusterSubunitMap.keySet()) { 332 333 // Obtain the cluster and the alignment in it 334 SubunitCluster cluster = clusters.get(key); 335 336 // For each subunit matching in the cluster 337 for (Entry<Integer, Integer> pair : clusterSubunitMap.get(key) 338 .entrySet()) { 339 340 int i = pair.getKey(); 341 int j = pair.getValue(); 342 343 // Apend atoms to the arrays 344 atomArray1.addAll(Arrays.asList(cluster 345 .getAlignedAtomsSubunit(i))); 346 atomArray2.addAll(Arrays.asList(cluster 347 .getAlignedAtomsSubunit(j))); 348 } 349 350 } 351 return new Pair<Atom[]>( 352 atomArray1.toArray(new Atom[atomArray1.size()]), 353 atomArray2.toArray(new Atom[atomArray2.size()])); 354 } 355 356 /** 357 * Returns the transformation matrix corresponding to the alignment of 358 * subunit matchings. 359 * <p> 360 * Utility method to cumulative calculate the alignment transformation. 361 * 362 * @param clusters 363 * List of SubunitClusters 364 * @param clusterSubunitMap 365 * map from cluster id to subunit matching 366 * @return transformation matrix 367 * @throws StructureException 368 */ 369 private static Matrix4d getTransformForClusterSubunitMap( 370 List<SubunitCluster> clusters, 371 Map<Integer, Map<Integer, Integer>> clusterSubunitMap) 372 throws StructureException { 373 374 Pair<Atom[]> pair = getAlignedAtomsForClusterSubunitMap(clusters, 375 clusterSubunitMap); 376 377 return SuperPositions.superpose(Calc.atomsToPoints(pair.getFirst()), 378 Calc.atomsToPoints(pair.getSecond())); 379 380 } 381}