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 - def __init__(self):
98 self._clear()
99 100 # Private methods 101
102 - def _clear(self):
103 self.reference_coords = None 104 self.coords = None 105 self.transformed_coords = None 106 self.rot = None 107 self.tran = None 108 self.rms = None 109 self.init_rms = None
110
111 - def _rms(self, coords1, coords2):
112 """Return rms deviations between coords1 and coords2.""" 113 diff = coords1 - coords2 114 l = coords1.shape[0] 115 return sqrt(sum(sum(diff * diff)) / l)
116 117 # Public methods 118
119 - def set(self, reference_coords, coords):
120 """Set the coordinates to be superimposed. 121 122 coords will be put on top of reference_coords. 123 124 - reference_coords: an NxDIM array 125 - coords: an NxDIM array 126 127 DIM is the dimension of the points, N is the number 128 of points to be superimposed. 129 """ 130 # clear everything from previous runs 131 self._clear() 132 # store cordinates 133 self.reference_coords = reference_coords 134 self.coords = coords 135 n = reference_coords.shape 136 m = coords.shape 137 if n != m or not(n[1] == m[1] == 3): 138 raise Exception("Coordinate number/dimension mismatch.") 139 self.n = n[0]
140
141 - def run(self):
142 """Superimpose the coordinate sets.""" 143 if self.coords is None or self.reference_coords is None: 144 raise Exception("No coordinates set.") 145 coords = self.coords 146 reference_coords = self.reference_coords 147 # center on centroid 148 av1 = sum(coords) / self.n 149 av2 = sum(reference_coords) / self.n 150 coords = coords - av1 151 reference_coords = reference_coords - av2 152 # correlation matrix 153 a = dot(transpose(coords), reference_coords) 154 u, d, vt = svd(a) 155 self.rot = transpose(dot(transpose(vt), transpose(u))) 156 # check if we have found a reflection 157 if det(self.rot) < 0: 158 vt[2] = -vt[2] 159 self.rot = transpose(dot(transpose(vt), transpose(u))) 160 self.tran = av2 - dot(av1, self.rot)
161
162 - def get_transformed(self):
163 """Get the transformed coordinate set.""" 164 if self.coords is None or self.reference_coords is None: 165 raise Exception("No coordinates set.") 166 if self.rot is None: 167 raise Exception("Nothing superimposed yet.") 168 if self.transformed_coords is None: 169 self.transformed_coords = dot(self.coords, self.rot) + self.tran 170 return self.transformed_coords
171
172 - def get_rotran(self):
173 """Right multiplying rotation matrix and translation.""" 174 if self.rot is None: 175 raise Exception("Nothing superimposed yet.") 176 return self.rot, self.tran
177
178 - def get_init_rms(self):
179 """Root mean square deviation of untransformed coordinates.""" 180 if self.coords is None: 181 raise Exception("No coordinates set yet.") 182 if self.init_rms is None: 183 self.init_rms = self._rms(self.coords, self.reference_coords) 184 return self.init_rms
185
186 - def get_rms(self):
187 """Root mean square deviation of superimposed coordinates.""" 188 if self.rms is None: 189 transformed_coords = self.get_transformed() 190 self.rms = self._rms(transformed_coords, self.reference_coords) 191 return self.rms
192 193 194 if __name__ == "__main__": 195 from Bio._utils import run_doctest 196 run_doctest(verbose=0) 197