Package Bio :: Package SVDSuperimposer
[hide private]
[frames] | no frames]

Source Code for Package Bio.SVDSuperimposer

  1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@vub.ac.be) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5  """Align on protein structure onto another using SVD alignment. 
  6   
  7  SVDSuperimposer finds the best rotation and translation to put 
  8  two point sets on top of each other (minimizing the RMSD). This is 
  9  eg. useful to superimpose crystal structures. SVD stands for singular 
 10  value decomposition, which is used in the algorithm. 
 11  """ 
 12   
 13  from __future__ import print_function 
 14   
 15  from numpy import dot, transpose, sqrt, array 
 16  from numpy.linalg import svd, det 
 17   
 18   
19 -class SVDSuperimposer(object):
20 """Class to run SVD alignment, 21 22 SVDSuperimposer finds the best rotation and translation to put 23 two point sets on top of each other (minimizing the RMSD). This is 24 eg. useful to superimpose crystal structures. 25 26 SVD stands for Singular Value Decomposition, which is used to calculate 27 the superposition. 28 29 Reference: 30 31 Matrix computations, 2nd ed. Golub, G. & Van Loan, CF., The Johns 32 Hopkins University Press, Baltimore, 1989 33 """
34 - def __init__(self):
35 self._clear()
36 37 # Private methods 38
39 - def _clear(self):
40 self.reference_coords = None 41 self.coords = None 42 self.transformed_coords = None 43 self.rot = None 44 self.tran = None 45 self.rms = None 46 self.init_rms = None
47
48 - def _rms(self, coords1, coords2):
49 """Return rms deviations between coords1 and coords2.""" 50 diff = coords1 - coords2 51 l = coords1.shape[0] 52 return sqrt(sum(sum(diff * diff)) / l)
53 54 # Public methods 55
56 - def set(self, reference_coords, coords):
57 """Set the coordinates to be superimposed. 58 59 coords will be put on top of reference_coords. 60 61 - reference_coords: an NxDIM array 62 - coords: an NxDIM array 63 64 DIM is the dimension of the points, N is the number 65 of points to be superimposed. 66 """ 67 # clear everything from previous runs 68 self._clear() 69 # store cordinates 70 self.reference_coords = reference_coords 71 self.coords = coords 72 n = reference_coords.shape 73 m = coords.shape 74 if n != m or not(n[1] == m[1] == 3): 75 raise Exception("Coordinate number/dimension mismatch.") 76 self.n = n[0]
77
78 - def run(self):
79 """Superimpose the coordinate sets.""" 80 if self.coords is None or self.reference_coords is None: 81 raise Exception("No coordinates set.") 82 coords = self.coords 83 reference_coords = self.reference_coords 84 # center on centroid 85 av1 = sum(coords) / self.n 86 av2 = sum(reference_coords) / self.n 87 coords = coords - av1 88 reference_coords = reference_coords - av2 89 # correlation matrix 90 a = dot(transpose(coords), reference_coords) 91 u, d, vt = svd(a) 92 self.rot = transpose(dot(transpose(vt), transpose(u))) 93 # check if we have found a reflection 94 if det(self.rot) < 0: 95 vt[2] = -vt[2] 96 self.rot = transpose(dot(transpose(vt), transpose(u))) 97 self.tran = av2 - dot(av1, self.rot)
98
99 - def get_transformed(self):
100 """Get the transformed coordinate set.""" 101 if self.coords is None or self.reference_coords is None: 102 raise Exception("No coordinates set.") 103 if self.rot is None: 104 raise Exception("Nothing superimposed yet.") 105 if self.transformed_coords is None: 106 self.transformed_coords = dot(self.coords, self.rot) + self.tran 107 return self.transformed_coords
108
109 - def get_rotran(self):
110 """Right multiplying rotation matrix and translation.""" 111 if self.rot is None: 112 raise Exception("Nothing superimposed yet.") 113 return self.rot, self.tran
114
115 - def get_init_rms(self):
116 """Root mean square deviation of untransformed coordinates.""" 117 if self.coords is None: 118 raise Exception("No coordinates set yet.") 119 if self.init_rms is None: 120 self.init_rms = self._rms(self.coords, self.reference_coords) 121 return self.init_rms
122
123 - def get_rms(self):
124 """Root mean square deviation of superimposed coordinates.""" 125 if self.rms is None: 126 transformed_coords = self.get_transformed() 127 self.rms = self._rms(transformed_coords, self.reference_coords) 128 return self.rms
129 130 131 if __name__ == "__main__": 132 133 # start with two coordinate sets (Nx3 arrays - float) 134 135 x = array([[51.65, -1.90, 50.07], 136 [50.40, -1.23, 50.65], 137 [50.68, -0.04, 51.54], 138 [50.22, -0.02, 52.85]], 'f') 139 140 y = array([[51.30, -2.99, 46.54], 141 [51.09, -1.88, 47.58], 142 [52.36, -1.20, 48.03], 143 [52.71, -1.18, 49.38]], 'f') 144 145 # start! 146 sup = SVDSuperimposer() 147 148 # set the coords 149 # y will be rotated and translated on x 150 sup.set(x, y) 151 152 # do the lsq fit 153 sup.run() 154 155 # get the rmsd 156 rms = sup.get_rms() 157 158 # get rotation (right multiplying!) and the translation 159 rot, tran = sup.get_rotran() 160 161 # rotate y on x 162 y_on_x1 = dot(y, rot) + tran 163 164 # same thing 165 y_on_x2 = sup.get_transformed() 166 167 print(y_on_x1) 168 print("") 169 print(y_on_x2) 170 print("") 171 print("%.2f" % rms) 172