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