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