diff --git a/modules/bindings/pymod/naccess.py b/modules/bindings/pymod/naccess.py new file mode 100644 index 0000000000000000000000000000000000000000..f18d8e9a6029ffd824393e8eb00239a7452d1c6f --- /dev/null +++ b/modules/bindings/pymod/naccess.py @@ -0,0 +1,231 @@ +''' +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 +