Commit naccess parser for surface area calculation

git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2625 5a81b35b-ba03-0410-adc8-b2c5c5119f08

modules/bindings/pymod/naccess.py

+'''
+Naccess module
+
+Author: Florian Kiefer
+
+This module is for calculating surface areas
+from OpenStructure using the external program naccess
+
+How To Use This Module:
+ 1. Import it (e.g. as "from ost.bindings import naccess")
+ 2. Use it (e.g. as "sasa = naccess.CalculateSurfaceArea(entity)")
+
+Requirement:
+ - naccess installed
+'''
+
+import tempfile
+import subprocess
+import os
+import re
+from ost import io
+from ost import mol
+from ost import settings
+from ost import geom
+
+## \brief Method to check if naccess executable is present
+#
+# \param naccess Explicit path to msms executable
+# \return         Path to the executable
+# \exception      FileNotFound if executable is not found
+def _GetExecutable(naccess_exe):
+  return settings.Locate('naccess', explicit_file_name=naccess_exe)
+
+## \brief Setup files for naccess calculation in temporary directory
+#
+# \param entity      EntityHandle or EntityView to calculate surface
+# \param selection   Calculate surface for subset of entity 
+# \return            array containing temporary directory, input filename for naccess and directory of the input file
+# \exception         RuntimeError if selection is not valid
+def _SetupFiles(entity, selection):
+  # create temporary directory
+  tmp_dir_name=tempfile.mkdtemp()
+
+  # select only heavy atoms if no_hydrogens is true
+  entity_view=entity.Select(selection)
+  if len(entity_view.atoms) > 50000:
+    raise RuntimeError, "Too much atoms for NACCESS (> 50 000)"
+  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.pdb")
+  tmp_file_base = os.path.join(tmp_dir_name,"entity")
+  io.SavePDB(entity_view, tmp_file_name)
+  return (tmp_dir_name, tmp_file_name, tmp_file_base)
+
+## \brief Reads Area file (.asa) and attach asa per atom to an entitiy
+#
+# \param entity   EntityHandle or EntityView for attaching sasa on atom level
+# \param file     Filename of area file
+# \param asa_atom Name of the float property for SASA
+ 
+
+def _ParseAsaFile(entity, file, asa_atom):
+  
+  asa_fh = open(file)
+  asa_lines = asa_fh.readlines()
+  asa_fh.close()
+  
+  for l in  asa_lines:
+    if l.startswith("ATOM"):
+      # get res_number, chain_id and atom name
+      atom_name = l[12:16]
+      chain_id = l[21]
+      res_number = l[22:27]
+      asa = l[54:63]
+      atom_name = atom_name.strip()
+      chain_id = chain_id.strip()
+      res_number = res_number.strip()
+      asa = asa.strip()
+      #print "res_number:", res_number
+      m=re.match(r'(?P<num>-?\d+)(?P<ins>\w)?', res_number)
+      di = m.groupdict()
+      
+      if di["ins"] == None:
+        resNum = mol.ResNum(int(di["num"]))
+      else:
+        resNum = mol.ResNum(int(di["num"]), di["ins"])
+      #print res_number, resNum.num, resNum.ins
+      a = entity.FindAtom(chain_id, resNum, atom_name)
+      a.SetFloatProp(asa_atom, float(asa))
+      
+## \brief Reads Area file (.rsa) and attach asa (absolute + relative) per residue to an entitiy
+#
+# \param entity   EntityHandle or EntityView for attaching sasa on atom level
+# \param file     Filename of .rsa file
+# \param asa_atom Name of the float property for absolute SASA
+# \param asa_atom Name of the float property for relative SASA
+# \exception         RuntimeError if residue names are not the same  
+
+def _ParseRsaFile(enti,file, asa_abs, asa_rel):
+  area_fh = open(file)
+  area_lines = area_fh.readlines()
+  area_fh.close()
+  # shift first line
+  area_lines = area_lines[4:]
+  
+  
+  for l in area_lines:
+    if l.startswith("RES"):
+      p = re.compile(r'\s+')
+      t = p.split(l)
+      #res_name, chain_id , res_number, abs_all, rel_all = t[1:6]
+      res_name = l[3:8]
+      res_name = res_name.strip()
+      chain_id = l[8:9]
+      res_number = l[9:14]
+      res_number= res_number.strip()
+      #print l[15:30]
+      abs_all, rel_all =l[15:28].strip().split()
+      m=re.match(r'(?P<num>-?\d+)(?P<ins>\w)?', res_number)
+      di = m.groupdict()
+      if di["ins"] == None:
+        resNum = mol.ResNum(int(di["num"]))
+      else:
+        resNum = mol.ResNum(int(di["num"]), di["ins"])
+      
+      res = enti.handle.FindResidue(chain_id, resNum)
+      #res = entity.FindResidue(chain_id, mol.ResNum(int(res_number)))
+      #print res_name, chain_id, res_number
+      if res_name == res.name:
+        res.SetFloatProp(asa_rel, float(rel_all) )
+        res.SetFloatProp(asa_abs, float(abs_all) )
+      else:
+        raise RuntimeError, "Residue Names are not the same for ResNumb: %s (%s vs %s)" % (res_number, res.name, res_name)
+      
+## \brief Method which recursively deletes a directory
+#
+# \warning This method removes also non-empty directories without asking, so
+#          be careful!
+def __CleanupFiles(dir_name):
+  import shutil
+  shutil.rmtree(dir_name)
+
+## \brief Method to run the MSMS surface calculation
+#
+# This method starts the external MSMS executable and returns the stdout of MSMS
+#
+# \param command          Command to execute
+# \return                 stdout of MSMS
+# \exception              CalledProcessError for non-zero return value
+def _RunNACCESS(command, temp_dir):
+  proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, cwd = temp_dir)
+  stdout_value, stderr_value = proc.communicate()
+
+  #check for successful completion of msms
+  if proc.returncode!=0:
+    print "WARNING: naccess error\n", stdout_value
+    raise subprocess.CalledProcessError(proc.returncode, command)
+
+  return stdout_value
+  
+
+## \brief Calculates analytical the solvent accessible surface
+#  area by using the external naccess program
+#
+# This method calculates the molecular surface areas by invoking the external
+# program naccess. First, it is checked if the naccess executable is present, then, 
+# the necessary files are prepared in a temporary directory and naccess is
+# executed. The last step is to remove the temporary directory.
+# 
+#
+# \param entity             OST entity to calculate surface
+# \param radius             Surface probe radius
+# \param include_hydrogens  Calculate surface including hydrogens
+# \param include_hetatm     Calculate surface including hetatms
+# \param include_water      Calculate surface including water
+# \param selection          Calculate surface for subset of entity
+# \param naccess  _exe      msms executable (full path to executable)
+# \param keep_files         Do not delete temporary files
+# \param asa_abs            Attaches per residue absolute SASA to specified FloatProp on residue level
+# \param asa_rel            Attaches per residue relative SASA to specified FloatProp on residue level
+# \param asa_atom           Attaches per atom SASA to specified FloatProp at atom level
+# \return                   absolute SASA calculated using asa_atom 
+def CalculateSurfaceArea(entity,  radius=1.4,  
+                         include_hydrogens=False, include_hetatm = False, 
+                         include_water = False, selection="",
+                         naccess_exe=None, keep_files=False , asa_abs= "asaAbs", asa_rel="asaRel", asa_atom="asaAtom"):
+  import re 
+  
+  # check if msms executable is specified
+  naccess_executable=_GetExecutable(naccess_exe)
+  # parse selection
+  
+  # setup files for msms
+  (naccess_data_dir, naccess_data_file,naccess_data_base )=_SetupFiles(entity, selection)
+
+  # set command line
+  command="%s %s -p %f " % \
+          (naccess_executable, naccess_data_file, radius)
+  if include_hydrogens:
+    command = "%s -w" % command
+  if include_water:
+    command = "%s -y" % command
+  if include_hetatm:
+    command = "%s -h" % command
+  #print command
+  stdout_value=_RunNACCESS(command, naccess_data_dir)
+  
+  new_asa= os.path.join(naccess_data_dir, "%s.asa" % naccess_data_base) 
+  _ParseAsaFile(entity, new_asa, asa_atom)
+    
+  new_rsa = os.path.join(naccess_data_dir, "%s.rsa" % naccess_data_base) 
+  _ParseRsaFile(entity, new_rsa, asa_abs, asa_rel)
+  
+  # Calculate Asa for all atoms:
+  sasa = 0.0 
+  for a in entity.atoms:
+    sasa += a.GetFloatProp(asa_atom, 0.0)
+    
+  
+  # first read in result_file
+  
+  # parse MSMS output
+  
+  # clean up
+  if not keep_files:
+    __CleanupFiles(naccess_data_dir)
+
+  return sasa
+