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  try: 
 16      from numpy import dot, transpose, sqrt 
 17      from numpy.linalg import svd, det 
 18  except ImportError: 
 19      from Bio import MissingPythonDependencyError 
 20      raise MissingPythonDependencyError( 
 21          "Install NumPy if you want to use Bio.SVDSuperimposer.") 
 22   
 23   
24 -class SVDSuperimposer(object):
25 """Class to run SVD alignment, 26 27 SVDSuperimposer finds the best rotation and translation to put 28 two point sets on top of each other (minimizing the RMSD). This is 29 eg. useful to superimpose crystal structures. 30 31 SVD stands for Singular Value Decomposition, which is used to calculate 32 the superposition. 33 34 Reference: 35 36 Matrix computations, 2nd ed. Golub, G. & Van Loan, CF., The Johns 37 Hopkins University Press, Baltimore, 1989 38 39 start with two coordinate sets (Nx3 arrays - float) 40 41 >>> from Bio.SVDSuperimposer import SVDSuperimposer 42 >>> from numpy import array, dot, set_printoptions 43 >>> 44 >>> x = array([[51.65, -1.90, 50.07], 45 ... [50.40, -1.23, 50.65], 46 ... [50.68, -0.04, 51.54], 47 ... [50.22, -0.02, 52.85]], 'f') 48 >>> 49 >>> y = array([[51.30, -2.99, 46.54], 50 ... [51.09, -1.88, 47.58], 51 ... [52.36, -1.20, 48.03], 52 ... [52.71, -1.18, 49.38]], 'f') 53 54 start 55 56 >>> sup = SVDSuperimposer() 57 58 set the coords y will be rotated and translated on x 59 60 >>> sup.set(x, y) 61 62 do the lsq fit 63 64 >>> sup.run() 65 66 get the rmsd 67 68 >>> rms = sup.get_rms() 69 70 get rotation (right multiplying!) and the translation 71 72 >>> rot, tran = sup.get_rotran() 73 74 rotate y on x 75 76 >>> y_on_x1 = dot(y, rot) + tran 77 78 same thing 79 80 >>> y_on_x2 = sup.get_transformed() 81 82 >>> set_printoptions(precision=2) 83 >>> print(y_on_x1) 84 [[ 5.17e+01 -1.90e+00 5.01e+01] 85 [ 5.04e+01 -1.23e+00 5.06e+01] 86 [ 5.07e+01 -4.16e-02 5.15e+01] 87 [ 5.02e+01 -1.94e-02 5.29e+01]] 88 >>> print(y_on_x2) 89 [[ 5.17e+01 -1.90e+00 5.01e+01] 90 [ 5.04e+01 -1.23e+00 5.06e+01] 91 [ 5.07e+01 -4.16e-02 5.15e+01] 92 [ 5.02e+01 -1.94e-02 5.29e+01]] 93 >>> print("%.2f" % rms) 94 0.00 95 96 """ 97
98 - def __init__(self):
99 self._clear()
100 101 # Private methods 102
103 - def _clear(self):
104 self.reference_coords = None 105 self.coords = None 106 self.transformed_coords = None 107 self.rot = None 108 self.tran = None 109 self.rms = None 110 self.init_rms = None
111
112 - def _rms(self, coords1, coords2):
113 """Return rms deviations between coords1 and coords2.""" 114 diff = coords1 - coords2 115 l = coords1.shape[0] 116 return sqrt(sum(sum(diff * diff)) / l)
117 118 # Public methods 119
120 - def set(self, reference_coords, coords):
121 """Set the coordinates to be superimposed. 122 123 coords will be put on top of reference_coords. 124 125 - reference_coords: an NxDIM array 126 - coords: an NxDIM array 127 128 DIM is the dimension of the points, N is the number 129 of points to be superimposed. 130 """ 131 # clear everything from previous runs 132 self._clear() 133 # store cordinates 134 self.reference_coords = reference_coords 135 self.coords = coords 136 n = reference_coords.shape 137 m = coords.shape 138 if n != m or not(n[1] == m[1] == 3): 139 raise Exception("Coordinate number/dimension mismatch.") 140 self.n = n[0]
141
142 - def run(self):
143 """Superimpose the coordinate sets.""" 144 if self.coords is None or self.reference_coords is None: 145 raise Exception("No coordinates set.") 146 coords = self.coords 147 reference_coords = self.reference_coords 148 # center on centroid 149 av1 = sum(coords) / self.n 150 av2 = sum(reference_coords) / self.n 151 coords = coords - av1 152 reference_coords = reference_coords - av2 153 # correlation matrix 154 a = dot(transpose(coords), reference_coords) 155 u, d, vt = svd(a) 156 self.rot = transpose(dot(transpose(vt), transpose(u))) 157 # check if we have found a reflection 158 if det(self.rot) < 0: 159 vt[2] = -vt[2] 160 self.rot = transpose(dot(transpose(vt), transpose(u))) 161 self.tran = av2 - dot(av1, self.rot)
162
163 - def get_transformed(self):
164 """Get the transformed coordinate set.""" 165 if self.coords is None or self.reference_coords is None: 166 raise Exception("No coordinates set.") 167 if self.rot is None: 168 raise Exception("Nothing superimposed yet.") 169 if self.transformed_coords is None: 170 self.transformed_coords = dot(self.coords, self.rot) + self.tran 171 return self.transformed_coords
172
173 - def get_rotran(self):
174 """Right multiplying rotation matrix and translation.""" 175 if self.rot is None: 176 raise Exception("Nothing superimposed yet.") 177 return self.rot, self.tran
178
179 - def get_init_rms(self):
180 """Root mean square deviation of untransformed coordinates.""" 181 if self.coords is None: 182 raise Exception("No coordinates set yet.") 183 if self.init_rms is None: 184 self.init_rms = self._rms(self.coords, self.reference_coords) 185 return self.init_rms
186
187 - def get_rms(self):
188 """Root mean square deviation of superimposed coordinates.""" 189 if self.rms is None: 190 transformed_coords = self.get_transformed() 191 self.rms = self._rms(transformed_coords, self.reference_coords) 192 return self.rms
193 194 195 if __name__ == "__main__": 196 from Bio._utils import run_doctest 197 run_doctest(verbose=0) 198