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