Ticket #2262: mlp-ghose-united-shifted.py

File mlp-ghose-united-shifted.py, 12.9 KB (added by Tom Goddard, 6 years ago)
Line 
1# vim: set expandtab shiftwidth=4 softtabstop=4:
2
3# === UCSF ChimeraX Copyright ===
4# Copyright 2016 Regents of the University of California.
5# All rights reserved. This software provided pursuant to a
6# license agreement containing restrictions on its disclosure,
7# duplication and use. For details see:
8# http://www.rbvi.ucsf.edu/chimerax/docs/licensing.html
9# This notice must be embedded in or attached to all copies,
10# including partial copies, of the software or any revisions
11# or derivations thereof.
12# === UCSF ChimeraX Copyright ===
13
14def mlp(session, atoms=None, method="fauchere", spacing=1.0, max_distance=5.0, nexp=3.0,
15 color=True, palette=None, range=None, map=False):
16 '''Display Molecular Lipophilic Potential for a single model.
17
18 Parameters
19 ----------
20 atoms : Atoms
21 Color surfaces for these atoms using MLP map. Only amino acid residues are used.
22 method : 'dubost','fauchere','brasseur','buckingham','type5'
23 Distance dependent function to use for calculation
24 spacing : float
25 Grid spacing, default 1 Angstrom.
26 max_distance : float
27 Maximum distance from atom to sum lipophilicity. Default 5 Angstroms.
28 nexp : float
29 The buckingham method uses this numerical exponent.
30 color : bool
31 Whether to color molecular surfaces. They are created if they don't yet exist.
32 palette : Colormap
33 Color palette for coloring surfaces.
34 Default is lipophilicity colormap (orange lipophilic, blue lipophobic).
35 range : 2-tuple of float
36 Range of lipophilicity values defining ends of color map. Default is -20,20
37 map : bool
38 Whether to open a volume model of lipophilicity values
39 '''
40 if atoms is None:
41 from chimerax.atomic import all_atoms
42 atoms = all_atoms(session)
43
44 from chimerax.atomic import Residue
45 patoms = atoms[atoms.residues.polymer_types == Residue.PT_AMINO]
46 if len(patoms) == 0:
47 from chimerax.core.errors import UserError
48 raise UserError('mlp: no amino acids specified')
49
50 if palette is None:
51 from chimerax.core.colors import BuiltinColormaps
52 cmap = BuiltinColormaps['lipophilicity']
53 else:
54 cmap = palette
55 if range is None and not cmap.values_specified:
56 range = (-20,20)
57
58 # Color surfaces by lipophilicity
59 if color:
60 # Compute surfaces if not already created
61 from chimerax.surface import surface
62 surfs = surface(session, patoms)
63 for s in surfs:
64 satoms = s.atoms
65 name = 'mlp ' + s.name.split(maxsplit=1)[0]
66 v = mlp_map(session, satoms, method, spacing, max_distance, nexp, name, open_map = map)
67 from chimerax.surface import color_surfaces_by_map_value
68 color_surfaces_by_map_value(satoms, map = v, palette = cmap, range = range)
69 else:
70 name = 'mlp map'
71 v = mlp_map(session, patoms, method, spacing, max_distance, nexp, name, open_map = map)
72
73
74def register_mlp_command(logger):
75 from chimerax.core.commands import register, CmdDesc, SaveFileNameArg, FloatArg, EnumOf, BoolArg, ColormapArg, ColormapRangeArg
76 from chimerax.atomic import AtomsArg
77 desc = CmdDesc(optional=[('atoms', AtomsArg)],
78 keyword=[('spacing', FloatArg),
79 ('max_distance', FloatArg),
80 ('method', EnumOf(['dubost','fauchere','brasseur','buckingham','type5'])),
81 ('nexp', FloatArg),
82 ('color', BoolArg),
83 ('palette', ColormapArg),
84 ('range', ColormapRangeArg),
85 ('map', BoolArg),
86 ],
87 synopsis='display molecular lipophilic potential for selected models')
88 register('mlp', desc, mlp, logger=logger)
89
90def mlp_map(session, atoms, method, spacing, max_dist, nexp, name, open_map):
91 data, bounds = calculatefimap(atoms, method, spacing, max_dist, nexp)
92
93 # m.pot is 1-dimensional if m.writedxfile() was called. Has indices in x,y,z order.
94 origin = tuple(xmin for xmin,xmax in bounds)
95 s = spacing
96 step = (s,s,s)
97 from chimerax.map.data import ArrayGridData
98 g = ArrayGridData(data, origin, step, name = name)
99 g.polar_values = True
100 from chimerax.map import volume_from_grid_data
101 v = volume_from_grid_data(g, session, open_model = open_map, show_dialog = open_map)
102 v.update_drawings() # Compute surface levels
103 v.set_parameters(surface_colors = [(0, 139/255, 139/255, 1), (184/255, 134/255, 11/255, 1)])
104 return v
105
106#
107# Code below is modified version of pyMLP, eliminating most the of code
108# (unneeded parsing PDB files, writing dx files, ...) and optimizing the calculation speed.
109#
110
111class Defaults(object):
112 """Constants"""
113
114 def __init__(self):
115 self.gridmargin = 10.0
116 self.fidatadefault = { #Default fi table
117 'ALA': {'CB': 0.5595, #fi : lipophilic atomic potential
118 'C': -0.0702,
119 'CA': -0.0971,
120 'O': 0.0067,
121 'N': -0.5549},
122 'ARG': {'C': -0.0702,
123 'CA': -0.0971,
124 'CB': 0.4112,
125 'CG': 0.4112,
126 'CD': 0.1016,
127 'CZ': 0.5442,
128 'N': -0.5549,
129 'NE': -0.0825,
130 'NH1': -0.5055,
131 'NH2': -0.5055,
132 'O': 0.0067},
133 'ASN': {'C': -0.0702,
134 'CA': -0.0971,
135 'CB': 0.1248,
136 'CG': -0.0702,
137 'N': -0.5549,
138 'ND2': -0.6285,
139 'O': 0.0067,
140 'OD1': 0.0067},
141 'ASP': {'C': -0.0702,
142 'CA': -0.0971,
143 'CB': 0.1248,
144 'CG': -0.0702,
145 'N': -0.5549,
146 'O': 0.0067,
147 'OD1': -0.3787,
148 'OD2': -0.3787},
149 'CYS': {'C': -0.0702,
150 'CA': -0.0971,
151 'CB': 0.1016,
152 'N': -0.5549,
153 'O': 0.0067,
154 'SG': 0.5710},
155 'GLN': {'C': -0.0702,
156 'CA': -0.0971,
157 'CB': 0.4112,
158 'CG': 0.1248,
159 'CD': -0.0702,
160 'N': -0.5549,
161 'NE2': -0.6285,
162 'O': 0.0067,
163 'OE1': 0.0067},
164 'GLU': {'C': -0.0702,
165 'CA': -0.0971,
166 'CB': 0.4112,
167 'CG': 0.1248,
168 'CD': -0.0702,
169 'N': -0.5549,
170 'O': 0.0067,
171 'OE1': -0.3787,
172 'OE2': -0.3787},
173 'GLY': {'C': -0.0702,
174 'CA': -0.1118,
175 'O': 0.0067,
176 'N': -0.5549},
177 'HIS': {'C': -0.0702,
178 'CA': -0.0971,
179 'CB': 0.1248,
180 'CG': 0.2661,
181 'CD2': 0.5785,
182 'CE1': 0.2043,
183 'N': -0.5549,
184 'ND1': -0.2060,
185 'NE2': -0.2060,
186 'O': 0.0067},
187 'HYP': {'C': -0.0702,
188 'CA': -0.0971,
189 'CB': 0.4112,
190 'CG': -0.0096,
191 'CD': 0.1016,
192 'N': -0.4813,
193 'O': 0.0067,
194 'OD1': -0.4003},
195 'ILE': {'C': -0.0702,
196 'CA': -0.0971,
197 'CB': 0.0585,
198 'CG1': 0.5462,
199 'CG2': 0.7620,
200 'CD1': 0.7620,
201 'N': -0.5549,
202 'O': 0.0067},
203 'LEU': {'C': -0.0702,
204 'CA': -0.0971,
205 'CB': 0.4112,
206 'CG': 0.0660,
207 'CD1': 0.7620,
208 'CD2': 0.7620,
209 'N': -0.5549,
210 'O': 0.0067},
211 'LYS': {'C': -0.0702,
212 'CA': -0.0971,
213 'CB': 0.4112,
214 'CG': 0.5462,
215 'CD': 0.5462,
216 'CE': 0.1016,
217 'NZ': -0.7335,
218 'N': -0.5549,
219 'O': 0.0067},
220 'MET': {'C': -0.0702,
221 'CA': -0.0971,
222 'CB': 0.4112,
223 'CG': 0.1016,
224 'CE': 0.2223,
225 'N': -0.5549,
226 'O': 0.0067,
227 'SD': 0.6206},
228 'MSE': {'C': -0.0702,
229 'CA': -0.0971,
230 'CB': 0.4112,
231 'CG': 0.1016,
232 'CE': 0.2223,
233 'N': -0.5549,
234 'O': 0.0067,
235 'SE': 0.6901},
236 'UNK': {'C': -0.0702,
237 'CA': -0.0971,
238 'N': -0.5549,
239 'O': 0.0067},
240 'ACE': {'C': -0.0702,
241 'CH3': 0.1299,
242 'O': 0.0067},
243 'NME': {'N': -0.5549,
244 'C': 0.2223},
245 'NH2': {'N': -0.6285},
246 'PCA': {'C': -0.0702,
247 'CA': -0.0971,
248 'CB': 0.4112,
249 'CG': 0.1248,
250 'CD': -0.0702,
251 'N': -0.5549,
252 'O': 0.0067,
253 'OE': 0.0067},
254 'PHE': {'C': -0.0702,
255 'CA': -0.0971,
256 'CB': 0.4112,
257 'CG': 0.1792,
258 'CD1': 0.3650,
259 'CD2': 0.3650,
260 'CE1': 0.3650,
261 'CE2': 0.3650,
262 'CZ': 0.3650,
263 'N': -0.5549,
264 'O': 0.0067},
265 'PRO': {'C': -0.0702,
266 'CA': -0.0971,
267 'CB': 0.4112,
268 'CG': 0.4112,
269 'CD': 0.1016,
270 'N': -0.4813,
271 'O': 0.0067},
272 'SER': {'C': -0.0702,
273 'CA': -0.0971,
274 'CB': 0.1016,
275 'N': -0.5549,
276 'O': 0.0067,
277 'OG': -0.4003},
278 'THR': {'C': -0.0702,
279 'CA': -0.0971,
280 'CB': 0.0086,
281 'CG2': 0.5595,
282 'N': -0.5549,
283 'O': 0.0067,
284 'OG1': -0.4003},
285 'TRP': {'C': -0.0702,
286 'CA': -0.0971,
287 'CB': 0.1248,
288 'CG': 0.1792,
289 'CD1': 0.5185,
290 'CD2': 0.1792,
291 'CE2': 0.1839,
292 'CE3': 0.3650,
293 'CH2': 0.3650,
294 'CZ2': 0.3650,
295 'CZ3': 0.3650,
296 'N': -0.5549,
297 'NE1': 0.0823,
298 'O': 0.0067},
299 'TYR': {'C': -0.0702,
300 'CA': -0.0971,
301 'CB': 0.4112,
302 'CG': 0.1792,
303 'CD1': 0.3650,
304 'CD2': 0.3650,
305 'CE1': 0.3650,
306 'CE2': 0.3650,
307 'CZ': 0.1839,
308 'N': -0.5549,
309 'O': 0.0067,
310 'OH': -0.0563},
311 'VAL': {'C': -0.0702,
312 'CA': -0.0971,
313 'CB': 0.0585,
314 'CG1': 0.7620,
315 'CG2': 0.7620,
316 'N': -0.5549,
317 'O': 0.0067}}
318
319def assignfi(fidata, atoms):
320 """assign fi parameters to each atom in the pdbfile"""
321 n = len(atoms)
322 from numpy import empty, float32
323 fi = empty((n,), float32)
324 resname = atoms.residues.names
325 aname = atoms.names
326 for i in range(n):
327 rname = resname[i]
328 rfidata = fidata.get(rname)
329 if rfidata:
330 fi[i]=rfidata.get(aname[i], 0)
331 return fi
332
333def _griddimcalc(listcoord, spacing, gridmargin):
334 """Determination of the grid dimension"""
335 coordmin = min(listcoord) - gridmargin
336 coordmax = max(listcoord) + gridmargin
337 adjustment = ((spacing - (coordmax - coordmin)) % spacing) / 2.
338 coordmin = coordmin - adjustment
339 coordmax = coordmax + adjustment
340 ngrid = int(round((coordmax - coordmin) / spacing))
341 return coordmin, coordmax, ngrid
342
343def calculatefimap(atoms, method, spacing, max_dist, nexp):
344 """Calculation loop"""
345
346 #grid settings in angstrom
347 gridmargin = Defaults().gridmargin
348 xyz = atoms.scene_coords
349 xmingrid, xmaxgrid, nxgrid = _griddimcalc(xyz[:,0], spacing, gridmargin)
350 ymingrid, ymaxgrid, nygrid = _griddimcalc(xyz[:,1], spacing, gridmargin)
351 zmingrid, zmaxgrid, nzgrid = _griddimcalc(xyz[:,2], spacing, gridmargin)
352 bounds = [[xmingrid, xmaxgrid],
353 [ymingrid, ymaxgrid],
354 [zmingrid, zmaxgrid]]
355 origin = (xmingrid, ymingrid, zmingrid)
356
357 fi_table = Defaults().fidatadefault
358 fi = assignfi(fi_table, atoms)
359
360 from numpy import zeros, float32
361 pot = zeros((nzgrid+1, nygrid+1, nxgrid+1), float32)
362 from ._mlp import mlp_sum
363 mlp_sum(xyz, fi, origin, spacing, max_dist, method, nexp, pot)
364
365 return pot, bounds
366
367def mlp_sum(xyz, fi, origin, spacing, max_dist, method, nexp, pot):
368 if method == 'dubost':
369 computemethod = _dubost
370 elif method == 'fauchere':
371 computemethod = _fauchere
372 elif method == 'brasseur':
373 computemethod = _brasseur
374 elif method == 'buckingham':
375 computemethod = _buckingham
376 elif method == 'type5':
377 computemethod = _type5
378 else:
379 raise ValueError('Unknown lipophilicity method %s\n' % computemethod)
380
381 from numpy import zeros, float32, empty, subtract, sqrt
382 grid_pt = empty((3,), float32)
383 dxyz = empty((len(xyz),3), float32)
384 dist = empty((len(xyz),), float32)
385 nz,ny,nx = pot.shape
386 x0,y0,z0 = origin
387 for k in range(nz):
388 grid_pt[2] = z0 + k * spacing
389 for j in range(ny):
390 grid_pt[1] = y0 + j * spacing
391 for i in range(nx):
392 #Evaluation of the distance between th grid point and each atoms
393 grid_pt[0] = x0 + i * spacing
394 subtract(xyz, grid_pt, dxyz)
395 dxyz *= dxyz
396 dist = dxyz[:,0]
397 dist += dxyz[:,1]
398 dist += dxyz[:,2]
399 sqrt(dist, dist)
400 pot[k,j,i] = computemethod(fi, dist, nexp)
401
402def _dubost(fi, d, n):
403 return (100 * fi / (1 + d)).sum()
404
405def _fauchere(fi, d, n):
406 from numpy import exp
407 return (100 * fi * exp(-d)).sum()
408
409def _brasseur(fi, d, n):
410 #3.1 division is there to remove any units in the equation
411 #3.1A is the average diameter of a water molecule (2.82 -> 3.2)
412 from numpy import exp
413 return (100 * fi * exp(-d/3.1)).sum()
414
415def _buckingham(fi, d, n):
416 return (100 * fi / (d**n)).sum()
417
418def _type5(fi, d, n):
419 from numpy import exp, sqrt
420 return (100 * fi * exp(-sqrt(d))).sum()