Source Code for Module Bio.PDB.NACCESS

  1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5   
  6  # NACCESS interface adapted from Bio/PDB/DSSP.py 
  7   
  8  import os 
  9  import tempfile 
 10  from Bio.PDB.PDBIO import PDBIO 
 11  from Bio.PDB.AbstractPropertyMap import AbstractResiduePropertyMap, AbstractAtomPropertyMap 
 12   
 13  """Interface for the program NACCESS. 
 14   
 15  See: http://wolf.bms.umist.ac.uk/naccess/ 
 16   
 17  errors likely to occur with the binary: 
 18  default values are often due to low default settings in accall.pars 
 19  - e.g. max cubes error: change in accall.pars and recompile binary 
 20   
 21  use naccess -y, naccess -h or naccess -w to include HETATM records 
 22  """ 
 23   
 24   
 25 -def run_naccess(model, pdb_file, probe_size = None, z_slice = None, 
 26                  naccess = 'naccess', temp_path = '/tmp/'): 
 27   
 28      # make temp directory; chdir to temp directory, 
 29      # as NACCESS writes to current working directory 
 30      tmp_path = tempfile.mktemp(dir = temp_path) 
 31      os.mkdir(tmp_path) 
 32      old_dir = os.getcwd() 
 33      os.chdir(tmp_path) 
 34   
 35      # file name must end with '.pdb' to work with NACCESS 
 36      # -> create temp file of existing pdb 
 37      #    or write model to temp file 
 38      tmp_pdb_file = tempfile.mktemp('.pdb', dir = tmp_path) 
 39      if pdb_file: 
 40          os.system('cp %s %s' % (pdb_file, tmp_pdb_file)) 
 41      else: 
 42          writer = PDBIO() 
 43          writer.set_structure(model.get_parent()) 
 44          writer.save(tmp_pdb_file) 
 45   
 46      # create the command line and run 
 47      # catch standard out & err 
 48      command = '%s %s ' % (naccess, tmp_pdb_file) 
 49      if probe_size: 
 50          command += '-p %s ' % probe_size 
 51      if z_slice: 
 52          command += '-z %s ' % z_slice 
 53      in_, out, err = os.popen3(command) 
 54      in_.close() 
 55      stdout = out.readlines() 
 56      out.close() 
 57      stderr = err.readlines() 
 58      err.close() 
 59   
 60      # get the output, then delete the temp directory 
 61      rsa_file = tmp_pdb_file[:-4] + '.rsa' 
 62      rf = open(rsa_file) 
 63      rsa_data = rf.readlines() 
 64      rf.close() 
 65      asa_file = tmp_pdb_file[:-4] + '.asa' 
 66      af = open(asa_file) 
 67      asa_data = af.readlines() 
 68      af.close() 
 69      os.chdir(old_dir) 
 70      os.system('rm -rf %s >& /dev/null' % tmp_path) 
 71      return rsa_data, asa_data 
 72   
 73   
 74 -def process_rsa_data(rsa_data): 
 75      # process the .rsa output file: residue level SASA data 
 76      naccess_rel_dict = {} 
 77      for line in rsa_data: 
 78          if line.startswith('RES'): 
 79              res_name = line[4:7] 
 80              chain_id = line[8] 
 81              resseq = int(line[9:13]) 
 82              icode = line[13] 
 83              res_id = (' ', resseq, icode) 
 84              naccess_rel_dict[(chain_id, res_id)] = { 
 85                  'res_name': res_name, 
 86                  'all_atoms_abs': float(line[16:22]), 
 87                  'all_atoms_rel': float(line[23:28]), 
 88                  'side_chain_abs': float(line[29:35]), 
 89                  'side_chain_rel': float(line[36:41]), 
 90                  'main_chain_abs': float(line[42:48]), 
 91                  'main_chain_rel': float(line[49:54]), 
 92                  'non_polar_abs': float(line[55:61]), 
 93                  'non_polar_rel': float(line[62:67]), 
 94                  'all_polar_abs': float(line[68:74]), 
 95                  'all_polar_rel': float(line[75:80]) } 
 96      return naccess_rel_dict 
 97   
 98   
 99 -def process_asa_data(rsa_data): 
100      # process the .asa output file: atomic level SASA data 
101      naccess_atom_dict = {} 
102      for line in rsa_data: 
103          atom_serial = line[6:11] 
104          full_atom_id = line[12:16] 
105          atom_id = full_atom_id.strip() 
106          altloc = line[16] 
107          resname = line[17:20] 
108          chainid = line[21] 
109          resseq = int(line[22:26]) 
110          icode = line[26] 
111          res_id = (' ', resseq, icode) 
112          id = (chainid, res_id, atom_id) 
113          asa = line[54:62]               # solvent accessibility in Angstrom^2 
114          vdw = line[62:68]               # van der waal radius 
115          naccess_atom_dict[id] = asa 
116      return naccess_atom_dict 
117   
118   
119 -class NACCESS(AbstractResiduePropertyMap): 
120   
121 -    def __init__(self, model, pdb_file = None, 
122                   naccess_binary = 'naccess', tmp_directory = '/tmp'): 
123          res_data, atm_data = run_naccess(model, pdb_file, naccess = naccess_binary, 
124                                           temp_path = tmp_directory) 
125          naccess_dict = process_rsa_data(res_data) 
126          res_list = [] 
127          property_dict={} 
128          property_keys=[] 
129          property_list=[] 
130          # Now create a dictionary that maps Residue objects to accessibility 
131          for chain in model: 
132              chain_id=chain.get_id() 
133              for res in chain: 
134                  res_id=res.get_id() 
135                  if (chain_id, res_id) in naccess_dict: 
136                      item = naccess_dict[(chain_id, res_id)] 
137                      res_name = item['res_name'] 
138                      assert (res_name == res.get_resname()) 
139                      property_dict[(chain_id, res_id)] = item 
140                      property_keys.append((chain_id, res_id)) 
141                      property_list.append((res, item)) 
142                      res.xtra["EXP_NACCESS"]=item 
143                  else: 
144                      pass 
145          AbstractResiduePropertyMap.__init__(self, property_dict, property_keys, 
146                  property_list) 
147   
148   
149 -class NACCESS_atomic(AbstractAtomPropertyMap): 
150   
151 -    def __init__(self, model, pdb_file = None, 
152                   naccess_binary = 'naccess', tmp_directory = '/tmp'): 
153          res_data, atm_data = run_naccess(model, pdb_file, naccess = naccess_binary, 
154                                           temp_path = tmp_directory) 
155          self.naccess_atom_dict = process_asa_data(atm_data) 
156          atom_list = [] 
157          property_dict={} 
158          property_keys=[] 
159          property_list=[] 
160          # Now create a dictionary that maps Atom objects to accessibility 
161          for chain in model: 
162              chain_id = chain.get_id() 
163              for residue in chain: 
164                  res_id = residue.get_id() 
165                  for atom in residue: 
166                      atom_id = atom.get_id() 
167                      full_id=(chain_id, res_id, atom_id) 
168                      if full_id in self.naccess_atom_dict: 
169                          asa = self.naccess_atom_dict[full_id] 
170                          property_dict[full_id]=asa 
171                          property_keys.append((full_id)) 
172                          property_list.append((atom, asa)) 
173                          atom.xtra['EXP_NACCESS']=asa 
174          AbstractAtomPropertyMap.__init__(self, property_dict, property_keys, 
175                  property_list) 
176   
177   
178  if __name__=="__main__": 
179      import sys 
180      from Bio.PDB import PDBParser 
181   
182      p=PDBParser() 
183      s=p.get_structure('X', sys.argv[1]) 
184      model=s[0] 
185   
186      n = NACCESS(model, sys.argv[1]) 
187      for e in n.get_iterator(): 
188          print e 
189