Package Bio :: Package PDB :: Module NACCESS
[hide private]
[frames] | no frames]

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