1
2
3
4
5
6
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
25
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
32
33
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
43
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
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
70
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
94
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]
108 vdw = line[62:68]
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
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
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
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