OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
msms.py
Go to the documentation of this file.
1 '''
2 MSMS module
3 
4 Author: Tobias Schmidt
5 
6 This module is for calculating MSMS surfaces as well as surface areas
7 (SESA, SASA) from OpenStructure using the external program MSMS.
8 
9 How To Use This Module:
10  1. Import it (e.g. as "from ost.bindings import msms")
11  2. Use it (e.g. as "surfaces_list = msms.CalculateSurface(entity)"
12  "(sesa,sasa) = msms.CalculateSurfaceArea(entity)")
13 
14 Requirement:
15  - MSMS installed
16 '''
17 
18 import tempfile
19 import subprocess
20 import os
21 
22 from ost import io
23 from ost import mol
24 from ost import settings
25 from ost import geom
26 
27 
28 class MsmsProcessError(Exception):
29  """
30  Python 2.4 and older do not include the CalledProcessError exception. This
31  class substitutes it.
32  """
33  def __init__(self, returncode,command):
34  self.returncode = returncode
35  self.command = command
36  def __str__(self):
37  return repr(self.returncode)
38 
39 
40 def GetVersion(msms_exe=None, msms_env=None):
41  """
42  Get version of MSMS executable
43  """
44  msms_executable = _GetExecutable(msms_exe, msms_env)
45  command = "%s" % (msms_executable)
46  proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE)
47  stdout_value, stderr_value = proc.communicate()
48 
49  version = ""
50  for l in stdout_value.decode().splitlines():
51  if l[0:4]=='MSMS':
52  version = l.split(' ')[1]
53  return version
54  if version=="":
55  LogWarning('Could not parse MSMS version string')
56  return
57 
58 
59 def _GetExecutable(msms_exe, msms_env):
60  """
61  Function to check if MSMS executable is present
62 
63  :param msms_exe: Explicit path to msms executable
64  :param msms_env: Environment variable pointing to msms executable
65  :returns: Path to the executable
66  :raises: :class:`~ost.FileNotFound` if executable is not found
67  """
68  return settings.Locate('msms', explicit_file_name=msms_exe,
69  env_name=msms_env)
70 
71 
72 def _SetupFiles(entity, selection):
73  """
74  Setup files for MSMS calculation in temporary directory
75 
76  :param entity: The entity for which the surface is to be calculated
77  :type entity: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityHandle`
78  :param selection: Calculate surface for subset of entity
79  :type selection: :class:`str`
80  :returns: tuple containing temporary directory and msms input file
81  :raises: :class:`RuntimeError` if selection is not valid
82  """
83  # create temporary directory
84  tmp_dir_name=tempfile.mkdtemp()
85 
86  # select only heavy atoms if no_hydrogens is true
87  entity_view=entity.Select(selection)
88  if not entity_view.IsValid():
89  raise RuntimeError("Could not create view for selection (%s)"%(selection))
90 
91  # write entity to tmp file
92  tmp_file_name=os.path.join(tmp_dir_name,"entity")
93  tmp_file_handle=open(tmp_file_name, 'w')
94  for a in entity_view.GetAtomList():
95  position=a.GetPos()
96  tmp_file_handle.write('%8.3f %8.3f %8.3f %4.2f\n' % (position[0],
97  position[1], position[2], a.radius))
98  tmp_file_handle.close()
99 
100  return (tmp_dir_name, tmp_file_name)
101 
102 
103 def _ParseAreaFile(entity, selection, file, asa_prop, esa_prop):
104  """
105  Reads Area file (-af) and attach sasa and sesa per atom to an entitiy
106 
107  :param entity: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
108  for attaching sasa and sesa on atom level
109  :param file: Filename of area file
110  :param asa_prop: Name of the float property for SASA
111  :param esa_prop: Name of the float property for SESA
112  :raises: :class:`RuntimeError` if number of atoms in file != number of atoms in entity
113  """
114  view=entity.Select(selection)
115  area_fh = open(file)
116  area_lines = area_fh.readlines()
117  area_fh.close()
118  # shift first line
119  area_lines = area_lines[1:]
120  if view.GetAtomCount() != len(area_lines):
121  raise RuntimeError("Atom count (%d) unequeal to number of atoms in area file (%d)" % (view.GetAtomCount(), len(area_lines)))
122  for l in area_lines:
123  atom_no, sesa, sasa = l.split()
124  a = view.atoms[int(atom_no)]
125  if asa_prop:
126  a.SetFloatProp(asa_prop, float(sasa))
127  if esa_prop:
128  a.SetFloatProp(esa_prop, float(sesa))
129 
130 
131 
132 def _CleanupFiles(dir_name):
133  """
134  Function which recursively deletes a directory and all the files contained
135  in it. *Warning*: This method removes also non-empty directories without
136  asking, so be careful!
137  """
138  import shutil
139  shutil.rmtree(dir_name)
140 
141 def _RunMSMS(command):
142  """
143  Run the MSMS surface calculation
144 
145  This functions starts the external MSMS executable and returns the stdout of
146  MSMS.
147 
148  :param command: Command to execute
149  :returns: stdout of MSMS
150  :raises: :class:`CalledProcessError` for non-zero return value
151  """
152  proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE)
153  stdout_value, stderr_value = proc.communicate()
154 
155  #check for successful completion of msms
156  if proc.returncode!=0:
157  print("WARNING: msms error\n", stdout_value.decode())
158  raise MsmsProcessError(proc.returncode, command)
159 
160  return stdout_value.decode()
161 
162 
163 
164 def CalculateSurfaceArea(entity, density=1.0, radius=1.5, all_surf=False,
165  no_hydrogens=False, no_hetatoms=False, no_waters=False,
166  selection='',
167  msms_exe=None, msms_env=None, keep_files=False,
168  attach_asa=None, attach_esa=None):
169  """
170  Calculates analytical solvent excluded and solvent accessible surface
171  area by using the external MSMS program.
172 
173  This method calculates the molecular surface areas by invoking the external
174  program MSMS. First, it is checked if the MSMS executable is present, then,
175  the necessary files are prepared in a temporary directory and MSMS is
176  executed. The last step is to remove the temporary directory.
177 
178 
179  :param entity: OST entity to calculate surface
180  :param density: Surface point density
181  :param radius: Surface probe radius
182  :param all_surf: Calculate surface area for all cavities (returns multiple
183  surfaces areas as a list)
184  :param no_hydrogens: Calculate surface only for hevy atoms
185  :param selection: Calculate surface for subset of entity
186  :param msms_exe: msms executable (full path to executable)
187  :param msms_env: msms environment variable
188  :param keep_files: Do not delete temporary files
189  :param attach_asa: Attaches per atom SASA to specified FloatProp at atom level
190  :param attach_esa: Attaches per atom SESA to specified FloatProp at atom level
191  :returns: Tuple of lists for (SES, SAS)
192  """
193  import re
194 
195  # check if msms executable is specified
196  msms_executable=_GetExecutable(msms_exe, msms_env)
197 
198  # parse selection
199  if no_hydrogens:
200  if selection!='':
201  selection+=" and "
202  selection+="ele!=H"
203 
204  if no_hetatoms:
205  if selection!='':
206  selection+=" and "
207  selection+="ishetatm=False"
208 
209  if no_waters:
210  if selection!='':
211  selection+=" and "
212  selection+="rname!=HOH"
213 
214  # setup files for msms
215  (msms_data_dir, msms_data_file)=_SetupFiles(entity, selection)
216 
217  # set command line
218  command="%s -if %s -of %s -density %s -probe_radius %s -surface ases" % \
219  (msms_executable, msms_data_file, msms_data_file, density, radius)
220  if all_surf:
221  command+=" -all"
222  if attach_asa != None or attach_esa != None:
223  command+=" -af %s" % os.path.join(msms_data_dir, "asa_atom")
224  # run msms
225  stdout_value=_RunMSMS(command)
226 
227  # add sesa and asa to entity if attach_asa is specified
228  if attach_asa != None or attach_esa != None:
229  _ParseAreaFile(entity, selection, os.path.join(msms_data_dir, "asa_atom.area"),
230  attach_asa, attach_esa)
231 
232  # parse MSMS output
233  msms_ases=[]
234  msms_asas=[]
235  data_paragraph=False
236  for line in stdout_value.splitlines():
237  if re.match('MSMS terminated normally', line):
238  data_paragraph=False
239  if data_paragraph:
240  (ses_,sas_)=line.split()[5:7]
241  msms_ases.append(float(ses_))
242  msms_asas.append(float(sas_))
243  if re.match(' Comp. probe_radius, reent, toric, contact SES SAS', line):
244  data_paragraph=True
245 
246  # clean up
247  if not keep_files:
248  _CleanupFiles(msms_data_dir)
249 
250  return (msms_ases, msms_asas)
251 
252 def CalculateSurfaceVolume(entity, density=1.0, radius=1.5, all_surf=False,
253  no_hydrogens=False, no_hetatoms=False, no_waters=False,
254  selection='',
255  msms_exe=None, msms_env=None, keep_files=False,
256  attach_asa=None, attach_esa=None):
257  """
258  Calculates the volume of the solvent excluded surface by using the external MSMS program.
259 
260  This method calculates the volume of the molecular surface by invoking the external
261  program MSMS. First, it is checked if the MSMS executable is present, then,
262  the necessary files are prepared in a temporary directory and MSMS is
263  executed. The last step is to remove the temporary directory.
264 
265 
266  :param entity: OST entity to calculate surface
267  :param density: Surface point density
268  :param radius: Surface probe radius
269  :param all_surf: Calculate surface area for all cavities (returns multiple
270  surfaces areas as a list)
271  :param no_hydrogens: Calculate surface only for hevy atoms
272  :param selection: Calculate surface for subset of entity
273  :param msms_exe: msms executable (full path to executable)
274  :param msms_env: msms environment variable
275  :param keep_files: Do not delete temporary files
276  :param attach_asa: Attaches per atom SASA to specified FloatProp at atom level
277  :param attach_esa: Attaches per atom SESA to specified FloatProp at atom level
278  :returns: Tuple of lists for (SES, SAS)
279  """
280  import re
281 
282  # check if msms executable is specified
283  msms_executable=_GetExecutable(msms_exe, msms_env)
284 
285  # parse selection
286  if no_hydrogens:
287  if selection!='':
288  selection+=" and "
289  selection+="ele!=H"
290 
291  if no_hetatoms:
292  if selection!='':
293  selection+=" and "
294  selection+="ishetatm=False"
295 
296  if no_waters:
297  if selection!='':
298  selection+=" and "
299  selection+="rname!=HOH"
300 
301  # setup files for msms
302  (msms_data_dir, msms_data_file)=_SetupFiles(entity, selection)
303 
304  # set command line
305  command="%s -if %s -of %s -density %s -probe_radius %s " % \
306  (msms_executable, msms_data_file, msms_data_file, density, radius)
307  if all_surf:
308  command+=" -all"
309  if attach_asa != None or attach_esa != None:
310  command+=" -af %s" % os.path.join(msms_data_dir, "asa_atom")
311  # run msms
312  stdout_value=_RunMSMS(command)
313 
314  # add sesa and asa to entity if attach_asa is specified
315  if attach_asa != None or attach_esa != None:
316  _ParseAreaFile(entity, selection, os.path.join(msms_data_dir, "asa_atom.area"),
317  attach_asa, attach_esa)
318 
319  # parse MSMS output
320  ses_volume=0
321  for line in stdout_value.splitlines():
322  if re.match(' Total ses_volume:', line):
323  ses_volume=float(line.split(':')[1])
324 
325  # clean up
326  if not keep_files:
327  _CleanupFiles(msms_data_dir)
328 
329  return ses_volume
330 
331 
332 def CalculateSurface(entity, density=1.0, radius=1.5, all_surf=False,
333  no_hydrogens=False, no_hetatoms=False, no_waters=False,
334  selection='',
335  msms_exe=None, msms_env=None, keep_files=False,
336  attach_asa=None, attach_esa=None):
337 
338  """
339  Calculates molecular surface by using the external MSMS program
340 
341  This method calculates a molecular surface by invoking the external program
342  MSMS. First, it is checked if the MSMS executable is present, then, the
343  necessary files are prepared in a temporary directory and MSMS is executed.
344  The last step is to remove the temporary directory.
345 
346 
347  :param entity: Entity for which the surface is to be calculated
348  :param density: Surface point density
349  :param radius: Surface probe radius
350  :param all_surf: Calculate surface for all cavities (returns multiple
351  surfaces as a list)
352  :param no_hydrogens: Calculate surface only for heavy atoms
353  :param selection: Calculate surface for subset of entity
354  :param msms_exe: msms executable (full path to executable)
355  :param msms_env: msms environment variable
356  :param keep_files: Do not delete temporary files
357  :param attach_asa: Attaches per atom SASA to specified FloatProp at atom level
358  :param attach_esa: Attaches per atom SESA to specified FloatProp at atom level
359  :returns: list of :class:`~ost.mol.SurfaceHandle` objects
360  """
361  import os
362  import re
363 
364  # check if msms executable is specified
365  msms_executable=_GetExecutable(msms_exe, msms_env)
366 
367  # parse selection
368  if no_hydrogens:
369  if selection!='':
370  selection+=" and "
371  selection+="ele!=H"
372 
373  if no_hetatoms:
374  if selection!='':
375  selection+=" and "
376  selection+="ishetatm=False"
377 
378  if no_waters:
379  if selection!='':
380  selection+=" and "
381  selection+="rname!=HOH"
382 
383  # setup files for msms
384  (msms_data_dir, msms_data_file)=_SetupFiles(entity, selection)
385 
386  # set command line
387  command="%s -if %s -of %s -density %s -probe_radius %s" % (msms_executable,
388  msms_data_file, msms_data_file, density, radius)
389  if all_surf:
390  command+=" -all"
391  if attach_asa != None or attach_esa != None:
392  command+=" -af %s" % os.path.join(msms_data_dir, "asa_atom")
393 
394  # run msms
395  stdout_value=_RunMSMS(command)
396 
397  # add sesa and asa to entity if attach_asa is specified
398  if attach_asa != None or attach_esa != None:
399  _ParseAreaFile(entity, selection, os.path.join(msms_data_dir, "asa_atom.area"),
400  attach_asa, attach_esa)
401 
402  # parse msms output
403  num_surf=0
404  for line in stdout_value.splitlines():
405  if re.search('RS component [0-9]+ identified',line):
406  num_surf=int(line.split()[2])
407 
408  # get surfaces
409  entity_sel = entity.Select(selection)
410  msms_surfaces=[]
411  s = io.LoadSurface(msms_data_file, "msms")
412  s.Attach(entity_sel, 3+radius)
413  msms_surfaces.append(s)
414  for n in range(1,num_surf+1):
415  filename=msms_data_file+'_'+str(n)
416  s = io.LoadSurface(filename, "msms")
417  s.Attach(entity_sel, 3+radius)
418  msms_surfaces.append(s)
419 
420  # clean up
421  if not keep_files:
422  _CleanupFiles(msms_data_dir)
423 
424  return msms_surfaces
425 
def CalculateSurfaceArea
Definition: msms.py:168
def CalculateSurface
Definition: msms.py:336
def GetVersion
Definition: msms.py:40
def CalculateSurfaceVolume
Definition: msms.py:256