''' MSMS module Author: Tobias Schmidt This module is for calculating MSMS surfaces as well as surface areas (SESA, SASA) from OpenStructure using the external program MSMS. How To Use This Module: 1. Import it (e.g. as "from ost.bindings import msms") 2. Use it (e.g. as "surfaces_list = msms.CalculateSurface(entity)" "(sesa,sasa) = msms.CalculateSurfaceArea(entity)") Requirement: - MSMS installed ''' import tempfile import subprocess import os from ost import io from ost import mol from ost import settings from ost import geom class MsmsProcessError(Exception): """ Python 2.4 and older do not include the CalledProcessError exception. This class substitutes it. """ def __init__(self, returncode,command): self.returncode = returncode self.command = command def __str__(self): return repr(self.returncode) def GetVersion(msms_exe=None, msms_env=None): """ Get version of MSMS executable """ msms_executable = _GetExecutable(msms_exe, msms_env) command = "%s" % (msms_executable) proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE) stdout_value, stderr_value = proc.communicate() version = "" for l in stdout_value.decode().splitlines(): if l[0:4]=='MSMS': version = l.split(' ')[1] return version if version=="": LogWarning('Could not parse MSMS version string') return def _GetExecutable(msms_exe, msms_env): """ Function to check if MSMS executable is present :param msms_exe: Explicit path to msms executable :param msms_env: Environment variable pointing to msms executable :returns: Path to the executable :raises: :class:`~ost.FileNotFound` if executable is not found """ return settings.Locate('msms', explicit_file_name=msms_exe, env_name=msms_env) def _SetupFiles(entity, selection): """ Setup files for MSMS calculation in temporary directory :param entity: The entity for which the surface is to be calculated :type entity: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityHandle` :param selection: Calculate surface for subset of entity :type selection: :class:`str` :returns: tuple containing temporary directory and msms input file :raises: :class:`RuntimeError` if selection is not valid """ # create temporary directory tmp_dir_name=tempfile.mkdtemp() # select only heavy atoms if no_hydrogens is true entity_view=entity.Select(selection) if not entity_view.IsValid(): raise RuntimeError("Could not create view for selection (%s)"%(selection)) # write entity to tmp file tmp_file_name=os.path.join(tmp_dir_name,"entity") tmp_file_handle=open(tmp_file_name, 'w') for a in entity_view.GetAtomList(): position=a.GetPos() tmp_file_handle.write('%8.3f %8.3f %8.3f %4.2f\n' % (position[0], position[1], position[2], a.radius)) tmp_file_handle.close() return (tmp_dir_name, tmp_file_name) def _ParseAreaFile(entity, selection, file, asa_prop, esa_prop): """ Reads Area file (-af) and attach sasa and sesa per atom to an entitiy :param entity: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView` for attaching sasa and sesa on atom level :param file: Filename of area file :param asa_prop: Name of the float property for SASA :param esa_prop: Name of the float property for SESA :raises: :class:`RuntimeError` if number of atoms in file != number of atoms in entity """ view=entity.Select(selection) area_fh = open(file) area_lines = area_fh.readlines() area_fh.close() # shift first line area_lines = area_lines[1:] if view.GetAtomCount() != len(area_lines): raise RuntimeError("Atom count (%d) unequeal to number of atoms in area file (%d)" % (view.GetAtomCount(), len(area_lines))) for l in area_lines: atom_no, sesa, sasa = l.split() a = view.atoms[int(atom_no)] if asa_prop: a.SetFloatProp(asa_prop, float(sasa)) if esa_prop: a.SetFloatProp(esa_prop, float(sesa)) def _CleanupFiles(dir_name): """ Function which recursively deletes a directory and all the files contained in it. *Warning*: This method removes also non-empty directories without asking, so be careful! """ import shutil shutil.rmtree(dir_name) def _RunMSMS(command): """ Run the MSMS surface calculation This functions starts the external MSMS executable and returns the stdout of MSMS. :param command: Command to execute :returns: stdout of MSMS :raises: :class:`CalledProcessError` for non-zero return value """ proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE) stdout_value, stderr_value = proc.communicate() #check for successful completion of msms if proc.returncode!=0: print("WARNING: msms error\n", stdout_value.decode()) raise MsmsProcessError(proc.returncode, command) return stdout_value.decode() def CalculateSurfaceArea(entity, density=1.0, radius=1.5, all_surf=False, no_hydrogens=False, no_hetatoms=False, no_waters=False, selection='', msms_exe=None, msms_env=None, keep_files=False, attach_asa=None, attach_esa=None): """ Calculates analytical solvent excluded and solvent accessible surface area by using the external MSMS program. This method calculates the molecular surface areas by invoking the external program MSMS. First, it is checked if the MSMS executable is present, then, the necessary files are prepared in a temporary directory and MSMS is executed. The last step is to remove the temporary directory. :param entity: OST entity to calculate surface :param density: Surface point density :param radius: Surface probe radius :param all_surf: Calculate surface area for all cavities (returns multiple surfaces areas as a list) :param no_hydrogens: Calculate surface only for hevy atoms :param selection: Calculate surface for subset of entity :param msms_exe: msms executable (full path to executable) :param msms_env: msms environment variable :param keep_files: Do not delete temporary files :param attach_asa: Attaches per atom SASA to specified FloatProp at atom level :param attach_esa: Attaches per atom SESA to specified FloatProp at atom level :returns: Tuple of lists for (SES, SAS) """ import re # check if msms executable is specified msms_executable=_GetExecutable(msms_exe, msms_env) # parse selection if no_hydrogens: if selection!='': selection+=" and " selection+="ele!=H" if no_hetatoms: if selection!='': selection+=" and " selection+="ishetatm=False" if no_waters: if selection!='': selection+=" and " selection+="rname!=HOH" # setup files for msms (msms_data_dir, msms_data_file)=_SetupFiles(entity, selection) # set command line command="%s -if %s -of %s -density %s -probe_radius %s -surface ases" % \ (msms_executable, msms_data_file, msms_data_file, density, radius) if all_surf: command+=" -all" if attach_asa != None or attach_esa != None: command+=" -af %s" % os.path.join(msms_data_dir, "asa_atom") # run msms stdout_value=_RunMSMS(command) # add sesa and asa to entity if attach_asa is specified if attach_asa != None or attach_esa != None: _ParseAreaFile(entity, selection, os.path.join(msms_data_dir, "asa_atom.area"), attach_asa, attach_esa) # parse MSMS output msms_ases=[] msms_asas=[] data_paragraph=False for line in stdout_value.splitlines(): if re.match('MSMS terminated normally', line): data_paragraph=False if data_paragraph: (ses_,sas_)=line.split()[5:7] msms_ases.append(float(ses_)) msms_asas.append(float(sas_)) if re.match(' Comp. probe_radius, reent, toric, contact SES SAS', line): data_paragraph=True # clean up if not keep_files: _CleanupFiles(msms_data_dir) return (msms_ases, msms_asas) def CalculateSurfaceVolume(entity, density=1.0, radius=1.5, all_surf=False, no_hydrogens=False, no_hetatoms=False, no_waters=False, selection='', msms_exe=None, msms_env=None, keep_files=False, attach_asa=None, attach_esa=None): """ Calculates the volume of the solvent excluded surface by using the external MSMS program. This method calculates the volume of the molecular surface by invoking the external program MSMS. First, it is checked if the MSMS executable is present, then, the necessary files are prepared in a temporary directory and MSMS is executed. The last step is to remove the temporary directory. :param entity: OST entity to calculate surface :param density: Surface point density :param radius: Surface probe radius :param all_surf: Calculate surface area for all cavities (returns multiple surfaces areas as a list) :param no_hydrogens: Calculate surface only for hevy atoms :param selection: Calculate surface for subset of entity :param msms_exe: msms executable (full path to executable) :param msms_env: msms environment variable :param keep_files: Do not delete temporary files :param attach_asa: Attaches per atom SASA to specified FloatProp at atom level :param attach_esa: Attaches per atom SESA to specified FloatProp at atom level :returns: Tuple of lists for (SES, SAS) """ import re # check if msms executable is specified msms_executable=_GetExecutable(msms_exe, msms_env) # parse selection if no_hydrogens: if selection!='': selection+=" and " selection+="ele!=H" if no_hetatoms: if selection!='': selection+=" and " selection+="ishetatm=False" if no_waters: if selection!='': selection+=" and " selection+="rname!=HOH" # setup files for msms (msms_data_dir, msms_data_file)=_SetupFiles(entity, selection) # set command line command="%s -if %s -of %s -density %s -probe_radius %s " % \ (msms_executable, msms_data_file, msms_data_file, density, radius) if all_surf: command+=" -all" if attach_asa != None or attach_esa != None: command+=" -af %s" % os.path.join(msms_data_dir, "asa_atom") # run msms stdout_value=_RunMSMS(command) # add sesa and asa to entity if attach_asa is specified if attach_asa != None or attach_esa != None: _ParseAreaFile(entity, selection, os.path.join(msms_data_dir, "asa_atom.area"), attach_asa, attach_esa) # parse MSMS output ses_volume=0 for line in stdout_value.splitlines(): if re.match(' Total ses_volume:', line): ses_volume=float(line.split(':')[1]) # clean up if not keep_files: _CleanupFiles(msms_data_dir) return ses_volume def CalculateSurface(entity, density=1.0, radius=1.5, all_surf=False, no_hydrogens=False, no_hetatoms=False, no_waters=False, selection='', msms_exe=None, msms_env=None, keep_files=False, attach_asa=None, attach_esa=None): """ Calculates molecular surface by using the external MSMS program This method calculates a molecular surface by invoking the external program MSMS. First, it is checked if the MSMS executable is present, then, the necessary files are prepared in a temporary directory and MSMS is executed. The last step is to remove the temporary directory. :param entity: Entity for which the surface is to be calculated :param density: Surface point density :param radius: Surface probe radius :param all_surf: Calculate surface for all cavities (returns multiple surfaces as a list) :param no_hydrogens: Calculate surface only for heavy atoms :param selection: Calculate surface for subset of entity :param msms_exe: msms executable (full path to executable) :param msms_env: msms environment variable :param keep_files: Do not delete temporary files :param attach_asa: Attaches per atom SASA to specified FloatProp at atom level :param attach_esa: Attaches per atom SESA to specified FloatProp at atom level :returns: list of :class:`~ost.mol.SurfaceHandle` objects """ import os import re # check if msms executable is specified msms_executable=_GetExecutable(msms_exe, msms_env) # parse selection if no_hydrogens: if selection!='': selection+=" and " selection+="ele!=H" if no_hetatoms: if selection!='': selection+=" and " selection+="ishetatm=False" if no_waters: if selection!='': selection+=" and " selection+="rname!=HOH" # setup files for msms (msms_data_dir, msms_data_file)=_SetupFiles(entity, selection) # set command line command="%s -if %s -of %s -density %s -probe_radius %s" % (msms_executable, msms_data_file, msms_data_file, density, radius) if all_surf: command+=" -all" if attach_asa != None or attach_esa != None: command+=" -af %s" % os.path.join(msms_data_dir, "asa_atom") # run msms stdout_value=_RunMSMS(command) # add sesa and asa to entity if attach_asa is specified if attach_asa != None or attach_esa != None: _ParseAreaFile(entity, selection, os.path.join(msms_data_dir, "asa_atom.area"), attach_asa, attach_esa) # parse msms output num_surf=0 for line in stdout_value.splitlines(): if re.search('RS component [0-9]+ identified',line): num_surf=int(line.split()[2]) # get surfaces entity_sel = entity.Select(selection) msms_surfaces=[] s = io.LoadSurface(msms_data_file, "msms") s.Attach(entity_sel, 3+radius) msms_surfaces.append(s) for n in range(1,num_surf+1): filename=msms_data_file+'_'+str(n) s = io.LoadSurface(filename, "msms") s.Attach(entity_sel, 3+radius) msms_surfaces.append(s) # clean up if not keep_files: _CleanupFiles(msms_data_dir) return msms_surfaces