/
NACCESS.py
198 lines (166 loc) · 6.95 KB
/
NACCESS.py
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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
# 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
"""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
"""
from __future__ import print_function
import os
import tempfile
import shutil
import subprocess
import warnings
from Bio.PDB.PDBIO import PDBIO
from Bio.PDB.AbstractPropertyMap import AbstractResiduePropertyMap, AbstractAtomPropertyMap
def run_naccess(model, pdb_file, probe_size=None, z_slice=None,
naccess='naccess', temp_path='/tmp/'):
# make temp directory;
tmp_path = tempfile.mkdtemp(dir=temp_path)
# file name must end with '.pdb' to work with NACCESS
# -> create temp file of existing pdb
# or write model to temp file
handle, tmp_pdb_file = tempfile.mkstemp('.pdb', dir=tmp_path)
os.close(handle)
if pdb_file:
pdb_file = os.path.abspath(pdb_file)
shutil.copy(pdb_file, tmp_pdb_file)
else:
writer = PDBIO()
writer.set_structure(model.get_parent())
writer.save(tmp_pdb_file)
# chdir to temp directory, as NACCESS writes to current working directory
old_dir = os.getcwd()
os.chdir(tmp_path)
# create the command line and run
# catch standard out & err
command = [naccess, tmp_pdb_file]
if probe_size:
command.extend(['-p', probe_size])
if z_slice:
command.extend(['-z', z_slice])
p = subprocess.Popen(command, universal_newlines=True,
stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()
os.chdir(old_dir)
rsa_file = tmp_pdb_file[:-4] + '.rsa'
asa_file = tmp_pdb_file[:-4] + '.asa'
# Alert user for errors
if err.strip():
warnings.warn(err)
if (not os.path.exists(rsa_file)) or (not os.path.exists(asa_file)):
raise Exception('NACCESS did not execute or finish properly.')
# get the output, then delete the temp directory
with open(rsa_file) as rf:
rsa_data = rf.readlines()
with open(asa_file) as af:
asa_data = af.readlines()
# shutil.rmtree(tmp_path, ignore_errors=True)
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:
full_atom_id = line[12:16]
atom_id = full_atom_id.strip()
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
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'):
"""Initialize the class."""
res_data, atm_data = run_naccess(model, pdb_file,
naccess=naccess_binary,
temp_path=tmp_directory)
naccess_dict = process_rsa_data(res_data)
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'):
"""Initialize the class."""
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)
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:
"""Initialize the class."""
print(e)