diff --git a/modules/bindings/pymod/naccess.py b/modules/bindings/pymod/naccess.py index 1612c503f9a83e0284e021384b82834aefba166c..df56c2af0944ee603f1467ec61431fa735ed3a47 100644 --- a/modules/bindings/pymod/naccess.py +++ b/modules/bindings/pymod/naccess.py @@ -30,16 +30,35 @@ def _GetExecutable(naccess_exe): """ return settings.Locate('naccess', explicit_file_name=naccess_exe) +def _CheckNaccessRoot(naccess_root): + """ + :return: True, if given directory contains "accall" binary and files + "vdw.radii" and "standard.data". + :param naccess_root: Path to naccess folder to check. + """ + accall_exe = os.path.join(naccess_root, "accall") + check = (os.path.exists(accall_exe) and os.access(accall_exe, os.X_OK) \ + and os.path.exists(os.path.join(naccess_root, "vdw.radii")) \ + and os.path.exists(os.path.join(naccess_root, "standard.data"))) + if not check: + LogWarning("NACCESS: Could not find required files to launch accall " \ + "directly in %s." % naccess_root) + return check + def _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms): """ Setup files for naccess calculation in temporary directory - :param entity: EntityHandle or EntityView to calculate surface - :param selection: Calculate surface for subset of entity - :param scratch_dir: Directory for temporary files (NACCESS is sensitive to "." in directory names - :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: array containing temporary directory, input filename for naccess and directory of the input file - :exception: RuntimeError if selection is not valid + :param entity: OST entity to calculate surface + :param selection: Calculate surface for subset of entity + :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: Tuple containing temporary directory, input filename for naccess and + directory of the input file + :exception: RuntimeError if selection is not valid or too many atoms """ # create temporary directory @@ -103,7 +122,8 @@ def _ParseAsaFile(entity, file, asa_atom): if(a.IsValid()): a.SetFloatProp(asa_atom, float(asa)) else: - print chain_id, resNum, atom_name + LogWarning("NACCESS: invalid asa entry %s %s %s" \ + % (chain_id, resNum, atom_name)) def _ParseRsaFile(enti,file, asa_abs, asa_rel): """ @@ -161,33 +181,62 @@ def __CleanupFiles(dir_name): import shutil shutil.rmtree(dir_name) + +def _RunACCALL(command, temp_dir, query): + """ + Fast method to run the Naccess surface calculation. + + This method starts the accall binary directly and pipes in the input provided + in *query*. This is faster than calling the "naccess" script since the script + has a constant overhead of roughly 1.3s in each call. + + :param command: Command to execute + :param temp_dir: Command is executed with this working directory + :param query: User input to pipe into *command* + :returns: stdout of command + :exception: CalledProcessError for non-zero return value + """ + proc = subprocess.Popen(command, stdout=subprocess.PIPE, + stderr=subprocess.PIPE, stdin=subprocess.PIPE, + cwd=temp_dir) + stdout_value, stderr_value = proc.communicate(query) + + # check for successful completion of naccess + if proc.returncode != 0: + LogWarning("WARNING: naccess error\n%s\n%s" % (stdout_value, stderr_value)) + raise subprocess.CalledProcessError(proc.returncode, command) + + return stdout_value + + def _RunNACCESS(command, temp_dir): """ - Method to run the MSMS surface calculation + Method to run the Naccess surface calculation. - This method starts the external MSMS executable and returns the stdout of MSMS + This method starts the external Naccess executable and returns the stdout. - :param command: Command to execute - :returns: stdout of MSMS - :exception: CalledProcessError for non-zero return value + :param command: Command to execute + :param temp_dir: Command is executed with this working directory + :returns: stdout of command + :exception: CalledProcessError for non-zero return value """ proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, cwd=temp_dir) stdout_value, stderr_value = proc.communicate() # check for successful completion of naccess - if proc.returncode!=0: - print "WARNING: naccess error\n", stdout_value + if proc.returncode != 0: + LogWarning("WARNING: naccess error\n%s" % stdout_value) raise subprocess.CalledProcessError(proc.returncode, command) return stdout_value - + 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", scratch_dir=None, + naccess_root=None, 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 @@ -206,6 +255,9 @@ def CalculateSurfaceArea(entity, radius=1.4, :param include_water: Calculate surface including water :param selection: Calculate surface for subset of entity :param naccess_exe: naccess executable (full path to executable) + :param naccess_root: Path to folder containing "accall" binary and + files "vdw.radii" and "standard.data". This is the + fastest way to call naccess! :param keep_files: If True, do not delete temporary files :param asa_abs: Attaches per residue absolute SASA to specified FloatProp on residue level @@ -223,38 +275,70 @@ def CalculateSurfaceArea(entity, radius=1.4, """ # check if naccess executable is specified - naccess_executable=_GetExecutable(naccess_exe) + if naccess_root and _CheckNaccessRoot(naccess_root): + # use faster, direct call to accall binary + fast_mode = True + else: + # get naccess executable + naccess_executable = _GetExecutable(naccess_exe) + # see if we can extract naccess_root from there (fallback to old mode) + naccess_root = os.path.dirname(naccess_executable) + fast_mode = _CheckNaccessRoot(naccess_root) # setup files for naccess (naccess_data_dir, naccess_data_file, naccess_data_base) \ = _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms) - # 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 - # execute naccess - stdout_value = _RunNACCESS(command, naccess_data_dir) - - # parse outout - new_asa= os.path.join(naccess_data_dir, "%s.asa" % naccess_data_base) - _ParseAsaFile(entity, new_asa, asa_atom) + try: + # call naccess + if fast_mode: + # cook up stdin query (same logic as naccess script) + query = "PDBFILE %s\n" \ + "VDWFILE %s\n" \ + "STDFILE %s\n" \ + "PROBE %f\n" \ + "ZSLICE 0.05\n" \ + % (naccess_data_file, os.path.join(naccess_root, "vdw.radii"), + os.path.join(naccess_root, "standard.data"), radius) + if include_hydrogens: + query += "HYDROGENS\n" + if include_water: + query += "WATERS\n" + if include_hetatm: + query += "HETATOMS\n" + # call it + command = os.path.join(naccess_root, "accall") + _RunACCALL(command, naccess_data_dir, query) + else: + LogWarning("NACCESS: Falling back to slower call to %s." \ + % naccess_executable) + # 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 + # execute naccess + _RunNACCESS(command, naccess_data_dir) - new_rsa = os.path.join(naccess_data_dir, "%s.rsa" % naccess_data_base) - _ParseRsaFile(entity, new_rsa, asa_abs, asa_rel) + # parse outout + 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) + + finally: + # clean up + if not keep_files: + __CleanupFiles(naccess_data_dir) # sum up Asa for all atoms sasa = 0.0 for a in entity.atoms: sasa += a.GetFloatProp(asa_atom, 0.0) - - # clean up - if not keep_files: - __CleanupFiles(naccess_data_dir) return sasa