1
2
3
4
5
6
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
27
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
34
35
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
45
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
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
72
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
96
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]
110 vdw = line[62:68]
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
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
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
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