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

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