1 """@namespace IMP.pmi.mmcif
2 @brief Support for the mmCIF file format.
4 IMP has basic support for writing out files in mmCIF format, for
5 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
6 mmCIF files are currently generated by creating an
7 IMP.pmi.mmcif.ProtocolOutput class, and attaching it to an
8 IMP.pmi.representation.Representation object, after which any
9 generated models and metadata are collected and output as mmCIF.
12 from __future__
import print_function
41 import ihm.representation
43 import ihm.cross_linkers
45 def _assign_id(obj, seen_objs, obj_by_id):
46 """Assign a unique ID to obj, and track all ids in obj_by_id."""
47 if obj
not in seen_objs:
48 if not hasattr(obj,
'id'):
50 obj.id = len(obj_by_id)
51 seen_objs[obj] = obj.id
53 obj.id = seen_objs[obj]
56 def _get_by_residue(p):
57 """Determine whether the given particle represents a specific residue
58 or a more coarse-grained object."""
62 class _ComponentMapper(object):
63 """Map a Particle to a component name"""
64 def __init__(self, prot):
67 self.name =
'cif-output'
68 self.o.dictionary_pdbs[self.name] = self.prot
69 self.o._init_dictchain(self.name, self.prot, multichar_chain=
True)
71 def __getitem__(self, p):
72 protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
76 class _AsymMapper(object):
77 """Map a Particle to an asym_unit"""
78 def __init__(self, simo, prot):
80 self._cm = _ComponentMapper(prot)
81 self._seen_ranges = {}
83 def __getitem__(self, p):
84 protname = self._cm[p]
85 return self.simo.asym_units[protname]
87 def get_feature(self, ps):
88 """Get an ihm.restraint.Feature that covers the given particles"""
96 rng = asym(rind, rind)
100 rng = asym(rinds[0], rinds[-1])
102 raise ValueError(
"Unsupported particle type %s" % str(p))
104 if len(rngs) > 0
and rngs[-1].asym == asym \
105 and rngs[-1].seq_id_range[1] == rng.seq_id_range[0] - 1:
106 rngs[-1].seq_id_range = (rngs[-1].seq_id_range[0],
113 if hrngs
in self._seen_ranges:
114 return self._seen_ranges[hrngs]
116 feat = ihm.restraint.ResidueFeature(rngs)
117 self._seen_ranges[hrngs] = feat
121 class _AllSoftware(object):
122 def __init__(self, system):
124 self.modeller_used = self.phyre2_used =
False
125 self.pmi = ihm.Software(
126 name=
"IMP PMI module",
127 version=IMP.pmi.__version__,
128 classification=
"integrative model building",
129 description=
"integrative model building",
130 location=
'https://integrativemodeling.org')
131 self.imp = ihm.Software(
132 name=
"Integrative Modeling Platform (IMP)",
133 version=IMP.__version__,
134 classification=
"integrative model building",
135 description=
"integrative model building",
136 location=
'https://integrativemodeling.org')
137 self.system.software.extend([self.pmi, self.imp])
139 def set_modeller_used(self, version, date):
140 if self.modeller_used:
142 self.modeller_used =
True
143 self.system.software.append(ihm.Software(
144 name=
'MODELLER', classification=
'comparative modeling',
145 description=
'Comparative modeling by satisfaction '
146 'of spatial restraints, build ' + date,
147 location=
'https://salilab.org/modeller/',
150 def set_phyre2_used(self):
153 self.phyre2_used =
True
154 self.system.software.append(ihm.Software(
155 name=
'Phyre2', classification=
'protein homology modeling',
156 description=
'Protein Homology/analogY Recognition '
159 location=
'http://www.sbg.bio.ic.ac.uk/~phyre2/'))
161 def _get_fragment_is_rigid(fragment):
162 """Determine whether a fragment is modeled rigidly"""
167 class _PDBFragment(ihm.representation.ResidueSegment):
168 """Record details about part of a PDB file used as input
170 def __init__(self, state, component, start, end, pdb_offset,
171 pdbname, chain, hier, asym_unit):
173 super(_PDBFragment, self).__init__(
174 asym_unit=asym_unit.pmi_range(start, end),
175 rigid=
None, primitive=
'sphere')
176 self.component, self.start, self.end, self.offset, self.pdbname \
177 = component, start, end, pdb_offset, pdbname
178 self.state, self.chain, self.hier = state, chain, hier
183 rigid = property(
lambda self: _get_fragment_is_rigid(self),
184 lambda self, val:
None)
186 def combine(self, other):
190 class _BeadsFragment(ihm.representation.FeatureSegment):
191 """Record details about beads used to represent part of a component."""
193 def __init__(self, state, component, start, end, count, hier, asym_unit):
194 super(_BeadsFragment, self).__init__(asym_unit=asym_unit(start, end),
195 rigid=
None, primitive=
'sphere', count=count)
196 self.state, self.component, self.hier \
197 = state, component, hier
199 rigid = property(
lambda self: _get_fragment_is_rigid(self),
200 lambda self, val:
None)
202 def combine(self, other):
204 if type(other) == type(self)
and \
205 other.asym_unit.seq_id_range[0] == self.asym_unit.seq_id_range[1] + 1:
206 self.asym_unit.seq_id_range = (self.asym_unit.seq_id_range[0],
207 other.asym_unit.seq_id_range[1])
208 self.count += other.count
212 class _AllModelRepresentations(object):
213 def __init__(self, simo):
217 self.fragments = OrderedDict()
218 self._all_representations = {}
220 def copy_component(self, state, name, original, asym_unit):
221 """Copy all representation for `original` in `state` to `name`"""
224 newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
226 for rep
in self.fragments:
227 if original
in self.fragments[rep]:
228 if name
not in self.fragments[rep]:
229 self.fragments[rep][name] = OrderedDict()
230 self.fragments[rep][name][state] = [copy_frag(f)
231 for f
in self.fragments[rep][original][state]]
234 first_state = list(self.fragments[rep][name].keys())[0]
235 if state
is first_state:
236 representation = self._all_representations[rep]
237 representation.extend(self.fragments[rep][name][state])
239 def add_fragment(self, state, representation, fragment):
240 """Add a model fragment."""
241 comp = fragment.component
242 id_rep = id(representation)
243 self._all_representations[id_rep] = representation
244 if id_rep
not in self.fragments:
245 self.fragments[id_rep] = OrderedDict()
246 if comp
not in self.fragments[id_rep]:
247 self.fragments[id_rep][comp] = OrderedDict()
248 if state
not in self.fragments[id_rep][comp]:
249 self.fragments[id_rep][comp][state] = []
250 fragments = self.fragments[id_rep][comp][state]
251 if len(fragments) == 0
or not fragments[-1].combine(fragment):
252 fragments.append(fragment)
255 first_state = list(self.fragments[id_rep][comp].keys())[0]
256 if state
is first_state:
257 representation.append(fragment)
260 class _AllDatasets(object):
261 """Track all datasets generated by PMI and add them to the ihm.System"""
262 def __init__(self, system):
265 self._datasets_by_state = {}
266 self._restraints_by_state = {}
268 def get_all_group(self, state):
269 """Get a DatasetGroup encompassing all datasets so far in this state"""
273 g = ihm.dataset.DatasetGroup(self._datasets_by_state.get(state, [])
274 + [r.dataset
for r
in self._restraints_by_state.get(state, [])
278 def add(self, state, dataset):
279 """Add a new dataset."""
280 self._datasets.append(dataset)
281 if state
not in self._datasets_by_state:
282 self._datasets_by_state[state] = []
283 self._datasets_by_state[state].append(dataset)
286 self.system.orphan_datasets.append(dataset)
290 """Add the dataset for a restraint"""
291 if state
not in self._restraints_by_state:
292 self._restraints_by_state[state] = []
293 self._restraints_by_state[state].append(restraint)
296 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
297 """Restrain to a set of cross links"""
300 _label_map = {
'wtDSS':
'DSS',
'scDSS':
'DSS',
'scEDC':
'EDC'}
301 _descriptor_map = {
'DSS': ihm.cross_linkers.dss,
302 'EDC': ihm.cross_linkers.edc}
304 def __init__(self, pmi_restraint):
305 self.pmi_restraint = pmi_restraint
306 label = self.pmi_restraint.label
308 label = self._label_map.get(label, label)
310 super(_CrossLinkRestraint, self).__init__(
311 dataset=self.pmi_restraint.dataset,
312 linker=self._get_chem_descriptor(label))
315 def _get_chem_descriptor(cls, label):
316 if label
not in cls._descriptor_map:
320 d = ihm.ChemDescriptor(label)
321 cls._descriptor_map[label] = d
322 return cls._descriptor_map[label]
324 def _set_psi_sigma(self, model):
327 if model.m != self.pmi_restraint.model:
329 for resolution
in self.pmi_restraint.sigma_dictionary:
330 statname =
'ISDCrossLinkMS_Sigma_%s_%s' % (resolution, self.label)
331 if model.stats
and statname
in model.stats:
332 sigma = float(model.stats[statname])
333 p = self.pmi_restraint.sigma_dictionary[resolution][0]
334 old_values.append((p, p.get_scale()))
336 for psiindex
in self.pmi_restraint.psi_dictionary:
337 statname =
'ISDCrossLinkMS_Psi_%s_%s' % (psiindex, self.label)
338 if model.stats
and statname
in model.stats:
339 psi = float(model.stats[statname])
340 p = self.pmi_restraint.psi_dictionary[psiindex][0]
341 old_values.append((p, p.get_scale()))
345 return list(reversed(old_values))
347 def add_fits_from_model_statfile(self, model):
349 old_values = self._set_psi_sigma(model)
353 for xl
in self.cross_links:
355 xl.fits[model] = ihm.restraint.CrossLinkFit(
356 psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
358 for p, scale
in old_values:
362 def __set_dataset(self, val):
363 self.pmi_restraint.dataset = val
364 dataset = property(
lambda self: self.pmi_restraint.dataset,
368 def get_asym_mapper_for_state(simo, state, asym_map):
369 asym = asym_map.get(state,
None)
371 asym = _AsymMapper(simo, state.prot)
372 asym_map[state] = asym
375 class _PMICrossLink(object):
378 psi = property(
lambda self: self.psi_p.get_scale(),
379 lambda self, val:
None)
380 sigma1 = property(
lambda self: self.sigma1_p.get_scale(),
381 lambda self, val:
None)
382 sigma2 = property(
lambda self: self.sigma2_p.get_scale(),
383 lambda self, val:
None)
386 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
390 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
394 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
395 def __init__(self, state, pmi_restraint, image_number, resolution,
396 pixel_size, image_resolution, projection_number,
398 self.pmi_restraint, self.image_number = pmi_restraint, image_number
399 super(_EM2DRestraint, self).__init__(
400 dataset=pmi_restraint.datasets[image_number],
401 assembly=state.modeled_assembly,
402 segment=
False, number_raw_micrographs=micrographs_number,
403 pixel_size_width=pixel_size, pixel_size_height=pixel_size,
404 image_resolution=image_resolution,
405 number_of_projections=projection_number)
408 def __get_dataset(self):
409 return self.pmi_restraint.datasets[self.image_number]
410 def __set_dataset(self, val):
411 self.pmi_restraint.datasets[self.image_number] = val
412 dataset = property(__get_dataset, __set_dataset)
414 def add_fits_from_model_statfile(self, model):
415 ccc = self._get_cross_correlation(model)
416 transform = self._get_transformation(model)
417 rot = transform.get_rotation()
418 rm = [[e
for e
in rot.get_rotation_matrix_row(i)]
for i
in range(3)]
419 self.fits[model] = ihm.restraint.EM2DRestraintFit(
420 cross_correlation_coefficient=ccc,
422 tr_vector=transform.get_translation())
424 def _get_transformation(self, model):
425 """Get the transformation that places the model on the image"""
426 stats = model.em2d_stats
or model.stats
427 prefix =
'ElectronMicroscopy2D_%s_Image%d' % (self.pmi_restraint.label,
428 self.image_number + 1)
429 r = [float(stats[prefix +
'_Rotation%d' % i])
for i
in range(4)]
430 t = [float(stats[prefix +
'_Translation%d' % i])
434 inv = model.transform.get_inverse()
438 def _get_cross_correlation(self, model):
439 """Get the cross correlation coefficient between the model projection
441 stats = model.em2d_stats
or model.stats
442 return float(stats[
'ElectronMicroscopy2D_%s_Image%d_CCC'
443 % (self.pmi_restraint.label,
444 self.image_number + 1)])
447 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
449 def __init__(self, simo, state, pmi_restraint, target_ps, densities):
450 self.pmi_restraint = pmi_restraint
451 super(_EM3DRestraint, self).__init__(
452 dataset=pmi_restraint.dataset,
453 assembly=self._get_assembly(densities, simo, state),
454 fitting_method=
'Gaussian mixture models',
455 number_of_gaussians=len(target_ps))
458 def __set_dataset(self, val):
459 self.pmi_restraint.dataset = val
460 dataset = property(
lambda self: self.pmi_restraint.dataset,
463 def _get_assembly(self, densities, simo, state):
464 """Get the Assembly that this restraint acts on"""
465 cm = _ComponentMapper(state.prot)
468 components[cm[d]] =
None
469 a = simo._get_subassembly(components,
470 name=
"EM subassembly",
471 description=
"All components that fit the EM map")
474 def add_fits_from_model_statfile(self, model):
475 ccc = self._get_cross_correlation(model)
476 self.fits[model] = ihm.restraint.EM3DRestraintFit(
477 cross_correlation_coefficient=ccc)
479 def _get_cross_correlation(self, model):
480 """Get the cross correlation coefficient between the model
482 if model.stats
is not None:
483 return float(model.stats[
'GaussianEMRestraint_%s_CCC'
484 % self.pmi_restraint.label])
487 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
489 def __init__(self, simo, state, pmi_restraint, geometric_object,
490 feature, distance, sigma):
491 self.pmi_restraint = pmi_restraint
492 super(_GeometricRestraint, self).__init__(
493 dataset=pmi_restraint.dataset,
494 geometric_object=geometric_object, feature=feature,
495 distance=distance, harmonic_force_constant = 1. / sigma,
499 def __set_dataset(self, val):
500 self.pmi_restraint.dataset = val
501 dataset = property(
lambda self: self.pmi_restraint.dataset,
505 class _ReplicaExchangeProtocolStep(ihm.protocol.Step):
506 def __init__(self, state, rex):
507 if rex.monte_carlo_sample_objects
is not None:
508 method =
'Replica exchange monte carlo'
510 method =
'Replica exchange molecular dynamics'
511 self.monte_carlo_temperature = rex.vars[
'monte_carlo_temperature']
512 self.replica_exchange_minimum_temperature = \
513 rex.vars[
'replica_exchange_minimum_temperature']
514 self.replica_exchange_maximum_temperature = \
515 rex.vars[
'replica_exchange_maximum_temperature']
516 super(_ReplicaExchangeProtocolStep, self).__init__(
517 assembly=state.modeled_assembly,
519 method=method, name=
'Sampling',
520 num_models_begin=
None,
521 num_models_end=rex.vars[
"number_of_frames"],
522 multi_scale=
True, multi_state=
False, ordered=
False)
525 class _ReplicaExchangeProtocolDumper(ihm.dumper.Dumper):
526 """Write IMP-specific information about replica exchange to mmCIF.
527 Note that IDs will have already been assigned by python-ihm's
528 standard modeling protocol dumper."""
529 def dump(self, system, writer):
530 with writer.loop(
"_imp_replica_exchange_protocol",
531 [
"protocol_id",
"step_id",
"monte_carlo_temperature",
532 "replica_exchange_minimum_temperature",
533 "replica_exchange_maximum_temperature"])
as l:
534 for p
in system._all_protocols():
536 if isinstance(s, _ReplicaExchangeProtocolStep):
537 self._dump_step(p, s, l)
539 def _dump_step(self, p, s, l):
540 l.write(protocol_id=p._id, step_id=s._id,
541 monte_carlo_temperature=s.monte_carlo_temperature,
542 replica_exchange_minimum_temperature= \
543 s.replica_exchange_minimum_temperature,
544 replica_exchange_maximum_temperature= \
545 s.replica_exchange_maximum_temperature)
548 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
549 category =
'_imp_replica_exchange_protocol'
551 """Read IMP-specific information about replica exchange from mmCIF."""
552 def __call__(self, protocol_id, step_id, monte_carlo_temperature,
553 replica_exchange_minimum_temperature,
554 replica_exchange_maximum_temperature):
555 p = self.sysr.protocols.get_by_id(protocol_id)
557 s = p.steps[int(step_id)-1]
559 s.__class__ = _ReplicaExchangeProtocolStep
560 s.monte_carlo_temperature = \
561 self.get_float(monte_carlo_temperature)
562 s.replica_exchange_minimum_temperature = \
563 self.get_float(replica_exchange_minimum_temperature)
564 s.replica_exchange_maximum_temperature = \
565 self.get_float(replica_exchange_maximum_temperature)
568 class _SimpleProtocolStep(ihm.protocol.Step):
569 def __init__(self, state, num_models_end, method):
570 super(_SimpleProtocolStep, self).__init__(
571 assembly=state.modeled_assembly,
573 method=method, name=
'Sampling',
574 num_models_begin=
None,
575 num_models_end=num_models_end,
576 multi_scale=
True, multi_state=
False, ordered=
False)
579 class _Chain(object):
580 """Represent a single chain in a Model"""
581 def __init__(self, pmi_chain_id, asym_unit):
582 self.pmi_chain_id, self.asym_unit = pmi_chain_id, asym_unit
586 def add(self, xyz, atom_type, residue_type, residue_index,
587 all_indexes, radius):
588 if atom_type
is None:
589 self.spheres.append((xyz, residue_type, residue_index,
590 all_indexes, radius))
592 self.atoms.append((xyz, atom_type, residue_type, residue_index,
593 all_indexes, radius))
594 orig_comp = property(
lambda self: self.comp)
596 class _TransformedChain(object):
597 """Represent a chain that is a transformed version of another"""
598 def __init__(self, orig_chain, asym_unit, transform):
599 self.orig_chain, self.asym_unit = orig_chain, asym_unit
600 self.transform = transform
602 def __get_spheres(self):
603 for (xyz, residue_type, residue_index, all_indexes,
604 radius)
in self.orig_chain.spheres:
605 yield (self.transform * xyz, residue_type, residue_index,
607 spheres = property(__get_spheres)
609 def __get_atoms(self):
610 for (xyz, atom_type, residue_type, residue_index, all_indexes,
611 radius)
in self.orig_chain.atoms:
612 yield (self.transform * xyz, atom_type, residue_type, residue_index,
614 atoms = property(__get_atoms)
616 entity = property(
lambda self: self.orig_chain.entity)
617 orig_comp = property(
lambda self: self.orig_chain.comp)
619 class _Excluder(object):
620 def __init__(self, component, simo):
621 self._seqranges = simo._exclude_coords.get(component, [])
623 def is_excluded(self, indexes):
624 """Return True iff the given sequence range is excluded."""
625 for seqrange
in self._seqranges:
626 if indexes[0] >= seqrange[0]
and indexes[-1] <= seqrange[1]:
630 class _Model(ihm.model.Model):
631 def __init__(self, prot, simo, protocol, assembly, representation):
632 super(_Model, self).__init__(assembly=assembly, protocol=protocol,
633 representation=representation)
634 self.simo = weakref.proxy(simo)
640 self.em2d_stats =
None
643 self._is_restrained =
True
646 self.m = prot.get_model()
647 o.dictionary_pdbs[name] = prot
648 o._init_dictchain(name, prot, multichar_chain=
True)
649 (particle_infos_for_pdb,
650 self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
652 self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
655 def all_chains(self, simo):
656 """Yield all chains, including transformed ones"""
658 for c
in self.chains:
660 chain_for_comp[c.comp] = c
661 for tc
in simo._transformed_components:
662 orig_chain = chain_for_comp.get(tc.original,
None)
664 asym = simo.asym_units[tc.name]
665 c = _TransformedChain(orig_chain, asym, tc.transform)
669 def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
670 entity_for_chain = {}
673 for protname, chain_id
in o.dictchain[name].items():
674 if protname
in simo.entities:
675 entity_for_chain[chain_id] = simo.entities[protname]
678 pn = protname.split(
'.')[0]
679 entity_for_chain[chain_id] = simo.entities[pn]
680 comp_for_chain[chain_id] = protname
684 correct_asym[chain_id] = simo.asym_units[protname]
691 for (xyz, atom_type, residue_type, chain_id, residue_index,
692 all_indexes, radius)
in particle_infos_for_pdb:
693 if chain
is None or chain.pmi_chain_id != chain_id:
694 chain = _Chain(chain_id, correct_asym[chain_id])
695 chain.entity = entity_for_chain[chain_id]
696 chain.comp = comp_for_chain[chain_id]
697 self.chains.append(chain)
698 excluder = _Excluder(chain.comp, simo)
699 if not excluder.is_excluded(all_indexes
if all_indexes
700 else [residue_index]):
701 chain.add(xyz, atom_type, residue_type, residue_index,
704 def parse_rmsf_file(self, fname, component):
705 self.rmsf[component] = rmsf = {}
706 with open(fname)
as fh:
708 resnum, blocknum, val = line.split()
709 rmsf[int(resnum)] = (int(blocknum), float(val))
711 def get_rmsf(self, component, indexes):
712 """Get the RMSF value for the given residue indexes."""
715 rmsf = self.rmsf[component]
716 blocknums = dict.fromkeys(rmsf[ind][0]
for ind
in indexes)
717 if len(blocknums) != 1:
718 raise ValueError(
"Residue indexes %s aren't all in the same block"
720 return rmsf[indexes[0]][1]
723 for chain
in self.all_chains(self.simo):
724 pmi_offset = chain.asym_unit.entity.pmi_offset
725 for atom
in chain.atoms:
726 (xyz, atom_type, residue_type, residue_index,
727 all_indexes, radius) = atom
728 pt = self.transform * xyz
729 yield ihm.model.Atom(asym_unit=chain.asym_unit,
730 seq_id=residue_index - pmi_offset,
731 atom_id=atom_type.get_string(),
733 x=pt[0], y=pt[1], z=pt[2])
735 def get_spheres(self):
736 for chain
in self.all_chains(self.simo):
737 pmi_offset = chain.asym_unit.entity.pmi_offset
738 for sphere
in chain.spheres:
739 (xyz, residue_type, residue_index,
740 all_indexes, radius) = sphere
741 if all_indexes
is None:
742 all_indexes = (residue_index,)
743 pt = self.transform * xyz
744 yield ihm.model.Sphere(asym_unit=chain.asym_unit,
745 seq_id_range=(all_indexes[0] - pmi_offset,
746 all_indexes[-1] - pmi_offset),
747 x=pt[0], y=pt[1], z=pt[2], radius=radius,
748 rmsf=self.get_rmsf(chain.orig_comp, all_indexes))
751 class _AllProtocols(object):
752 def __init__(self, simo):
753 self.simo = weakref.proxy(simo)
755 self.protocols = OrderedDict()
757 def add_protocol(self, state):
758 """Add a new Protocol"""
759 if state
not in self.protocols:
760 self.protocols[state] = []
761 p = ihm.protocol.Protocol()
762 self.simo.system.orphan_protocols.append(p)
763 self.protocols[state].append(p)
765 def add_step(self, step, state):
766 """Add a ProtocolStep to the last Protocol of the given State"""
767 if state
not in self.protocols:
768 self.add_protocol(state)
769 protocol = self.get_last_protocol(state)
770 if len(protocol.steps) == 0:
771 step.num_models_begin = 0
773 step.num_models_begin = protocol.steps[-1].num_models_end
774 protocol.steps.append(step)
775 step.id = len(protocol.steps)
777 step.dataset_group = self.simo.all_datasets.get_all_group(state)
779 def add_postproc(self, step, state):
780 """Add a postprocessing step to the last protocol"""
781 protocol = self.get_last_protocol(state)
782 if len(protocol.analyses) == 0:
783 protocol.analyses.append(ihm.analysis.Analysis())
784 protocol.analyses[-1].steps.append(step)
786 def get_last_protocol(self, state):
787 """Return the most recently-added _Protocol"""
788 return self.protocols[state][-1]
791 class _AllStartingModels(object):
792 def __init__(self, simo):
796 self.models = OrderedDict()
799 def add_pdb_fragment(self, fragment):
800 """Add a starting model PDB fragment."""
801 comp = fragment.component
802 state = fragment.state
803 if comp
not in self.models:
804 self.models[comp] = OrderedDict()
805 if state
not in self.models[comp]:
806 self.models[comp][state] = []
807 models = self.models[comp][state]
808 if len(models) == 0 \
809 or models[-1].fragments[0].pdbname != fragment.pdbname:
810 model = self._add_model(fragment)
814 models[-1].fragments.append(weakref.proxy(fragment))
816 pmi_offset = models[-1].asym_unit.entity.pmi_offset
817 sid_begin = min(fragment.start + fragment.offset - pmi_offset,
818 models[-1].asym_unit.seq_id_range[0])
819 sid_end = max(fragment.end + fragment.offset - pmi_offset,
820 models[-1].asym_unit.seq_id_range[1])
821 models[-1].asym_unit = fragment.asym_unit.asym(sid_begin, sid_end)
822 fragment.starting_model = models[-1]
824 def _add_model(self, f):
825 parser = ihm.metadata.PDBParser()
826 r = parser.parse_file(f.pdbname)
828 self.simo._add_dataset(r[
'dataset'])
830 templates = r[
'templates'].get(f.chain, [])
833 self.simo.system.locations.append(t.alignment_file)
835 self.simo._add_dataset(t.dataset)
836 source = r[
'entity_source'].get(f.chain)
838 f.asym_unit.entity.source = source
839 pmi_offset = f.asym_unit.entity.pmi_offset
841 asym_unit=f.asym_unit.asym.pmi_range(f.start, f.end),
842 dataset=r[
'dataset'], asym_id=f.chain,
843 templates=templates, offset=f.offset + pmi_offset,
844 metadata=r[
'metadata'],
845 software=r[
'software'][0]
if r[
'software']
else None,
846 script_file=r[
'script'])
847 m.fragments = [weakref.proxy(f)]
851 class _StartingModel(ihm.startmodel.StartingModel):
852 def get_seq_dif(self):
856 pmi_offset = self.asym_unit.entity.pmi_offset
857 mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
859 for f
in self.fragments:
861 residue_indexes=list(range(f.start - f.offset,
862 f.end - f.offset + 1)))
863 for a
in mh.get_ihm_atoms(sel.get_selected_particles(),
864 f.offset - pmi_offset):
866 self._seq_dif = mh._seq_dif
869 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
870 """Post processing using AnalysisReplicaExchange0 macro"""
872 def __init__(self, rex, num_models_begin):
875 for fname
in self.get_all_stat_files():
876 with open(fname)
as fh:
877 num_models_end += len(fh.readlines())
878 super(_ReplicaExchangeAnalysisPostProcess, self).__init__(
879 feature=
'RMSD', num_models_begin=num_models_begin,
880 num_models_end=num_models_end)
882 def get_stat_file(self, cluster_num):
883 return os.path.join(self.rex._outputdir,
"cluster.%d" % cluster_num,
886 def get_all_stat_files(self):
887 for i
in range(self.rex._number_of_clusters):
888 yield self.get_stat_file(i)
891 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
892 """Ensemble generated using AnalysisReplicaExchange0 macro"""
894 num_models_deposited =
None
896 def __init__(self, pp, cluster_num, model_group, num_deposit):
897 with open(pp.get_stat_file(cluster_num))
as fh:
898 num_models = len(fh.readlines())
899 super(_ReplicaExchangeAnalysisEnsemble, self).__init__(
900 num_models=num_models,
901 model_group=model_group, post_process=pp,
902 clustering_feature=pp.feature,
903 name=model_group.name)
904 self.cluster_num = cluster_num
905 self.num_models_deposited = num_deposit
907 def get_rmsf_file(self, component):
908 return os.path.join(self.post_process.rex._outputdir,
909 'cluster.%d' % self.cluster_num,
910 'rmsf.%s.dat' % component)
912 def load_rmsf(self, model, component):
913 fname = self.get_rmsf_file(component)
914 if os.path.exists(fname):
915 model.parse_rmsf_file(fname, component)
917 def get_localization_density_file(self, fname):
918 return os.path.join(self.post_process.rex._outputdir,
919 'cluster.%d' % self.cluster_num,
922 def load_localization_density(self, state, fname, select_tuple, asym_units):
923 fullpath = self.get_localization_density_file(fname)
924 if os.path.exists(fullpath):
925 details =
"Localization density for %s %s" \
926 % (fname, self.model_group.name)
927 local_file = ihm.location.OutputFileLocation(fullpath,
929 for s
in select_tuple:
930 if isinstance(s, tuple)
and len(s) == 3:
931 asym = asym_units[s[2]].pmi_range(s[0], s[1])
934 den = ihm.model.LocalizationDensity(file=local_file,
936 self.densities.append(den)
938 def load_all_models(self, simo, state):
939 stat_fname = self.post_process.get_stat_file(self.cluster_num)
941 with open(stat_fname)
as fh:
942 stats = ast.literal_eval(fh.readline())
944 rmf_file = os.path.join(os.path.dirname(stat_fname),
945 "%d.rmf3" % model_num)
946 for c
in state.all_modeled_components:
948 state._pmi_object.set_coordinates_from_rmf(c, rmf_file, 0,
949 force_rigid_update=
True)
953 if model_num >= self.num_models_deposited:
957 def _get_precision(self):
958 precfile = os.path.join(self.post_process.rex._outputdir,
959 "precision.%d.%d.out" % (self.cluster_num,
961 if not os.path.exists(precfile):
964 r = re.compile(
'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
966 with open(precfile)
as fh:
970 return float(m.group(1))
972 precision = property(
lambda self: self._get_precision(),
973 lambda self, val:
None)
975 class _SimpleEnsemble(ihm.model.Ensemble):
976 """Simple manually-created ensemble"""
978 num_models_deposited =
None
980 def __init__(self, pp, model_group, num_models, drmsd,
981 num_models_deposited, ensemble_file):
982 super(_SimpleEnsemble, self).__init__(
983 model_group=model_group, post_process=pp, num_models=num_models,
984 file=ensemble_file, precision=drmsd, name=model_group.name,
985 clustering_feature=
'dRMSD')
986 self.num_models_deposited = num_models_deposited
988 def load_localization_density(self, state, component, asym, local_file):
989 den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
990 self.densities.append(den)
993 class _EntityMapper(dict):
994 """Handle mapping from IMP components (without copy number) to CIF entities.
995 Multiple components may map to the same entity if they share sequence."""
996 def __init__(self, system):
997 super(_EntityMapper, self).__init__()
998 self._sequence_dict = {}
1000 self.system = system
1002 def add(self, component_name, sequence, offset):
1003 def entity_seq(sequence):
1006 return [
'UNK' if s ==
'X' else s
for s
in sequence]
1009 if sequence
not in self._sequence_dict:
1012 d = component_name.split(
"@")[0].split(
".")[0]
1013 entity = Entity(entity_seq(sequence), description=d,
1015 self.system.entities.append(entity)
1016 self._sequence_dict[sequence] = entity
1017 self[component_name] = self._sequence_dict[sequence]
1020 class _TransformedComponent(object):
1021 def __init__(self, name, original, transform):
1022 self.name, self.original, self.transform = name, original, transform
1025 class _SimpleRef(object):
1026 """Class with similar interface to weakref.ref, but keeps a strong ref"""
1027 def __init__(self, ref):
1033 class _State(ihm.model.State):
1034 """Representation of a single state in the system."""
1036 def __init__(self, pmi_object, po):
1041 self._pmi_object = weakref.proxy(pmi_object)
1042 if hasattr(pmi_object,
'state'):
1045 self._pmi_state = _SimpleRef(pmi_object.state)
1047 self._pmi_state = weakref.ref(pmi_object)
1049 old_name = self.name
1050 super(_State, self).__init__(experiment_type=
'Fraction of bulk')
1051 self.name = old_name
1055 self.modeled_assembly = ihm.Assembly(
1056 name=
"Modeled assembly",
1057 description=self.get_postfixed_name(
1058 "All components modeled by IMP"))
1059 po.system.orphan_assemblies.append(self.modeled_assembly)
1061 self.all_modeled_components = []
1064 return hash(self._pmi_state())
1065 def __eq__(self, other):
1066 return self._pmi_state() == other._pmi_state()
1068 def add_model_group(self, group):
1072 def get_prefixed_name(self, name):
1073 """Prefix the given name with the state name, if available."""
1075 return self.short_name +
' ' + name
1079 return name[0].upper() + name[1:]
if name
else ''
1081 def get_postfixed_name(self, name):
1082 """Postfix the given name with the state name, if available."""
1084 return "%s in state %s" % (name, self.short_name)
1088 short_name = property(
lambda self: self._pmi_state().short_name)
1089 long_name = property(
lambda self: self._pmi_state().long_name)
1090 def __get_name(self):
1091 return self._pmi_state().long_name
1092 def __set_name(self, val):
1093 self._pmi_state().long_name = val
1094 name = property(__get_name, __set_name)
1098 """A single entity in the system.
1099 This functions identically to the base ihm.Entity class, but it
1100 allows identifying residues by either the IHM numbering scheme
1101 (seq_id, which is always contiguous starting at 1) or by the PMI
1102 scheme (which is similar, but may not start at 1 if the sequence in
1103 the FASTA file has one or more N-terminal gaps"""
1104 def __init__(self, sequence, pmi_offset, *args, **keys):
1107 self.pmi_offset = pmi_offset
1108 super(Entity, self).__init__(sequence, *args, **keys)
1111 """Return a single IHM residue indexed using PMI numbering"""
1112 return self.residue(res_id - self.pmi_offset)
1115 """Return a range of IHM residues indexed using PMI numbering"""
1116 off = self.pmi_offset
1117 return self(res_id_begin - off, res_id_end - off)
1121 """A single asymmetric unit in the system.
1122 This functions identically to the base ihm.AsymUnit class, but it
1123 allows identifying residues by either the IHM numbering scheme
1124 (seq_id, which is always contiguous starting at 1) or by the PMI
1125 scheme (which is similar, but may not start at 1 if the sequence in
1126 the FASTA file has one or more N-terminal gaps"""
1128 def __init__(self, entity, *args, **keys):
1129 super(AsymUnit, self).__init__(
1130 entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1133 """Return a single IHM residue indexed using PMI numbering"""
1134 return self.residue(res_id - self.entity.pmi_offset)
1137 """Return a range of IHM residues indexed using PMI numbering"""
1138 off = self.entity.pmi_offset
1139 return self(res_id_begin - off, res_id_end - off)
1143 """Class to encode a modeling protocol as mmCIF.
1145 IMP has basic support for writing out files in mmCIF format, for
1146 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1147 After creating an instance of this class, attach it to an
1148 IMP.pmi.representation.Representation object. After this, any
1149 generated models and metadata are automatically collected and
1152 def __init__(self, fh):
1154 self.system = ihm.System()
1155 self._state_group = ihm.model.StateGroup()
1156 self.system.state_groups.append(self._state_group)
1159 self._state_ensemble_offset = 0
1160 self._each_metadata = []
1161 self._file_datasets = []
1162 self._main_script = os.path.abspath(sys.argv[0])
1165 loc = ihm.location.WorkflowFileLocation(path=self._main_script,
1166 details=
"The main integrative modeling script")
1167 self.system.locations.append(loc)
1170 self.__asym_states = {}
1171 self._working_directory = os.getcwd()
1173 "Default representation")
1174 self.entities = _EntityMapper(self.system)
1176 self.asym_units = {}
1177 self._all_components = {}
1178 self.all_modeled_components = []
1179 self._transformed_components = []
1180 self.sequence_dict = {}
1183 self._xy_plane = ihm.geometry.XYPlane()
1184 self._xz_plane = ihm.geometry.XZPlane()
1185 self._z_axis = ihm.geometry.ZAxis()
1186 self._center_origin = ihm.geometry.Center(0,0,0)
1187 self._identity_transform = ihm.geometry.Transformation.identity()
1190 self._exclude_coords = {}
1192 self.all_representations = _AllModelRepresentations(self)
1193 self.all_protocols = _AllProtocols(self)
1194 self.all_datasets = _AllDatasets(self.system)
1195 self.all_starting_models = _AllStartingModels(self)
1197 self.all_software = _AllSoftware(self.system)
1200 """Create a new Representation and return it. This can be
1201 passed to add_model(), add_bead_element() or add_pdb_element()."""
1202 r = ihm.representation.Representation()
1203 self.system.orphan_representations.append(r)
1207 """Don't record coordinates for the given domain.
1208 Coordinates for the given domain (specified by a component name
1209 and a 2-element tuple giving the start and end residue numbers)
1210 will be excluded from the mmCIF file. This can be used to exclude
1211 parts of the structure that weren't well resolved in modeling.
1212 Any bead or residue that lies wholly within this range will be
1213 excluded. Multiple ranges for a given component can be excluded
1214 by calling this method multiple times."""
1215 if component
not in self._exclude_coords:
1216 self._exclude_coords[component] = []
1217 self._exclude_coords[component].append(seqrange)
1219 def _is_excluded(self, component, start, end):
1220 """Return True iff this chunk of sequence should be excluded"""
1221 for seqrange
in self._exclude_coords.get(component, ()):
1222 if start >= seqrange[0]
and end <= seqrange[1]:
1225 def _add_state(self, state):
1226 """Create a new state and return a pointer to it."""
1227 self._state_ensemble_offset = len(self.system.ensembles)
1228 s = _State(state, self)
1229 self._state_group.append(s)
1230 self._last_state = s
1233 def _get_chain_for_component(self, name, output):
1234 """Get the chain ID for a component, if any."""
1236 if name
in self.asym_units:
1237 return self.asym_units[name]._id
1242 def _get_assembly_comps(self, assembly):
1243 """Get the names of the components in the given assembly"""
1247 comps[ca.details] =
None
1251 """Make a new component that's a transformed copy of another.
1252 All representation for the existing component is copied to the
1254 assembly_comps = self._get_assembly_comps(state.modeled_assembly)
1255 if name
in assembly_comps:
1256 raise ValueError(
"Component %s already exists" % name)
1257 elif original
not in assembly_comps:
1258 raise ValueError(
"Original component %s does not exist" % original)
1259 self.create_component(state, name,
True)
1260 self.add_component_sequence(state, name, self.sequence_dict[original])
1261 self._transformed_components.append(_TransformedComponent(
1262 name, original, transform))
1263 self.all_representations.copy_component(state, name, original,
1264 self.asym_units[name])
1266 def create_component(self, state, name, modeled, asym_name=None):
1267 if asym_name
is None:
1269 new_comp = name
not in self._all_components
1270 self._all_components[name] =
None
1272 state.all_modeled_components.append(name)
1275 self.asym_units[asym_name] =
None
1276 self.all_modeled_components.append(name)
1278 def add_component_sequence(self, state, name, seq, asym_name=None):
1279 if asym_name
is None:
1281 def get_offset(seq):
1283 for i
in range(len(seq)):
1286 seq, offset = get_offset(seq)
1287 if name
in self.sequence_dict:
1288 if self.sequence_dict[name] != seq:
1289 raise ValueError(
"Sequence mismatch for component %s" % name)
1291 self.sequence_dict[name] = seq
1292 self.entities.add(name, seq, offset)
1293 if asym_name
in self.asym_units:
1294 if self.asym_units[asym_name]
is None:
1296 entity = self.entities[name]
1297 asym =
AsymUnit(entity, details=asym_name)
1298 self.system.asym_units.append(asym)
1299 self.asym_units[asym_name] = asym
1300 state.modeled_assembly.append(self.asym_units[asym_name])
1302 def _add_restraint_model_fits(self):
1303 """Add fits to restraints for all known models"""
1304 for group, m
in self.system._all_models():
1305 if m._is_restrained:
1306 for r
in self.system.restraints:
1307 if hasattr(r,
'add_fits_from_model_statfile'):
1308 r.add_fits_from_model_statfile(m)
1311 """Write out all information to the file.
1312 Information can be written in any format supported by
1313 the ihm library (typically this is 'mmCIF' or 'BCIF')."""
1315 ihm.dumper.write(self.fh, [self.system], format,
1316 dumpers=[_ReplicaExchangeProtocolDumper])
1319 """Do any final processing on the class hierarchy.
1320 After calling this method, the `system` member (an instance
1321 of `ihm.System`) completely reproduces the PMI modeling, and
1322 can be written out to an mmCIF file with `ihm.dumper.write`,
1323 and/or modified using the ihm API."""
1324 self._add_restraint_model_fits()
1327 self.system.software.extend(m
for m
in self._metadata
1328 if isinstance(m, ihm.Software))
1329 self.system.citations.extend(m
for m
in self._metadata
1330 if isinstance(m, ihm.Citation))
1331 self.system.locations.extend(m
for m
in self._metadata
1332 if isinstance(m, ihm.location.FileLocation))
1335 all_repos = [m
for m
in self._metadata
1336 if isinstance(m, ihm.location.Repository)]
1337 self.system.update_locations_in_repositories(all_repos)
1339 def add_pdb_element(self, state, name, start, end, offset, pdbname,
1340 chain, hier, representation=
None):
1341 if self._is_excluded(name, start, end):
1343 if representation
is None:
1344 representation = self.default_representation
1345 asym = self.asym_units[name]
1346 p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1348 self.all_representations.add_fragment(state, representation, p)
1349 self.all_starting_models.add_pdb_fragment(p)
1351 def add_bead_element(self, state, name, start, end, num, hier,
1352 representation=
None):
1353 if self._is_excluded(name, start, end):
1355 if representation
is None:
1356 representation = self.default_representation
1357 asym = self.asym_units[name]
1358 pmi_offset = asym.entity.pmi_offset
1359 b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1361 self.all_representations.add_fragment(state, representation, b)
1363 def get_cross_link_group(self, pmi_restraint):
1364 r = _CrossLinkRestraint(pmi_restraint)
1365 self.system.restraints.append(r)
1366 self._add_restraint_dataset(r)
1369 def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1370 if c1
not in self._all_components
or c2
not in self._all_components:
1376 e1 = self.entities[c1]
1377 e2 = self.entities[c2]
1378 xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1379 residue2=e2.pmi_residue(r2))
1380 rsr.experimental_cross_links.append([xl])
1383 def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1386 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1387 d = ihm.restraint.UpperBoundDistanceRestraint(length)
1389 if _get_by_residue(p1)
and _get_by_residue(p2):
1390 cls = _ResidueCrossLink
1392 cls = _FeatureCrossLink
1393 xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1396 xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1397 rsr.cross_links.append(xl)
1399 def add_replica_exchange(self, state, rex):
1404 step = _ReplicaExchangeProtocolStep(state, rex)
1405 step.software = self.all_software.pmi
1406 self.all_protocols.add_step(step, state)
1408 def _add_simple_dynamics(self, num_models_end, method):
1410 state = self._last_state
1411 self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1414 def _add_protocol(self):
1416 state = self._last_state
1417 self.all_protocols.add_protocol(state)
1419 def _add_dataset(self, dataset):
1420 return self.all_datasets.add(self._last_state, dataset)
1422 def _add_restraint_dataset(self, restraint):
1423 return self.all_datasets.add_restraint(self._last_state, restraint)
1425 def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1427 state = self._last_state
1428 pp = ihm.analysis.ClusterStep(
'RMSD', num_models_begin, num_models_end)
1429 self.all_protocols.add_postproc(pp, state)
1432 def _add_no_postprocessing(self, num_models):
1434 state = self._last_state
1435 pp = ihm.analysis.EmptyStep()
1436 pp.num_models_begin = pp.num_models_end = num_models
1437 self.all_protocols.add_postproc(pp, state)
1440 def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1441 num_models_deposited, localization_densities,
1443 """Add an ensemble generated by ad hoc methods (not using PMI).
1444 This is currently only used by the Nup84 system."""
1446 state = self._last_state
1447 group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1448 state.add_model_group(group)
1450 self.system.locations.append(ensemble_file)
1451 e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1453 self.system.ensembles.append(e)
1454 for c
in state.all_modeled_components:
1455 den = localization_densities.get(c,
None)
1457 e.load_localization_density(state, c, self.asym_units[c], den)
1461 """Point a previously-created ensemble to an 'all-models' file.
1462 This could be a trajectory such as DCD, an RMF, or a multimodel
1464 self.system.locations.append(location)
1466 ind = i + self._state_ensemble_offset
1467 self.system.ensembles[ind].file = location
1469 def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1475 protocol = self.all_protocols.get_last_protocol(state)
1476 num_models = protocol.steps[-1].num_models_end
1477 pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1478 pp.software = self.all_software.pmi
1479 self.all_protocols.add_postproc(pp, state)
1480 for i
in range(rex._number_of_clusters):
1481 group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1482 'cluster %d' % (i + 1)))
1483 state.add_model_group(group)
1485 e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1486 self.system.ensembles.append(e)
1488 for fname, stuple
in sorted(density_custom_ranges.items()):
1489 e.load_localization_density(state, fname, stuple,
1491 for stats
in e.load_all_models(self, state):
1492 m = self.add_model(group)
1495 m.name =
'Best scoring model'
1498 for c
in state.all_modeled_components:
1501 def _get_subassembly(self, comps, name, description):
1502 """Get an Assembly consisting of the given components.
1503 `compdict` is a dictionary of the components to add, where keys
1504 are the component names and values are the sequence ranges (or
1505 None to use all residues in the component)."""
1507 for comp, seqrng
in comps.items():
1508 a = self.asym_units[comp]
1509 asyms.append(a
if seqrng
is None else a(*seqrng))
1511 a = ihm.Assembly(asyms, name=name, description=description)
1514 def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1516 """Add a basic FoXS fit. This is largely intended for use from the
1518 assembly = self._get_subassembly({comp:seqrange},
1519 name=
"SAXS subassembly",
1520 description=
"All components that fit SAXS data")
1521 r = ihm.restraint.SASRestraint(dataset, assembly, segment=
False,
1522 fitting_method=
'FoXS', fitting_atom_type=
'Heavy atoms',
1523 multi_state=
False, radius_of_gyration=rg, details=details)
1524 r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1525 self.system.restraints.append(r)
1526 self._add_restraint_dataset(r)
1528 def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1529 image_resolution, projection_number,
1530 micrographs_number):
1531 r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1532 image_resolution, projection_number,
1534 self.system.restraints.append(r)
1535 self._add_restraint_dataset(r)
1537 def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1539 r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1540 self.system.restraints.append(r)
1541 self._add_restraint_dataset(r)
1543 def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1544 sigma, pmi_restraint):
1545 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1546 sigma, pmi_restraint, self._xy_plane)
1548 def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1549 sigma, pmi_restraint):
1550 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1551 sigma, pmi_restraint, self._xz_plane)
1553 def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1554 sigma, pmi_restraint):
1555 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1556 sigma, pmi_restraint, self._z_axis)
1558 def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1559 sigma, pmi_restraint, geom):
1560 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1561 r = _GeometricRestraint(self, state, pmi_restraint, geom,
1562 asym.get_feature(ps),
1563 ihm.restraint.LowerUpperBoundDistanceRestraint(
1564 lower_bound, upper_bound),
1566 self.system.restraints.append(r)
1567 self._add_restraint_dataset(r)
1569 def _get_membrane(self, tor_R, tor_r, tor_th):
1570 """Get an object representing a half-torus membrane"""
1571 if not hasattr(self,
'_seen_membranes'):
1572 self._seen_membranes = {}
1575 membrane_id = tuple(int(x * 100.)
for x
in (tor_R, tor_r, tor_th))
1576 if membrane_id
not in self._seen_membranes:
1577 m = ihm.geometry.HalfTorus(center=self._center_origin,
1578 transformation=self._identity_transform,
1579 major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1580 inner=
True, name=
'Membrane')
1581 self._seen_membranes[membrane_id] = m
1582 return self._seen_membranes[membrane_id]
1584 def add_membrane_surface_location_restraint(
1585 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1586 self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1587 pmi_restraint, ihm.restraint.UpperBoundDistanceRestraint(0.))
1589 def add_membrane_exclusion_restraint(
1590 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1591 self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1592 pmi_restraint, ihm.restraint.LowerBoundDistanceRestraint(0.))
1594 def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1595 sigma, pmi_restraint, rsr):
1596 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1597 r = _GeometricRestraint(self, state, pmi_restraint,
1598 self._get_membrane(tor_R, tor_r, tor_th),
1599 asym.get_feature(ps), rsr, sigma)
1600 self.system.restraints.append(r)
1601 self._add_restraint_dataset(r)
1603 def add_model(self, group, assembly=None, representation=None):
1604 state = self._last_state
1605 if representation
is None:
1606 representation = self.default_representation
1607 protocol = self.all_protocols.get_last_protocol(state)
1608 m = _Model(state.prot, self, protocol,
1609 assembly
if assembly
else state.modeled_assembly,
1614 def _update_locations(self, filelocs):
1615 """Update FileLocation to point to a parent repository, if any"""
1616 all_repos = [m
for m
in self._metadata
1617 if isinstance(m, ihm.location.Repository)]
1618 for fileloc
in filelocs:
1619 ihm.location.Repository._update_in_repos(fileloc, all_repos)
1621 _metadata = property(
lambda self:
1622 itertools.chain.from_iterable(self._each_metadata))
1625 def read(fh, model_class=ihm.model.Model, format='mmCIF', handlers=[]):
1626 """Read data from the mmCIF file handle `fh`.
1627 This is a simple wrapper around `ihm.reader.read()`, which also
1628 reads PMI-specific information from the mmCIF or BinaryCIF file."""
1629 return ihm.reader.read(fh, model_class=model_class, format=format,
1630 handlers=[_ReplicaExchangeProtocolHandler]+handlers)
1634 """Extract metadata from an EM density GMM file."""
1637 """Extract metadata from `filename`.
1638 @return a dict with key `dataset` pointing to the GMM file and
1639 `number_of_gaussians` to the number of GMMs (or None)"""
1640 l = ihm.location.InputFileLocation(filename,
1641 details=
"Electron microscopy density map, "
1642 "represented as a Gaussian Mixture Model (GMM)")
1645 l._allow_duplicates =
True
1646 d = ihm.dataset.EMDensityDataset(l)
1647 ret = {
'dataset':d,
'number_of_gaussians':
None}
1649 with open(filename)
as fh:
1651 if line.startswith(
'# data_fn: '):
1652 p = ihm.metadata.MRCParser()
1653 fn = line[11:].rstrip(
'\r\n')
1654 dataset = p.parse_file(os.path.join(
1655 os.path.dirname(filename), fn))[
'dataset']
1656 ret[
'dataset'].parents.append(dataset)
1657 elif line.startswith(
'# ncenters: '):
1658 ret[
'number_of_gaussians'] = int(line[12:])
Select non water and non hydrogen atoms.
def create_representation
Create a new Representation and return it.
def create_transformed_component
Make a new component that's a transformed copy of another.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def exclude_coordinates
Don't record coordinates for the given domain.
A decorator to associate a particle with a part of a protein/DNA/RNA.
Class to encode a modeling protocol as mmCIF.
def flush
Write out all information to the file.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def parse_file
Extract metadata from filename.
def pmi_range
Return a range of IHM residues indexed using PMI numbering.
Extract metadata from an EM density GMM file.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def pmi_residue
Return a single IHM residue indexed using PMI numbering.
void read_pdb(TextInput input, int model, Hierarchy h)
A single asymmetric unit in the system.
A single entity in the system.
Representation of the system.
def set_ensemble_file
Point a previously-created ensemble to an 'all-models' file.
Classes to represent data structures used in mmCIF.
def pmi_range
Return a range of IHM residues indexed using PMI numbering.
void add_restraint(RMF::FileHandle fh, Restraint *hs)
static bool get_is_setup(Model *m, ParticleIndex pi)
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def finalize
Do any final processing on the class hierarchy.
Base class for capturing a modeling protocol.
Basic utilities for handling cryo-electron microscopy 3D density maps.
Class for easy writing of PDBs, RMFs, and stat files.
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Classes for writing output files and processing them.
A decorator for a residue.
Basic functionality that is expected to be used by a wide variety of IMP users.
def pmi_residue
Return a single IHM residue indexed using PMI numbering.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
def read
Read data from the mmCIF file handle fh.
Functionality for loading, creating, manipulating and scoring atomic structures.
Hierarchies get_leaves(const Selection &h)
Select hierarchy particles identified by the biological name.
Select all ATOM and HETATM records with the given chain ids.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...