text: # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
  # This code is part of the Biopython distribution and governed by its
  # license.  Please see the LICENSE file that should have been included
  # as part of this package.
  
  # NACCESS interface adapted from Bio/PDB/DSSP.py
  
  import os, sys, tempfile
  from Bio.PDB.PDBIO import PDBIO
  from Bio.PDB.AbstractPropertyMap import AbstractResiduePropertyMap, AbstractAtomPropertyMap
  
  """Interface for the program NACCESS.
  
  See: http://wolf.bms.umist.ac.uk/naccess/
  
  errors likely to occur with the binary:
  default values are often due to low default settings in accall.pars
  - e.g. max cubes error: change in accall.pars and recompile binary
  
  use naccess -y, naccess -h or naccess -w to include HETATM records
  """
  
  def run_naccess(model, pdb_file, probe_size = None, z_slice = None, \
                  naccess = 'naccess', temp_path = '/tmp/'):
      
      # make temp directory; chdir to temp directory, 
      # as NACCESS writes to current working directory
      tmp_path = tempfile.mktemp(dir = temp_path)
      os.mkdir(tmp_path)
      old_dir = os.getcwd()
      os.chdir(tmp_path)
      
      # file name must end with '.pdb' to work with NACCESS
      # -> create temp file of existing pdb
      #    or write model to temp file
      tmp_pdb_file = tempfile.mktemp('.pdb', dir = tmp_path)
      if pdb_file:
          os.system('cp %s %s' % (pdb_file, tmp_pdb_file))
      else:
          writer = PDBIO()
          writer.set_structure(model.get_parent())
          writer.save(tmp_pdb_file)
  
      # create the command line and run
      # catch standard out & err
      command = '%s %s ' % (naccess, tmp_pdb_file)
      if probe_size:
          command += '-p %s ' % probe_size
      if z_slice:
          command += '-z %s ' % z_slice
      in_, out, err = os.popen3(command)
      in_.close()
      stdout = out.readlines()
      out.close()
      stderr = err.readlines()
      err.close()
  
      # get the output, then delete the temp directory
      rsa_file = tmp_pdb_file[:-4] + '.rsa'
      rf = open(rsa_file)
      rsa_data = rf.readlines()
      rf.close()
      asa_file = tmp_pdb_file[:-4] + '.asa'
      af = open(asa_file)
      asa_data = af.readlines()
      af.close()
      os.chdir(old_dir)
      os.system('rm -rf %s >& /dev/null' % tmp_path)
      return rsa_data, asa_data
  
  def process_rsa_data(rsa_data):
      # process the .rsa output file: residue level SASA data
      naccess_rel_dict = {}
      for line in rsa_data:
          if line.startswith('RES'):
              res_name = line[4:7]
              chain_id = line[8]
              resseq = int(line[9:13])
              icode = line[13]
              res_id = (' ', resseq, icode)
              naccess_rel_dict[(chain_id, res_id)] = { \
                  'res_name': res_name,
                  'all_atoms_abs': float(line[16:22]),
                  'all_atoms_rel': float(line[23:28]),
                  'side_chain_abs': float(line[29:35]),
                  'side_chain_rel': float(line[36:41]),
                  'main_chain_abs': float(line[42:48]),
                  'main_chain_rel': float(line[49:54]),
                  'non_polar_abs': float(line[55:61]),
                  'non_polar_rel': float(line[62:67]),
                  'all_polar_abs': float(line[68:74]),
                  'all_polar_rel': float(line[75:80]) } 
      return naccess_rel_dict
  
  def process_asa_data(rsa_data):
      # process the .asa output file: atomic level SASA data
      naccess_atom_dict = {}
      for line in rsa_data:
          atom_serial = line[6:11]
          full_atom_id = line[12:16]
          atom_id = full_atom_id.strip()
          altloc = line[16]
          resname = line[17:20]
          chainid = line[21]
          resseq = int(line[22:26])
          icode = line[26]
          res_id = (' ', resseq, icode)
          id = (chainid, res_id, atom_id)
          asa = line[54:62]               # solvent accessibility in Angstrom^2
          vdw = line[62:68]               # van der waal radius
          naccess_atom_dict[id] = asa
      return naccess_atom_dict
  
  
  class NACCESS(AbstractResiduePropertyMap):
      
      def __init__(self, model, pdb_file = None,
                   naccess_binary = 'naccess', tmp_directory = '/tmp'):
          res_data, atm_data = run_naccess(model, pdb_file, naccess = naccess_binary,
                                           temp_path = tmp_directory)
          naccess_dict = process_rsa_data(res_data)
          res_list = []
          property_dict={}
          property_keys=[]
          property_list=[]
          # Now create a dictionary that maps Residue objects to accessibility
          for chain in model:
              chain_id=chain.get_id()
              for res in chain:
                  res_id=res.get_id()
                  if (chain_id, res_id) in naccess_dict:
                      item = naccess_dict[(chain_id, res_id)]
                      res_name = item['res_name']
                      assert (res_name == res.get_resname())
                      property_dict[(chain_id, res_id)] = item
                      property_keys.append((chain_id, res_id))
                      property_list.append((res, item))
                      res.xtra["EXP_NACCESS"]=item
                  else:
                      pass
          AbstractResiduePropertyMap.__init__(self, property_dict, property_keys, 
                  property_list)
  
  class NACCESS_atomic(AbstractAtomPropertyMap):
  
      def __init__(self, model, pdb_file = None,
                   naccess_binary = 'naccess', tmp_directory = '/tmp'):
          res_data, atm_data = run_naccess(model, pdb_file, naccess = naccess_binary,
                                           temp_path = tmp_directory)
          self.naccess_atom_dict = process_asa_data(atm_data)
          atom_list = []
          property_dict={}
          property_keys=[]
          property_list=[]
          # Now create a dictionary that maps Atom objects to accessibility
          for chain in model:
              chain_id = chain.get_id()
              for residue in chain:
                  res_id = residue.get_id()
                  for atom in residue:
                      atom_id = atom.get_id()
                      full_id=(chain_id, res_id, atom_id)
                      if full_id in self.naccess_atom_dict:
                          asa = self.naccess_atom_dict[full_id]
                          property_dict[full_id]=asa
                          property_keys.append((full_id))
                          property_list.append((atom, asa))
                          atom.xtra['EXP_NACCESS']=asa
          AbstractAtomPropertyMap.__init__(self, property_dict, property_keys, 
                  property_list)
  
  
  if __name__=="__main__":
      
      import sys
      from Bio.PDB import PDBParser
      
      p=PDBParser()
      s=p.get_structure('X', sys.argv[1])
      model=s[0]
  
      n = NACCESS(model, sys.argv[1])
      for e in n.get_iterator():
          print e
  
  