Ticket #2262: mlp-orig.py

File mlp-orig.py, 12.2 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.63, #fi : lipophilic atomic potential
118 'C': -0.54,
119 'CA': 0.02,
120 'O': -0.68,
121 'N': -0.44},
122 'ARG': {'C': -0.54,
123 'CA': 0.02,
124 'CB': 0.45,
125 'CD': 0.45,
126 'CG': 0.45,
127 'CZ': -0.54,
128 'N': -0.44,
129 'NE': -0.55,
130 'NH1': -0.11,
131 'NH2': -0.83,
132 'O': -0.68},
133 'ASN': {'C': -0.54,
134 'CA': 0.02,
135 'CB': 0.02,
136 'CG': 0.45,
137 'N': -0.44,
138 'ND2': -0.11,
139 'O': -0.68,
140 'OD1': -0.68},
141 'ASP': {'C': -0.54,
142 'CA': 0.02,
143 'CB': 0.45,
144 'CG': 0.54,
145 'N': -0.44,
146 'O': -0.68,
147 'OD1': -0.68,
148 'OD2': 0.53},
149 'CYS': {'C': -0.54,
150 'CA': 0.02,
151 'CB': 0.45,
152 'N': -0.44,
153 'O': -0.68,
154 'SG': 0.27},
155 'GLN': {'C': -0.54,
156 'CA': 0.02,
157 'CB': 0.45,
158 'CD': -0.54,
159 'CG': 0.45,
160 'N': -0.44,
161 'NE2': -0.11,
162 'O': -0.68,
163 'OE1': -0.68},
164 'GLU': {'C': -0.54,
165 'CA': 0.02,
166 'CB': 0.45,
167 'CD': -0.54,
168 'CG': 0.45,
169 'N': -0.44,
170 'O': -0.68,
171 'OE1': -0.68,
172 'OE2': 0.53},
173 'GLY': {'C': -0.54,
174 'CA': 0.45,
175 'O': -0.68,
176 'N': -0.55},
177 'HIS': {'C': -0.54,
178 'CA': 0.02,
179 'CB': 0.45,
180 'CD2': 0.31,
181 'CE1': 0.31,
182 'CG': 0.09,
183 'N': -0.44,
184 'ND1': -0.56,
185 'NE2': -0.80,
186 'O': -0.68},
187 'HYP': {'C': -0.54,
188 'CA': 0.02,
189 'CB': 0.45,
190 'CD': 0.45,
191 'CG': 0.02,
192 'N': -0.92,
193 'O': -0.68,
194 'OD1': -0.93},
195 'ILE': {'C': -0.54,
196 'CA': 0.02,
197 'CB': 0.02,
198 'CD': 0.63,
199 'CD1': 0.63,
200 'CG1': 0.45,
201 'CG2': 0.63,
202 'N': -0.44,
203 'O': -0.68},
204 'LEU': {'C': -0.54,
205 'CA': 0.02,
206 'CB': 0.45,
207 'CD1': 0.63,
208 'CD2': 0.63,
209 'CG': 0.02,
210 'N': -0.44,
211 'O': -0.68},
212 'LYS': {'C': -0.54,
213 'CA': 0.02,
214 'CB': 0.45,
215 'CD': 0.45,
216 'CE': 0.45,
217 'CG': 0.45,
218 'N': -0.44,
219 'NZ': -1.08,
220 'O': -0.68},
221 'MET': {'C': -0.54,
222 'CA': 0.02,
223 'CB': 0.45,
224 'CE': 0.63,
225 'CG': 0.45,
226 'N': -0.44,
227 'O': -0.68,
228 'SD': -0.30},
229 'PCA': {'C': -0.54,
230 'CA': 0.02,
231 'CB': 0.45,
232 'CD': -0.54,
233 'CG': 0.45,
234 'N': 1.52,
235 'O': -0.68,
236 'OE': -0.68},
237 'PHE': {'C': -0.54,
238 'CA': 0.02,
239 'CB': 0.45,
240 'CD1': 0.31,
241 'CD2': 0.31,
242 'CE1': 0.31,
243 'CE2': 0.31,
244 'CG': 0.09,
245 'CZ': 0.31,
246 'N': -0.44,
247 'O': -0.68},
248 'PRO': {'C': -0.54,
249 'CA': 0.02,
250 'CB': 0.45,
251 'CD': 0.45,
252 'CG': 0.45,
253 'N': -0.92,
254 'O': -0.68},
255 'SER': {'C': -0.54,
256 'CA': 0.02,
257 'CB': 0.45,
258 'N': -0.44,
259 'O': -0.68,
260 'OG': -0.99},
261 'THR': {'C': -0.54,
262 'CA': 0.02,
263 'CB': 0.02,
264 'CG2': 0.63,
265 'N': -0.44,
266 'O': -0.68,
267 'OG1': -0.93},
268 'TRP': {'C': -0.54,
269 'CA': 0.02,
270 'CB': 0.45,
271 'CD1': 0.31,
272 'CD2': 0.24,
273 'CE2': 0.24,
274 'CE3': 0.31,
275 'CG': 0.09,
276 'CH2': 0.31,
277 'CZ2': 0.31,
278 'CZ3': 0.31,
279 'N': -0.44,
280 'NE1': -0.55,
281 'O': -0.68},
282 'TYR': {'C': -0.54,
283 'CA': 0.02,
284 'CB': 0.45,
285 'CD1': 0.31,
286 'CD2': 0.31,
287 'CE1': 0.31,
288 'CE2': 0.31,
289 'CG': 0.09,
290 'CZ': 0.09,
291 'N': -0.44,
292 'O': -0.68,
293 'OH': -0.17},
294 'VAL': {'C': -0.54,
295 'CA': 0.02,
296 'CB': 0.02,
297 'CG1': 0.63,
298 'CG2': 0.63,
299 'N': -0.44,
300 'O': -0.68}}
301
302def assignfi(fidata, atoms):
303 """assign fi parameters to each atom in the pdbfile"""
304 n = len(atoms)
305 from numpy import empty, float32
306 fi = empty((n,), float32)
307 resname = atoms.residues.names
308 aname = atoms.names
309 for i in range(n):
310 rname = resname[i]
311 rfidata = fidata.get(rname)
312 if rfidata:
313 fi[i]=rfidata.get(aname[i], 0)
314 return fi
315
316def _griddimcalc(listcoord, spacing, gridmargin):
317 """Determination of the grid dimension"""
318 coordmin = min(listcoord) - gridmargin
319 coordmax = max(listcoord) + gridmargin
320 adjustment = ((spacing - (coordmax - coordmin)) % spacing) / 2.
321 coordmin = coordmin - adjustment
322 coordmax = coordmax + adjustment
323 ngrid = int(round((coordmax - coordmin) / spacing))
324 return coordmin, coordmax, ngrid
325
326def calculatefimap(atoms, method, spacing, max_dist, nexp):
327 """Calculation loop"""
328
329 #grid settings in angstrom
330 gridmargin = Defaults().gridmargin
331 xyz = atoms.scene_coords
332 xmingrid, xmaxgrid, nxgrid = _griddimcalc(xyz[:,0], spacing, gridmargin)
333 ymingrid, ymaxgrid, nygrid = _griddimcalc(xyz[:,1], spacing, gridmargin)
334 zmingrid, zmaxgrid, nzgrid = _griddimcalc(xyz[:,2], spacing, gridmargin)
335 bounds = [[xmingrid, xmaxgrid],
336 [ymingrid, ymaxgrid],
337 [zmingrid, zmaxgrid]]
338 origin = (xmingrid, ymingrid, zmingrid)
339
340 fi_table = Defaults().fidatadefault
341 fi = assignfi(fi_table, atoms)
342
343 from numpy import zeros, float32
344 pot = zeros((nzgrid+1, nygrid+1, nxgrid+1), float32)
345 from ._mlp import mlp_sum
346 mlp_sum(xyz, fi, origin, spacing, max_dist, method, nexp, pot)
347
348 return pot, bounds
349
350def mlp_sum(xyz, fi, origin, spacing, max_dist, method, nexp, pot):
351 if method == 'dubost':
352 computemethod = _dubost
353 elif method == 'fauchere':
354 computemethod = _fauchere
355 elif method == 'brasseur':
356 computemethod = _brasseur
357 elif method == 'buckingham':
358 computemethod = _buckingham
359 elif method == 'type5':
360 computemethod = _type5
361 else:
362 raise ValueError('Unknown lipophilicity method %s\n' % computemethod)
363
364 from numpy import zeros, float32, empty, subtract, sqrt
365 grid_pt = empty((3,), float32)
366 dxyz = empty((len(xyz),3), float32)
367 dist = empty((len(xyz),), float32)
368 nz,ny,nx = pot.shape
369 x0,y0,z0 = origin
370 for k in range(nz):
371 grid_pt[2] = z0 + k * spacing
372 for j in range(ny):
373 grid_pt[1] = y0 + j * spacing
374 for i in range(nx):
375 #Evaluation of the distance between th grid point and each atoms
376 grid_pt[0] = x0 + i * spacing
377 subtract(xyz, grid_pt, dxyz)
378 dxyz *= dxyz
379 dist = dxyz[:,0]
380 dist += dxyz[:,1]
381 dist += dxyz[:,2]
382 sqrt(dist, dist)
383 pot[k,j,i] = computemethod(fi, dist, nexp)
384
385def _dubost(fi, d, n):
386 return (100 * fi / (1 + d)).sum()
387
388def _fauchere(fi, d, n):
389 from numpy import exp
390 return (100 * fi * exp(-d)).sum()
391
392def _brasseur(fi, d, n):
393 #3.1 division is there to remove any units in the equation
394 #3.1A is the average diameter of a water molecule (2.82 -> 3.2)
395 from numpy import exp
396 return (100 * fi * exp(-d/3.1)).sum()
397
398def _buckingham(fi, d, n):
399 return (100 * fi / (d**n)).sum()
400
401def _type5(fi, d, n):
402 from numpy import exp, sqrt
403 return (100 * fi * exp(-sqrt(d))).sum()