Skip to content
Snippets Groups Projects
Commit 84f56a2b authored by Gerardo Tauriello's avatar Gerardo Tauriello
Browse files

SCHWED-2349: clean-up of Naccess binding.

parent dd1c5b6c
No related branches found
No related tags found
No related merge requests found
''' '''
Naccess module Naccess module
Author: Florian Kiefer Author: Florian Kiefer, Gerardo Tauriello (cleanup/speedup)
This module is for calculating surface areas This module is for calculating surface areas
from OpenStructure using the external program naccess from OpenStructure using the external program naccess
...@@ -18,16 +18,13 @@ import tempfile ...@@ -18,16 +18,13 @@ import tempfile
import subprocess import subprocess
import os import os
import re import re
from ost import io from ost import io, mol, settings, LogWarning
from ost import mol
from ost import settings
from ost import geom
def _GetExecutable(naccess_exe): def _GetExecutable(naccess_exe):
""" """
Method to check if naccess executable is present Method to check if naccess executable is present
:param naccess: Explicit path to msms executable :param naccess: Explicit path to naccess executable
:returns: Path to the executable :returns: Path to the executable
:exception: FileNotFound if executable is not found :exception: FileNotFound if executable is not found
""" """
...@@ -46,14 +43,14 @@ def _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms): ...@@ -46,14 +43,14 @@ def _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms):
""" """
# create temporary directory # create temporary directory
tmp_dir_name="" tmp_dir_name = ""
if scratch_dir!=None: if scratch_dir != None:
tmp_dir_name=tempfile.mkdtemp(dir=scratch_dir) tmp_dir_name = tempfile.mkdtemp(dir=scratch_dir)
else: else:
tmp_dir_name=tempfile.mkdtemp() tmp_dir_name = tempfile.mkdtemp()
# select only heavy atoms if no_hydrogens is true # select as specified by user
entity_view=entity.Select(selection) entity_view = entity.Select(selection)
if len(entity_view.atoms) > max_number_of_atoms: if len(entity_view.atoms) > max_number_of_atoms:
raise RuntimeError, "Too much atoms for NACCESS (> %s)" % max_number_of_atoms raise RuntimeError, "Too much atoms for NACCESS (> %s)" % max_number_of_atoms
if not entity_view.IsValid(): if not entity_view.IsValid():
...@@ -61,7 +58,7 @@ def _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms): ...@@ -61,7 +58,7 @@ def _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms):
# write entity to tmp file # write entity to tmp file
tmp_file_name = "entity.pdb" tmp_file_name = "entity.pdb"
tmp_file_base = os.path.join(tmp_dir_name,"entity") tmp_file_base = os.path.join(tmp_dir_name, "entity")
io.SavePDB(entity_view, os.path.join(tmp_dir_name, tmp_file_name)) io.SavePDB(entity_view, os.path.join(tmp_dir_name, tmp_file_name))
return (tmp_dir_name, tmp_file_name, tmp_file_base) return (tmp_dir_name, tmp_file_name, tmp_file_base)
...@@ -171,10 +168,11 @@ def _RunNACCESS(command, temp_dir): ...@@ -171,10 +168,11 @@ def _RunNACCESS(command, temp_dir):
:returns: stdout of MSMS :returns: stdout of MSMS
:exception: CalledProcessError for non-zero return value :exception: CalledProcessError for non-zero return value
""" """
proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, cwd = temp_dir) proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE,
cwd=temp_dir)
stdout_value, stderr_value = proc.communicate() stdout_value, stderr_value = proc.communicate()
#check for successful completion of msms # check for successful completion of naccess
if proc.returncode!=0: if proc.returncode!=0:
print "WARNING: naccess error\n", stdout_value print "WARNING: naccess error\n", stdout_value
raise subprocess.CalledProcessError(proc.returncode, command) raise subprocess.CalledProcessError(proc.returncode, command)
...@@ -184,8 +182,10 @@ def _RunNACCESS(command, temp_dir): ...@@ -184,8 +182,10 @@ def _RunNACCESS(command, temp_dir):
def CalculateSurfaceArea(entity, radius=1.4, def CalculateSurfaceArea(entity, radius=1.4,
include_hydrogens=False, include_hetatm = False, include_hydrogens=False, include_hetatm = False,
include_water = False, selection="", include_water=False, selection="", naccess_exe=None,
naccess_exe=None, keep_files=False , asa_abs= "asaAbs", asa_rel="asaRel", asa_atom="asaAtom", scratch_dir = None, max_number_of_atoms=50000): keep_files=False, asa_abs= "asaAbs", asa_rel="asaRel",
asa_atom="asaAtom", scratch_dir=None,
max_number_of_atoms=50000):
""" """
Calculates analytical the solvent accessible surface area by using the Calculates analytical the solvent accessible surface area by using the
external naccess program external naccess program
...@@ -203,56 +203,55 @@ def CalculateSurfaceArea(entity, radius=1.4, ...@@ -203,56 +203,55 @@ def CalculateSurfaceArea(entity, radius=1.4,
:param include_water: Calculate surface including water :param include_water: Calculate surface including water
:param selection: Calculate surface for subset of entity :param selection: Calculate surface for subset of entity
:param naccess_exe: naccess executable (full path to executable) :param naccess_exe: naccess executable (full path to executable)
:param keep_files: Do not delete temporary files :param keep_files: If True, do not delete temporary files
:param asa_abs: Attaches per residue absolute SASA to specified FloatProp on residue level :param asa_abs: Attaches per residue absolute SASA to specified
:param asa_rel: Attaches per residue relative SASA to specified FloatProp on residue level FloatProp on residue level
:param asa_atom: Attaches per atom SASA to specified FloatProp at atom level :param asa_rel: Attaches per residue relative SASA to specified
:param scratch_dir: Directory for temporary files (NACCESS is sensitive to "." in directory names FloatProp on residue level
:param max_number_of_atoms: Max Number of atoms in the entity (i.e. is limited in the default NACCESS version to 50 000) :param asa_atom: Attaches per atom SASA to specified FloatProp at
atom level
:param scratch_dir: Scratch directory. A subfolder for temporary files
is created in there. If not specified, a default
directory is used (see :func:`tempfile.mkdtemp`).
:param max_number_of_atoms: Max Number of atoms in the entity (i.e. is limited
in the default NACCESS version to 50 000)
:returns: absolute SASA calculated using asa_atom :returns: absolute SASA calculated using asa_atom
""" """
import re
# check if msms executable is specified # check if naccess executable is specified
naccess_executable=_GetExecutable(naccess_exe) naccess_executable=_GetExecutable(naccess_exe)
# parse selection
# setup files for msms # setup files for naccess
(naccess_data_dir, naccess_data_file,naccess_data_base )=_SetupFiles(entity, selection, scratch_dir, max_number_of_atoms) (naccess_data_dir, naccess_data_file, naccess_data_base) \
= _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms)
# set command line # set command line
command="%s %s -p %f " % \ command = "%s %s -p %f " % \
(naccess_executable, naccess_data_file, radius) (naccess_executable, naccess_data_file, radius)
if include_hydrogens: if include_hydrogens:
command = "%s -w" % command command = "%s -w" % command
if include_water: if include_water:
command = "%s -y" % command command = "%s -y" % command
if include_hetatm: if include_hetatm:
command = "%s -h" % command command = "%s -h" % command
#print command # execute naccess
stdout_value=_RunNACCESS(command, naccess_data_dir) stdout_value = _RunNACCESS(command, naccess_data_dir)
# parse outout
new_asa= os.path.join(naccess_data_dir, "%s.asa" % naccess_data_base) new_asa= os.path.join(naccess_data_dir, "%s.asa" % naccess_data_base)
_ParseAsaFile(entity, new_asa, asa_atom) _ParseAsaFile(entity, new_asa, asa_atom)
new_rsa = os.path.join(naccess_data_dir, "%s.rsa" % naccess_data_base) new_rsa = os.path.join(naccess_data_dir, "%s.rsa" % naccess_data_base)
_ParseRsaFile(entity, new_rsa, asa_abs, asa_rel) _ParseRsaFile(entity, new_rsa, asa_abs, asa_rel)
# Calculate Asa for all atoms: # sum up Asa for all atoms
sasa = 0.0 sasa = 0.0
for a in entity.atoms: for a in entity.atoms:
sasa += a.GetFloatProp(asa_atom, 0.0) sasa += a.GetFloatProp(asa_atom, 0.0)
# first read in result_file
# parse MSMS output
# clean up # clean up
if not keep_files: if not keep_files:
__CleanupFiles(naccess_data_dir) __CleanupFiles(naccess_data_dir)
return sasa return sasa
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment