Package Pmv :: Module Ribbon
[hide private]
[frames] | no frames]

Source Code for Module Pmv.Ribbon

  1  ############################################################################# 
  2  # 
  3  # Author: Sophie Coon, Michel F. SANNER 
  4  # 
  5  # Copyright: M. Sanner TSRI 2000 
  6  # 
  7  ############################################################################# 
  8   
  9  # 
 10  # $Header: /opt/cvs/python/packages/share1.5/Pmv/Ribbon.py,v 1.12 2007/05/14 23:30:34 vareille Exp $ 
 11  # 
 12  # $Id: Ribbon.py,v 1.12 2007/05/14 23:30:34 vareille Exp $ 
 13  # 
 14   
 15  import Numeric 
 16   
17 -def cross( b, c):
18 a = [0.,0.,0.] 19 a[0]=b[1]*c[2]-c[1]*b[2] 20 a[1]=b[2]*c[0]-c[2]*b[0] 21 a[2]=b[0]*c[1]-c[0]*b[1] 22 return a
23
24 -def normalize(v):
25 import math 26 norm = 1.0 / math.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]) 27 v[0] = v[0]*norm 28 v[1] = v[1]*norm 29 v[2] = v[2]*norm 30 return v
31
32 -def bspline( v1, v2, v3, t ):
33 v4 = [0.,0.,0.,0.] 34 frac3 = 0.5 * t*t 35 frac1 = 0.5 * (1.-t) * (1.-t) 36 frac2 = 1. - (frac1 + frac3) 37 for i in (0,1,2,3): 38 v4[i] = frac1 * v1[i] + frac2 * v2[i] + frac3 * v3[i] 39 return v4
40
41 -def interpPoints( ctrl, nchord):
42 43 assert ctrl.shape[1]==4 44 tinc = 1./float(nchord-1) 45 46 smooth = Numeric.zeros( (((ctrl.shape[0]-2)*nchord+2),4), 'f') 47 if smooth is None: return None 48 49 # calculate spline segments 50 # first point in rib 51 smooth[0] = ctrl[0] 52 53 # splines go midpoint-to-midpoint 54 nb = 1 55 for pt in range(1,ctrl.shape[0]-1): 56 t=0.0 57 for i in range(nchord): 58 smooth[nb] = bspline( ctrl[pt-1], ctrl[pt], ctrl[pt+1], t) 59 nb = nb+1 60 t = t + tinc 61 # last point in rib 62 smooth[nb] = ctrl[-1] 63 64 return smooth
65 66 67 """ 68 ctrl: ctrl points 4D 69 nchord: how many interpolations per ctrl pt 70 71 splining from Larry Andrews 7-Nov-1988 72 """
73 -def ribdrw( ctrl, nrib, nchord ):
74 75 tinc = 1./nchord 76 smooth = Numeric.zeros( (nrib, ((ctrl.shape[0]-2)*nchord+2), 4 ), 'd') 77 78 # calculate spline segments 79 for irib in range(nrib): 80 81 # first point in rib 82 smooth[irib][0] = ctrl[0][irib] 83 84 # splines go midpoint-to-midpoint 85 nb = 1 86 for pt in range(1,ctrl.shape[0]-1): 87 t = 0.0 88 for i in range(nchord): 89 smooth[irib][nb] = bspline( ctrl[pt-1][irib], ctrl[pt][irib], 90 ctrl[pt+1][irib], t ) 91 t = t + tinc 92 nb = nb+1 93 94 # last point in rib 95 smooth[irib][nb] = ctrl[-1][irib]; 96 97 return smooth
98 99
100 -def ribbon2D( nrib, ribwid, nchords, offset, natoms, coords, isHelix, 101 off_c = 0.5):
102 """ 103 Generate ctrl points for protein ribbon, based on ideas on 104 Carson & Bugg, J.Molec.Graphics 4,121-122 (1986) 105 106 Ctrl points for Bspline are generated along a line passing 107 through each CA and along the average of the two peptide planes 108 109 nrib number of strands in ribbon 110 ribwid total ribbon width 111 nchord number of chords/residue 112 offset amount to offset ctrl points away from CA positions 113 natom number of atoms stored in arrays 114 coords coordinates of CA and O atoms of all residues in strand 115 isHelix set to 1 for residues in alpha helices 116 off_c this parameter has been added to account for DNA/RNA in the coil 117 """ 118 119 coords = Numeric.array(coords) 120 121 # C Strand separation and half 122 nres = natoms/2; 123 124 ## ctrl = Numeric.zeros( (nres, nrib, 4), 'd') 125 # MS changes 126 ctrl = Numeric.zeros( (nres+1, nrib, 4), 'd') 127 128 E = Numeric.zeros( (3,), 'd') 129 F = Numeric.zeros( (3,), 'd') 130 G = Numeric.zeros( (3,), 'd') 131 H = Numeric.zeros( (3,), 'd') 132 133 if natoms<=0: return None 134 if nrib > 1: drib = ribwid / (nrib-1) 135 rib2 = (nrib+1)*0.5 136 137 for nr in range(nres): 138 i=2*nr 139 if nr<nres-1: 140 # A is vector CAi to CAi+1 141 A = coords[i+2] - coords[i] 142 143 # B is vector CAi to Oi 144 B = coords[i+1] - coords[i] 145 146 # C = A x B; D = C x A 147 C = Numeric.array(cross(A,B)) 148 D = Numeric.array(cross(C,A)) 149 D = normalize(D) 150 151 if i==0: 152 E[:] = D[:] # First peptide, no previous one to average with 153 P = Numeric.zeros( (3,), 'd') # No offset for first CA 154 else: 155 # Not first, ribbon cross vector is average of peptide plane 156 # with previous one 157 if Numeric.dot(D,G)<0.0: B = -D 158 else: B = D 159 E = G+B 160 # Offset is along bisector of CA-CA-CA vectors A (H is Ai-1) 161 P = H-A 162 P = normalize(P) 163 else: 164 E[:] = G[:] # Last one, just use last plane 165 P = Numeric.zeros( (3,), 'd') # No offset for last CA 166 167 E = normalize(E) # Normalise vector E 168 169 # MS oct 11 2001. add a first ctrl point per rib the is before the 170 # first CA atom. Use -0.5*A as the offset along the chain 171 if i==0: 172 # Generate ctrl points 173 if not isHelix[nr]: 174 P = coords[i] - 0.5*A 175 for j in range(nrib): 176 fr=(float(j+1)-rib2)*drib 177 F = fr*E 178 ctrl[0][j][:3] = P+F 179 ctrl[0][j][3] = i+2; 180 181 else: 182 # is secondary structure starts with an Helix 183 # The first point in an helix is particular ... 184 for j in range(nrib): 185 fr=(float(j+1)-rib2)*drib 186 F = fr*E 187 ctrl[0][j][:3] = coords[i]-0.1 *A +F 188 ctrl[0][j][3] = i+2; 189 190 # Generate ctrl points 191 if not isHelix[nr]: # control points for none helical parts 192 P = coords[i] + off_c*A 193 else: # control point for helices 194 # increasing the offset increases the radius of the helix 195 # increasing the factor multiplying A shifts the per/residue 196 # point down the helix 197 P = (2*offset) * P + 0.75*A 198 P = coords[i]+P 199 200 for j in range(nrib): 201 fr=(float(j+1)-rib2)*drib 202 F = fr*E 203 # MS changed. 204 ctrl[nr+1][j][:3] = P+F 205 ctrl[nr+1][j][3] = i+2; 206 207 # Store things for next residue 208 G[:] = E[:] 209 H[:] = A[:] 210 # end loop over residues 211 212 smooth = ribdrw(ctrl, nrib, nchords) 213 return smooth#, ctrl
214 215 if __name__=='__main__': 216 217 cao = Numeric.array( 218 (( 16.966999, 12.784000, 4.338000 ), (15.268000, 13.825000, 5.594000), 219 ( 13.856000, 11.469000, 6.066000 ), (14.993000, 9.862000, 7.443000), 220 ( 13.660000, 10.707000, 9.787000 ), (11.393000, 11.308000, 10.185000), 221 ( 10.646000, 8.991000, 11.408000 ), (11.659000, 8.296000, 13.491000), 222 ( 9.448000, 9.034000, 15.012000 ), (9.490000, 7.519000, 16.819000), 223 ( 8.673000, 5.314000, 15.279000 ), (8.726000, 4.858000, 12.923000), 224 ( 8.912000, 2.083000, 13.258000 ), (7.670000, 2.031000, 11.245000), 225 ( 5.145000, 2.209000, 12.453000 ), (4.664000, 3.268000, 10.343000), 226 ( 5.598000, 5.767000, 11.082000 ), (6.052000, 5.933000, 8.744000), 227 ( 8.496000, 4.609000, 8.837000 ), (7.878000, 3.778000, 6.651000), 228 ( 6.500000, 1.584000, 7.565000 ), (5.213000, 2.016000, 5.557000), 229 ( 3.545000, 3.935000, 6.751000 ), (3.536000, 5.001000, 4.617000), 230 ( 5.929000, 6.358000, 5.055000 ), (6.136000, 6.072000, 2.653000), 231 ( 7.331000, 3.607000, 2.791000 ), (6.240000, 3.144000, 0.684000), 232 ( 3.782000, 2.599000, 1.742000 ), (2.947000, 3.817000, -0.189000), 233 ( 2.890000, 6.285000, 1.126000 ), (3.200000, 7.147000, -1.103000), 234 ( 5.895000, 6.489000, -1.213000 ), (6.228000, 5.901000, -3.507000), 235 ( 4.933000, 3.431000, -3.326000 ), (4.988000, 3.755000, -5.687000), 236 ( 2.792000, 5.376000, -5.797000 ), (3.260000, 7.045000, -7.422000), 237 ( 5.366000, 8.191000, -6.018000 ), (5.535000, 10.510000, -5.730000), 238 ( 3.767000, 10.609000, -3.513000 ), (5.947000, 10.757000, -2.523000), 239 ( 6.143000, 13.513000, -2.696000 ), (5.485000, 13.061000, -0.382000), 240 ( 8.114000, 13.103000, 0.500000 ), (7.036000, 13.682000, 2.540000), 241 ( 6.614000, 16.316999, 1.913000 ), (4.782000, 16.166000, 3.495000), 242 ( 3.074000, 14.894000, 1.756000 ), (2.315000, 13.523000, 3.578000), 243 ( 4.180000, 11.549000, 3.187000 ), (4.227000, 11.252000, 5.547000), 244 ( 5.879000, 13.502000, 6.026000 ), (4.528000, 13.422000, 8.025000), 245 ( 2.691000, 15.221000, 7.194000 ), (0.947000, 14.112000, 8.468000), 246 ( 0.715000, 12.045000, 6.657000 ), (0.286000, 10.632000, 8.545000), 247 ( 2.986000, 9.994000, 8.950000 ), (3.766000, 9.715000, 11.186000), 248 ( 4.769000, 12.336000, 11.360000 ), (7.037000, 12.750000, 11.954000), 249 ( 8.140000, 11.694000, 9.635000 ), (7.581000, 13.949000, 8.944000), 250 ( 10.280000, 14.760000, 8.823000 ), (11.971000, 13.583000, 7.552000), 251 ( 12.552000, 15.877000, 6.036000 ), (13.168000, 18.006001, 6.945000), 252 ( 15.930000, 17.454000, 6.941000 ), (17.097000, 16.660000, 4.970000), 253 ( 18.635000, 18.861000, 4.738000 ), (20.593000, 17.742001, 3.945000), 254 ( 21.452000, 16.969000, 6.513000 ), (20.138000, 15.023000, 5.878000), 255 ( 22.018999, 13.242000, 7.020000 ), (21.868999, 11.387000, 8.435000), 256 ( 21.936001, 12.911000, 10.809000 ), (20.357000, 14.317000, 11.948000), 257 ( 18.504000, 12.312000, 12.298000 ), (19.533001, 11.718000, 14.362000), 258 ( 17.924000, 13.421000, 15.877000 ), (16.652000, 11.368000, 16.033001), 259 ( 17.334000, 10.956000, 18.691000 ), (15.434000, 9.550000, 19.166000), 260 ( 13.564000, 11.573000, 18.836000 ), (11.720000, 11.040000, 17.427999), 261 ( 13.257000, 10.745000, 15.081000 ), (14.930000, 9.862000, 13.568000), 262 ( 15.445000, 7.667000, 15.246000 ), (16.093000, 5.705000, 14.039000), 263 ( 13.512000, 5.395000, 12.878000 ), (13.733000, 6.929000, 11.026000))) 264 265 nchords=10 266 nrib=2 267 nres=46 268 natoms=2*nres 269 offset=1.2 270 width=1.5 271 272 import profile 273 # profile.run("smooth, ctrl = ribbon2D( nrib, width, nchords, offset, natoms, cao, [0]*46 );") 274 profile.run("smooth = ribbon2D( nrib, width, nchords, offset, natoms, cao, [0]*46 );") 275 276 from DejaVu import Viewer 277 vi = Viewer() 278 from DejaVu.IndexedPolygons import IndexedPolygons 279 f = [] 280 n = smooth.shape[1] 281 f = map( lambda x: (x, x+1, x+n+1, x+n), range(smooth.shape[1]-1)) 282 v = Numeric.array(Numeric.reshape(smooth, (-1,4))[:,:3]) 283 p = IndexedPolygons('sheet2D', vertices = v, faces = f, protected=True,) 284 if self.vf.userpref['sharpColorBoundariesForMsms']['value'] == 'blur': 285 p.Set(inheritSharpColorBoundaries=False, sharpColorBoundaries=False,) 286 p.replace = False 287 vi.AddObject(p) 288 289 from DejaVu.Spheres import Spheres 290 ctrl.shape = (-1,4) 291 p = Spheres('ctrl', centers = ctrl[:, :3], radii = 0.6, protected=True) 292 p.replace = False 293 vi.AddObject(p) 294 295 from DejaVu.Spheres import Spheres 296 s = Spheres('sph', vertices = v, quality=5, protected=True) 297 s.radius = 0.2 298 vi.AddObject(s) 299 300 from MolKit import Read 301 mol = Read('../1crn.pdb') 302 allatoms = mol.chains.residues.atoms 303 caat = allatoms.get('CA') 304 coords = caat.coords 305 from DejaVu.Spheres import Spheres 306 p = Spheres('CA', centers = coords, radii = 0.6, protected=True) 307 p.replace = False 308 vi.AddObject(p) 309 310 ## from BioChem.protein import Protein 311 ## from BioChem.pdbParser import PdbParser 312 ## mol = Protein() 313 ## mol.read('/tsri/pdb/struct/1b1g.pdb', PdbParser()) 314 ## l = lambda x:x.name=='CA' or x.name=='O' 315 ## cao1 = mol.chains.residues.atoms.get(l).coords 316 ## natoms = len(cao1)/2 317 ## profile.run("smooth1 = ribbon2D( nrib, width, nchords, offset, natoms, cao1 );") 318 ## f1 = [] 319 ## n1 = smooth1.shape[1] 320 ## f1 = map( lambda x: (x, x+1, x+n+1, x+n), range(smooth1.shape[1]-1)) 321 ## v1 = Numeric.array(Numeric.reshape(smooth1, (-1,4))[:,:3]) 322 ## p1 = IndexedPolygons('sheet2D', vertices = v1, faces = f1) 323 ## vi.AddObject(p1) 324