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