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  from __future__ import print_function 
  9   
 10  import os 
 11  import tempfile 
 12  import shutil 
 13  import subprocess 
 14  from Bio.PDB.PDBIO import PDBIO 
 15  from Bio.PDB.AbstractPropertyMap import AbstractResiduePropertyMap, AbstractAtomPropertyMap 
 16   
 17  """Interface for the program NACCESS. 
 18   
 19  See: http://wolf.bms.umist.ac.uk/naccess/ 
 20   
 21  errors likely to occur with the binary: 
 22  default values are often due to low default settings in accall.pars 
 23  - e.g. max cubes error: change in accall.pars and recompile binary 
 24   
 25  use naccess -y, naccess -h or naccess -w to include HETATM records 
 26  """ 
 27   
 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 # get the output, then delete the temp directory 66 rsa_file = tmp_pdb_file[:-4] + '.rsa' 67 with open(rsa_file) as rf: 68 rsa_data = rf.readlines() 69 asa_file = tmp_pdb_file[:-4] + '.asa' 70 with open(asa_file) as af: 71 asa_data = af.readlines() 72 73 shutil.rmtree(tmp_path, ignore_errors=True) 74 return rsa_data, asa_data
75 76
77 -def process_rsa_data(rsa_data):
78 # process the .rsa output file: residue level SASA data 79 naccess_rel_dict = {} 80 for line in rsa_data: 81 if line.startswith('RES'): 82 res_name = line[4:7] 83 chain_id = line[8] 84 resseq = int(line[9:13]) 85 icode = line[13] 86 res_id = (' ', resseq, icode) 87 naccess_rel_dict[(chain_id, res_id)] = { 88 'res_name': res_name, 89 'all_atoms_abs': float(line[16:22]), 90 'all_atoms_rel': float(line[23:28]), 91 'side_chain_abs': float(line[29:35]), 92 'side_chain_rel': float(line[36:41]), 93 'main_chain_abs': float(line[42:48]), 94 'main_chain_rel': float(line[49:54]), 95 'non_polar_abs': float(line[55:61]), 96 'non_polar_rel': float(line[62:67]), 97 'all_polar_abs': float(line[68:74]), 98 'all_polar_rel': float(line[75:80])} 99 return naccess_rel_dict
100 101
102 -def process_asa_data(rsa_data):
103 # process the .asa output file: atomic level SASA data 104 naccess_atom_dict = {} 105 for line in rsa_data: 106 atom_serial = line[6:11] 107 full_atom_id = line[12:16] 108 atom_id = full_atom_id.strip() 109 altloc = line[16] 110 resname = line[17:20] 111 chainid = line[21] 112 resseq = int(line[22:26]) 113 icode = line[26] 114 res_id = (' ', resseq, icode) 115 id = (chainid, res_id, atom_id) 116 asa = line[54:62] # solvent accessibility in Angstrom^2 117 vdw = line[62:68] # van der waal radius 118 naccess_atom_dict[id] = asa 119 return naccess_atom_dict
120 121
122 -class NACCESS(AbstractResiduePropertyMap):
123
124 - def __init__(self, model, pdb_file=None, 125 naccess_binary='naccess', tmp_directory='/tmp'):
126 res_data, atm_data = run_naccess(model, pdb_file, 127 naccess=naccess_binary, 128 temp_path=tmp_directory) 129 naccess_dict = process_rsa_data(res_data) 130 res_list = [] 131 property_dict = {} 132 property_keys = [] 133 property_list = [] 134 # Now create a dictionary that maps Residue objects to accessibility 135 for chain in model: 136 chain_id = chain.get_id() 137 for res in chain: 138 res_id = res.get_id() 139 if (chain_id, res_id) in naccess_dict: 140 item = naccess_dict[(chain_id, res_id)] 141 res_name = item['res_name'] 142 assert (res_name == res.get_resname()) 143 property_dict[(chain_id, res_id)] = item 144 property_keys.append((chain_id, res_id)) 145 property_list.append((res, item)) 146 res.xtra["EXP_NACCESS"] = item 147 else: 148 pass 149 AbstractResiduePropertyMap.__init__(self, property_dict, property_keys, 150 property_list)
151 152
153 -class NACCESS_atomic(AbstractAtomPropertyMap):
154
155 - def __init__(self, model, pdb_file=None, 156 naccess_binary='naccess', tmp_directory='/tmp'):
157 res_data, atm_data = run_naccess(model, pdb_file, 158 naccess=naccess_binary, 159 temp_path=tmp_directory) 160 self.naccess_atom_dict = process_asa_data(atm_data) 161 property_dict = {} 162 property_keys = [] 163 property_list = [] 164 # Now create a dictionary that maps Atom objects to accessibility 165 for chain in model: 166 chain_id = chain.get_id() 167 for residue in chain: 168 res_id = residue.get_id() 169 for atom in residue: 170 atom_id = atom.get_id() 171 full_id = (chain_id, res_id, atom_id) 172 if full_id in self.naccess_atom_dict: 173 asa = self.naccess_atom_dict[full_id] 174 property_dict[full_id] = asa 175 property_keys.append((full_id)) 176 property_list.append((atom, asa)) 177 atom.xtra['EXP_NACCESS'] = asa 178 AbstractAtomPropertyMap.__init__(self, property_dict, 179 property_keys, property_list)
180 181 182 if __name__ == "__main__": 183 import sys 184 from Bio.PDB import PDBParser 185 186 p = PDBParser() 187 s = p.get_structure('X', sys.argv[1]) 188 model = s[0] 189 190 n = NACCESS(model, sys.argv[1]) 191 for e in n: 192 print(e) 193