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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
|
# 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.
"""mmCIF parsers"""
from __future__ import print_function
from string import ascii_letters
import numpy
import warnings
from Bio.File import as_handle
from Bio._py3k import range
from Bio.PDB.MMCIF2Dict import MMCIF2Dict
from Bio.PDB.StructureBuilder import StructureBuilder
from Bio.PDB.PDBExceptions import PDBConstructionException
from Bio.PDB.PDBExceptions import PDBConstructionWarning
class MMCIFParser(object):
"""Parse a mmCIF file and return a Structure object."""
def __init__(self, structure_builder=None, QUIET=False):
"""Create a PDBParser object.
The mmCIF parser calls a number of standard methods in an aggregated
StructureBuilder object. Normally this object is instanciated by the
MMCIParser object itself, but if the user provides his/her own
StructureBuilder object, the latter is used instead.
Arguments:
- structure_builder - an optional user implemented StructureBuilder class.
- QUIET - Evaluated as a Boolean. If true, warnings issued in constructing
the SMCRA data will be suppressed. If false (DEFAULT), they will be shown.
These warnings might be indicative of problems in the mmCIF file!
"""
if structure_builder is not None:
self._structure_builder = structure_builder
else:
self._structure_builder = StructureBuilder()
# self.header = None
# self.trailer = None
self.line_counter = 0
self.build_structure = None
self.QUIET = bool(QUIET)
# Public methods
def get_structure(self, structure_id, filename):
"""Return the structure.
Arguments:
- structure_id - string, the id that will be used for the structure
- filename - name of the mmCIF file OR an open filehandle
"""
with warnings.catch_warnings():
if self.QUIET:
warnings.filterwarnings("ignore", category=PDBConstructionWarning)
self._mmcif_dict = MMCIF2Dict(filename)
self._build_structure(structure_id)
return self._structure_builder.get_structure()
# Private methods
def _build_structure(self, structure_id):
mmcif_dict = self._mmcif_dict
atom_id_list = mmcif_dict["_atom_site.label_atom_id"]
residue_id_list = mmcif_dict["_atom_site.label_comp_id"]
try:
element_list = mmcif_dict["_atom_site.type_symbol"]
except KeyError:
element_list = None
seq_id_list = mmcif_dict["_atom_site.label_seq_id"]
chain_id_list = mmcif_dict["_atom_site.auth_asym_id"]
x_list = [float(x) for x in mmcif_dict["_atom_site.Cartn_x"]]
y_list = [float(x) for x in mmcif_dict["_atom_site.Cartn_y"]]
z_list = [float(x) for x in mmcif_dict["_atom_site.Cartn_z"]]
alt_list = mmcif_dict["_atom_site.label_alt_id"]
icode_list = mmcif_dict["_atom_site.pdbx_PDB_ins_code"]
b_factor_list = mmcif_dict["_atom_site.B_iso_or_equiv"]
occupancy_list = mmcif_dict["_atom_site.occupancy"]
fieldname_list = mmcif_dict["_atom_site.group_PDB"]
try:
serial_list = [int(n) for n in mmcif_dict["_atom_site.pdbx_PDB_model_num"]]
except KeyError:
# No model number column
serial_list = None
except ValueError:
# Invalid model number (malformed file)
raise PDBConstructionException("Invalid model number")
try:
aniso_u11 = mmcif_dict["_atom_site.aniso_U[1][1]"]
aniso_u12 = mmcif_dict["_atom_site.aniso_U[1][2]"]
aniso_u13 = mmcif_dict["_atom_site.aniso_U[1][3]"]
aniso_u22 = mmcif_dict["_atom_site.aniso_U[2][2]"]
aniso_u23 = mmcif_dict["_atom_site.aniso_U[2][3]"]
aniso_u33 = mmcif_dict["_atom_site.aniso_U[3][3]"]
aniso_flag = 1
except KeyError:
# no anisotropic B factors
aniso_flag = 0
# if auth_seq_id is present, we use this.
# Otherwise label_seq_id is used.
if "_atom_site.auth_seq_id" in mmcif_dict:
seq_id_list = mmcif_dict["_atom_site.auth_seq_id"]
else:
seq_id_list = mmcif_dict["_atom_site.label_seq_id"]
# Now loop over atoms and build the structure
current_chain_id = None
current_residue_id = None
structure_builder = self._structure_builder
structure_builder.init_structure(structure_id)
structure_builder.init_seg(" ")
# Historically, Biopython PDB parser uses model_id to mean array index
# so serial_id means the Model ID specified in the file
current_model_id = -1
current_serial_id = -1
for i in range(0, len(atom_id_list)):
# set the line_counter for 'ATOM' lines only and not
# as a global line counter found in the PDBParser()
# this number should match the '_atom_site.id' index in the MMCIF
structure_builder.set_line_counter(i)
x = x_list[i]
y = y_list[i]
z = z_list[i]
resname = residue_id_list[i]
chainid = chain_id_list[i]
altloc = alt_list[i]
if altloc == ".":
altloc = " "
int_resseq = int(seq_id_list[i])
icode = icode_list[i]
if icode == "?":
icode = " "
name = atom_id_list[i]
# occupancy & B factor
try:
tempfactor = float(b_factor_list[i])
except ValueError:
raise PDBConstructionException("Invalid or missing B factor")
try:
occupancy = float(occupancy_list[i])
except ValueError:
raise PDBConstructionException("Invalid or missing occupancy")
fieldname = fieldname_list[i]
if fieldname == "HETATM":
hetatm_flag = "H"
else:
hetatm_flag = " "
resseq = (hetatm_flag, int_resseq, icode)
if serial_list is not None:
# model column exists; use it
serial_id = serial_list[i]
if current_serial_id != serial_id:
# if serial changes, update it and start new model
current_serial_id = serial_id
current_model_id += 1
structure_builder.init_model(current_model_id, current_serial_id)
current_chain_id = None
current_residue_id = None
else:
# no explicit model column; initialize single model
structure_builder.init_model(current_model_id)
if current_chain_id != chainid:
current_chain_id = chainid
structure_builder.init_chain(current_chain_id)
current_residue_id = None
if current_residue_id != resseq:
current_residue_id = resseq
structure_builder.init_residue(resname, hetatm_flag, int_resseq, icode)
coord = numpy.array((x, y, z), 'f')
element = element_list[i] if element_list else None
structure_builder.init_atom(name, coord, tempfactor, occupancy, altloc,
name, element=element)
if aniso_flag == 1:
u = (aniso_u11[i], aniso_u12[i], aniso_u13[i],
aniso_u22[i], aniso_u23[i], aniso_u33[i])
mapped_anisou = [float(x) for x in u]
anisou_array = numpy.array(mapped_anisou, 'f')
structure_builder.set_anisou(anisou_array)
# Now try to set the cell
try:
a = float(mmcif_dict["_cell.length_a"])
b = float(mmcif_dict["_cell.length_b"])
c = float(mmcif_dict["_cell.length_c"])
alpha = float(mmcif_dict["_cell.angle_alpha"])
beta = float(mmcif_dict["_cell.angle_beta"])
gamma = float(mmcif_dict["_cell.angle_gamma"])
cell = numpy.array((a, b, c, alpha, beta, gamma), 'f')
spacegroup = mmcif_dict["_symmetry.space_group_name_H-M"]
spacegroup = spacegroup[1:-1] # get rid of quotes!!
if spacegroup is None:
raise Exception
structure_builder.set_symmetry(spacegroup, cell)
except Exception:
pass # no cell found, so just ignore
class FastMMCIFParser(object):
"""Parse an MMCIF file and return a Structure object."""
def __init__(self, structure_builder=None, QUIET=False):
"""Create a FastMMCIFParser object.
The mmCIF parser calls a number of standard methods in an aggregated
StructureBuilder object. Normally this object is instanciated by the
parser object itself, but if the user provides his/her own
StructureBuilder object, the latter is used instead.
The main difference between this class and the regular MMCIFParser is
that only 'ATOM' and 'HETATM' lines are parsed here. Use if you are
interested only in coordinate information.
Arguments:
- structure_builder - an optional user implemented StructureBuilder class.
- QUIET - Evaluated as a Boolean. If true, warnings issued in constructing
the SMCRA data will be suppressed. If false (DEFAULT), they will be shown.
These warnings might be indicative of problems in the mmCIF file!
"""
if structure_builder is not None:
self._structure_builder = structure_builder
else:
self._structure_builder = StructureBuilder()
self.line_counter = 0
self.build_structure = None
self.QUIET = bool(QUIET)
# Public methods
def get_structure(self, structure_id, filename):
"""Return the structure.
Arguments:
- structure_id - string, the id that will be used for the structure
- filename - name of the mmCIF file OR an open filehandle
"""
with warnings.catch_warnings():
if self.QUIET:
warnings.filterwarnings("ignore", category=PDBConstructionWarning)
with as_handle(filename) as handle:
self._build_structure(structure_id, handle)
return self._structure_builder.get_structure()
# Private methods
def _build_structure(self, structure_id, filehandle):
# Read only _atom_site. and atom_site_anisotrop entries
read_atom, read_aniso = False, False
_fields, _records = [], []
_anisof, _anisors = [], []
for line in filehandle:
if line.startswith('_atom_site.'):
read_atom = True
_fields.append(line.strip())
elif line.startswith('_atom_site_anisotrop.'):
read_aniso = True
_anisof.append(line.strip())
elif read_atom and line.startswith('#'):
read_atom = False
elif read_aniso and line.startswith('#'):
read_aniso = False
elif read_atom:
_records.append(line.strip())
elif read_aniso:
_anisors.append(line.strip())
# Dumping the shlex module here since this particular
# category should be rather straightforward.
# Quite a performance boost..
_record_tbl = zip(*map(str.split, _records))
_anisob_tbl = zip(*map(str.split, _anisors))
mmcif_dict = dict(zip(_fields, _record_tbl))
mmcif_dict.update(dict(zip(_anisof, _anisob_tbl)))
# Build structure object
atom_id_list = mmcif_dict["_atom_site.label_atom_id"]
residue_id_list = mmcif_dict["_atom_site.label_comp_id"]
try:
element_list = mmcif_dict["_atom_site.type_symbol"]
except KeyError:
element_list = None
seq_id_list = mmcif_dict["_atom_site.label_seq_id"]
chain_id_list = mmcif_dict["_atom_site.auth_asym_id"]
x_list = [float(x) for x in mmcif_dict["_atom_site.Cartn_x"]]
y_list = [float(x) for x in mmcif_dict["_atom_site.Cartn_y"]]
z_list = [float(x) for x in mmcif_dict["_atom_site.Cartn_z"]]
alt_list = mmcif_dict["_atom_site.label_alt_id"]
icode_list = mmcif_dict["_atom_site.pdbx_PDB_ins_code"]
b_factor_list = mmcif_dict["_atom_site.B_iso_or_equiv"]
occupancy_list = mmcif_dict["_atom_site.occupancy"]
fieldname_list = mmcif_dict["_atom_site.group_PDB"]
try:
serial_list = [int(n) for n in mmcif_dict["_atom_site.pdbx_PDB_model_num"]]
except KeyError:
# No model number column
serial_list = None
except ValueError:
# Invalid model number (malformed file)
raise PDBConstructionException("Invalid model number")
try:
aniso_u11 = mmcif_dict["_atom_site.aniso_U[1][1]"]
aniso_u12 = mmcif_dict["_atom_site.aniso_U[1][2]"]
aniso_u13 = mmcif_dict["_atom_site.aniso_U[1][3]"]
aniso_u22 = mmcif_dict["_atom_site.aniso_U[2][2]"]
aniso_u23 = mmcif_dict["_atom_site.aniso_U[2][3]"]
aniso_u33 = mmcif_dict["_atom_site.aniso_U[3][3]"]
aniso_flag = 1
except KeyError:
# no anisotropic B factors
aniso_flag = 0
# if auth_seq_id is present, we use this.
# Otherwise label_seq_id is used.
if "_atom_site.auth_seq_id" in mmcif_dict:
seq_id_list = mmcif_dict["_atom_site.auth_seq_id"]
else:
seq_id_list = mmcif_dict["_atom_site.label_seq_id"]
# Now loop over atoms and build the structure
current_chain_id = None
current_residue_id = None
structure_builder = self._structure_builder
structure_builder.init_structure(structure_id)
structure_builder.init_seg(" ")
# Historically, Biopython PDB parser uses model_id to mean array index
# so serial_id means the Model ID specified in the file
current_model_id = -1
current_serial_id = -1
for i in range(0, len(atom_id_list)):
# set the line_counter for 'ATOM' lines only and not
# as a global line counter found in the PDBParser()
# this number should match the '_atom_site.id' index in the MMCIF
structure_builder.set_line_counter(i)
x = x_list[i]
y = y_list[i]
z = z_list[i]
resname = residue_id_list[i]
chainid = chain_id_list[i]
altloc = alt_list[i]
if altloc == ".":
altloc = " "
int_resseq = int(seq_id_list[i])
icode = icode_list[i]
if icode == "?":
icode = " "
name = atom_id_list[i].strip('"') # Remove occasional " from quoted atom names (e.g. xNA)
# occupancy & B factor
try:
tempfactor = float(b_factor_list[i])
except ValueError:
raise PDBConstructionException("Invalid or missing B factor")
try:
occupancy = float(occupancy_list[i])
except ValueError:
raise PDBConstructionException("Invalid or missing occupancy")
fieldname = fieldname_list[i]
if fieldname == "HETATM":
hetatm_flag = "H"
else:
hetatm_flag = " "
resseq = (hetatm_flag, int_resseq, icode)
if serial_list is not None:
# model column exists; use it
serial_id = serial_list[i]
if current_serial_id != serial_id:
# if serial changes, update it and start new model
current_serial_id = serial_id
current_model_id += 1
structure_builder.init_model(current_model_id, current_serial_id)
current_chain_id = None
current_residue_id = None
else:
# no explicit model column; initialize single model
structure_builder.init_model(current_model_id)
if current_chain_id != chainid:
current_chain_id = chainid
structure_builder.init_chain(current_chain_id)
current_residue_id = None
if current_residue_id != resseq:
current_residue_id = resseq
structure_builder.init_residue(resname, hetatm_flag, int_resseq, icode)
coord = numpy.array((x, y, z), 'f')
element = element_list[i] if element_list else None
structure_builder.init_atom(name, coord, tempfactor, occupancy, altloc,
name, element=element)
if aniso_flag == 1:
u = (aniso_u11[i], aniso_u12[i], aniso_u13[i],
aniso_u22[i], aniso_u23[i], aniso_u33[i])
mapped_anisou = [float(x) for x in u]
anisou_array = numpy.array(mapped_anisou, 'f')
structure_builder.set_anisou(anisou_array)
if __name__ == "__main__":
import sys
if len(sys.argv) != 2:
print("Usage: python MMCIFparser.py filename")
raise SystemExit
filename = sys.argv[1]
p = MMCIFParser()
structure = p.get_structure("test", filename)
for model in structure.get_list():
print(model)
for chain in model.get_list():
print(chain)
print("Found %d residues." % len(chain.get_list()))
|
