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)
947 if os.path.exists(rmf_file):
948 rh = RMF.open_rmf_file_read_only(rmf_file)
949 system = state._pmi_object.system
955 if model_num >= self.num_models_deposited:
959 def _get_precision(self):
960 precfile = os.path.join(self.post_process.rex._outputdir,
961 "precision.%d.%d.out" % (self.cluster_num,
963 if not os.path.exists(precfile):
966 r = re.compile(
'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
968 with open(precfile)
as fh:
972 return float(m.group(1))
974 precision = property(
lambda self: self._get_precision(),
975 lambda self, val:
None)
977 class _SimpleEnsemble(ihm.model.Ensemble):
978 """Simple manually-created ensemble"""
980 num_models_deposited =
None
982 def __init__(self, pp, model_group, num_models, drmsd,
983 num_models_deposited, ensemble_file):
984 super(_SimpleEnsemble, self).__init__(
985 model_group=model_group, post_process=pp, num_models=num_models,
986 file=ensemble_file, precision=drmsd, name=model_group.name,
987 clustering_feature=
'dRMSD')
988 self.num_models_deposited = num_models_deposited
990 def load_localization_density(self, state, component, asym, local_file):
991 den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
992 self.densities.append(den)
995 class _EntityMapper(dict):
996 """Handle mapping from IMP components (without copy number) to CIF entities.
997 Multiple components may map to the same entity if they share sequence."""
998 def __init__(self, system):
999 super(_EntityMapper, self).__init__()
1000 self._sequence_dict = {}
1002 self.system = system
1004 def add(self, component_name, sequence, offset):
1005 def entity_seq(sequence):
1008 return [
'UNK' if s ==
'X' else s
for s
in sequence]
1011 if sequence
not in self._sequence_dict:
1014 d = component_name.split(
"@")[0].split(
".")[0]
1015 entity = Entity(entity_seq(sequence), description=d,
1017 self.system.entities.append(entity)
1018 self._sequence_dict[sequence] = entity
1019 self[component_name] = self._sequence_dict[sequence]
1022 class _TransformedComponent(object):
1023 def __init__(self, name, original, transform):
1024 self.name, self.original, self.transform = name, original, transform
1027 class _SimpleRef(object):
1028 """Class with similar interface to weakref.ref, but keeps a strong ref"""
1029 def __init__(self, ref):
1035 class _State(ihm.model.State):
1036 """Representation of a single state in the system."""
1038 def __init__(self, pmi_object, po):
1043 self._pmi_object = weakref.proxy(pmi_object)
1044 if hasattr(pmi_object,
'state'):
1047 self._pmi_state = _SimpleRef(pmi_object.state)
1049 self._pmi_state = weakref.ref(pmi_object)
1051 old_name = self.name
1052 super(_State, self).__init__(experiment_type=
'Fraction of bulk')
1053 self.name = old_name
1057 self.modeled_assembly = ihm.Assembly(
1058 name=
"Modeled assembly",
1059 description=self.get_postfixed_name(
1060 "All components modeled by IMP"))
1061 po.system.orphan_assemblies.append(self.modeled_assembly)
1063 self.all_modeled_components = []
1066 return hash(self._pmi_state())
1067 def __eq__(self, other):
1068 return self._pmi_state() == other._pmi_state()
1070 def add_model_group(self, group):
1074 def get_prefixed_name(self, name):
1075 """Prefix the given name with the state name, if available."""
1077 return self.short_name +
' ' + name
1081 return name[0].upper() + name[1:]
if name
else ''
1083 def get_postfixed_name(self, name):
1084 """Postfix the given name with the state name, if available."""
1086 return "%s in state %s" % (name, self.short_name)
1090 short_name = property(
lambda self: self._pmi_state().short_name)
1091 long_name = property(
lambda self: self._pmi_state().long_name)
1092 def __get_name(self):
1093 return self._pmi_state().long_name
1094 def __set_name(self, val):
1095 self._pmi_state().long_name = val
1096 name = property(__get_name, __set_name)
1100 """A single entity in the system.
1101 This functions identically to the base ihm.Entity class, but it
1102 allows identifying residues by either the IHM numbering scheme
1103 (seq_id, which is always contiguous starting at 1) or by the PMI
1104 scheme (which is similar, but may not start at 1 if the sequence in
1105 the FASTA file has one or more N-terminal gaps"""
1106 def __init__(self, sequence, pmi_offset, *args, **keys):
1109 self.pmi_offset = pmi_offset
1110 super(Entity, self).__init__(sequence, *args, **keys)
1113 """Return a single IHM residue indexed using PMI numbering"""
1114 return self.residue(res_id - self.pmi_offset)
1117 """Return a range of IHM residues indexed using PMI numbering"""
1118 off = self.pmi_offset
1119 return self(res_id_begin - off, res_id_end - off)
1123 """A single asymmetric unit in the system.
1124 This functions identically to the base ihm.AsymUnit class, but it
1125 allows identifying residues by either the IHM numbering scheme
1126 (seq_id, which is always contiguous starting at 1) or by the PMI
1127 scheme (which is similar, but may not start at 1 if the sequence in
1128 the FASTA file has one or more N-terminal gaps"""
1130 def __init__(self, entity, *args, **keys):
1131 super(AsymUnit, self).__init__(
1132 entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1135 """Return a single IHM residue indexed using PMI numbering"""
1136 return self.residue(res_id - self.entity.pmi_offset)
1139 """Return a range of IHM residues indexed using PMI numbering"""
1140 off = self.entity.pmi_offset
1141 return self(res_id_begin - off, res_id_end - off)
1145 """Class to encode a modeling protocol as mmCIF.
1147 IMP has basic support for writing out files in mmCIF format, for
1148 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1149 After creating an instance of this class, attach it to an
1150 IMP.pmi.representation.Representation object. After this, any
1151 generated models and metadata are automatically collected and
1154 def __init__(self, fh):
1156 self.system = ihm.System()
1157 self._state_group = ihm.model.StateGroup()
1158 self.system.state_groups.append(self._state_group)
1161 self._state_ensemble_offset = 0
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(name=name)
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()
1326 def add_pdb_element(self, state, name, start, end, offset, pdbname,
1327 chain, hier, representation=
None):
1328 if self._is_excluded(name, start, end):
1330 if representation
is None:
1331 representation = self.default_representation
1332 asym = self.asym_units[name]
1333 p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1335 self.all_representations.add_fragment(state, representation, p)
1336 self.all_starting_models.add_pdb_fragment(p)
1338 def add_bead_element(self, state, name, start, end, num, hier,
1339 representation=
None):
1340 if self._is_excluded(name, start, end):
1342 if representation
is None:
1343 representation = self.default_representation
1344 asym = self.asym_units[name]
1345 pmi_offset = asym.entity.pmi_offset
1346 b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1348 self.all_representations.add_fragment(state, representation, b)
1350 def get_cross_link_group(self, pmi_restraint):
1351 r = _CrossLinkRestraint(pmi_restraint)
1352 self.system.restraints.append(r)
1353 self._add_restraint_dataset(r)
1356 def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1357 if c1
not in self._all_components
or c2
not in self._all_components:
1363 e1 = self.entities[c1]
1364 e2 = self.entities[c2]
1365 xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1366 residue2=e2.pmi_residue(r2))
1367 rsr.experimental_cross_links.append([xl])
1370 def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1373 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1374 d = ihm.restraint.UpperBoundDistanceRestraint(length)
1376 if _get_by_residue(p1)
and _get_by_residue(p2):
1377 cls = _ResidueCrossLink
1379 cls = _FeatureCrossLink
1380 xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1383 xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1384 rsr.cross_links.append(xl)
1386 def add_replica_exchange(self, state, rex):
1391 step = _ReplicaExchangeProtocolStep(state, rex)
1392 step.software = self.all_software.pmi
1393 self.all_protocols.add_step(step, state)
1395 def _add_simple_dynamics(self, num_models_end, method):
1397 state = self._last_state
1398 self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1401 def _add_protocol(self):
1403 state = self._last_state
1404 self.all_protocols.add_protocol(state)
1406 def _add_dataset(self, dataset):
1407 return self.all_datasets.add(self._last_state, dataset)
1409 def _add_restraint_dataset(self, restraint):
1410 return self.all_datasets.add_restraint(self._last_state, restraint)
1412 def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1414 state = self._last_state
1415 pp = ihm.analysis.ClusterStep(
'RMSD', num_models_begin, num_models_end)
1416 self.all_protocols.add_postproc(pp, state)
1419 def _add_no_postprocessing(self, num_models):
1421 state = self._last_state
1422 pp = ihm.analysis.EmptyStep()
1423 pp.num_models_begin = pp.num_models_end = num_models
1424 self.all_protocols.add_postproc(pp, state)
1427 def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1428 num_models_deposited, localization_densities,
1430 """Add an ensemble generated by ad hoc methods (not using PMI).
1431 This is currently only used by the Nup84 system."""
1433 state = self._last_state
1434 group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1435 state.add_model_group(group)
1437 self.system.locations.append(ensemble_file)
1438 e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1440 self.system.ensembles.append(e)
1441 for c
in state.all_modeled_components:
1442 den = localization_densities.get(c,
None)
1444 e.load_localization_density(state, c, self.asym_units[c], den)
1448 """Point a previously-created ensemble to an 'all-models' file.
1449 This could be a trajectory such as DCD, an RMF, or a multimodel
1451 self.system.locations.append(location)
1453 ind = i + self._state_ensemble_offset
1454 self.system.ensembles[ind].file = location
1456 def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1462 protocol = self.all_protocols.get_last_protocol(state)
1463 num_models = protocol.steps[-1].num_models_end
1464 pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1465 pp.software = self.all_software.pmi
1466 self.all_protocols.add_postproc(pp, state)
1467 for i
in range(rex._number_of_clusters):
1468 group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1469 'cluster %d' % (i + 1)))
1470 state.add_model_group(group)
1472 e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1473 self.system.ensembles.append(e)
1475 for fname, stuple
in sorted(density_custom_ranges.items()):
1476 e.load_localization_density(state, fname, stuple,
1478 for stats
in e.load_all_models(self, state):
1479 m = self.add_model(group)
1482 m.name =
'Best scoring model'
1485 for c
in state.all_modeled_components:
1488 def _get_subassembly(self, comps, name, description):
1489 """Get an Assembly consisting of the given components.
1490 `compdict` is a dictionary of the components to add, where keys
1491 are the component names and values are the sequence ranges (or
1492 None to use all residues in the component)."""
1494 for comp, seqrng
in comps.items():
1495 a = self.asym_units[comp]
1496 asyms.append(a
if seqrng
is None else a(*seqrng))
1498 a = ihm.Assembly(asyms, name=name, description=description)
1501 def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1503 """Add a basic FoXS fit. This is largely intended for use from the
1505 assembly = self._get_subassembly({comp:seqrange},
1506 name=
"SAXS subassembly",
1507 description=
"All components that fit SAXS data")
1508 r = ihm.restraint.SASRestraint(dataset, assembly, segment=
False,
1509 fitting_method=
'FoXS', fitting_atom_type=
'Heavy atoms',
1510 multi_state=
False, radius_of_gyration=rg, details=details)
1511 r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1512 self.system.restraints.append(r)
1513 self._add_restraint_dataset(r)
1515 def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1516 image_resolution, projection_number,
1517 micrographs_number):
1518 r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1519 image_resolution, projection_number,
1521 self.system.restraints.append(r)
1522 self._add_restraint_dataset(r)
1524 def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1526 r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1527 self.system.restraints.append(r)
1528 self._add_restraint_dataset(r)
1530 def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1531 sigma, pmi_restraint):
1532 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1533 sigma, pmi_restraint, self._xy_plane)
1535 def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1536 sigma, pmi_restraint):
1537 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1538 sigma, pmi_restraint, self._xz_plane)
1540 def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1541 sigma, pmi_restraint):
1542 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1543 sigma, pmi_restraint, self._z_axis)
1545 def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1546 sigma, pmi_restraint, geom):
1547 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1548 r = _GeometricRestraint(self, state, pmi_restraint, geom,
1549 asym.get_feature(ps),
1550 ihm.restraint.LowerUpperBoundDistanceRestraint(
1551 lower_bound, upper_bound),
1553 self.system.restraints.append(r)
1554 self._add_restraint_dataset(r)
1556 def _get_membrane(self, tor_R, tor_r, tor_th):
1557 """Get an object representing a half-torus membrane"""
1558 if not hasattr(self,
'_seen_membranes'):
1559 self._seen_membranes = {}
1562 membrane_id = tuple(int(x * 100.)
for x
in (tor_R, tor_r, tor_th))
1563 if membrane_id
not in self._seen_membranes:
1564 m = ihm.geometry.HalfTorus(center=self._center_origin,
1565 transformation=self._identity_transform,
1566 major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1567 inner=
True, name=
'Membrane')
1568 self._seen_membranes[membrane_id] = m
1569 return self._seen_membranes[membrane_id]
1571 def add_membrane_surface_location_restraint(
1572 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1573 self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1574 pmi_restraint, ihm.restraint.UpperBoundDistanceRestraint(0.))
1576 def add_membrane_exclusion_restraint(
1577 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1578 self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1579 pmi_restraint, ihm.restraint.LowerBoundDistanceRestraint(0.))
1581 def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1582 sigma, pmi_restraint, rsr):
1583 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1584 r = _GeometricRestraint(self, state, pmi_restraint,
1585 self._get_membrane(tor_R, tor_r, tor_th),
1586 asym.get_feature(ps), rsr, sigma)
1587 self.system.restraints.append(r)
1588 self._add_restraint_dataset(r)
1590 def add_model(self, group, assembly=None, representation=None):
1591 state = self._last_state
1592 if representation
is None:
1593 representation = self.default_representation
1594 protocol = self.all_protocols.get_last_protocol(state)
1595 m = _Model(state.prot, self, protocol,
1596 assembly
if assembly
else state.modeled_assembly,
1602 def read(fh, model_class=ihm.model.Model, format='mmCIF', handlers=[]):
1603 """Read data from the mmCIF file handle `fh`.
1604 This is a simple wrapper around `ihm.reader.read()`, which also
1605 reads PMI-specific information from the mmCIF or BinaryCIF file."""
1606 return ihm.reader.read(fh, model_class=model_class, format=format,
1607 handlers=[_ReplicaExchangeProtocolHandler]+handlers)
1611 """Extract metadata from an EM density GMM file."""
1614 """Extract metadata from `filename`.
1615 @return a dict with key `dataset` pointing to the GMM file and
1616 `number_of_gaussians` to the number of GMMs (or None)"""
1617 l = ihm.location.InputFileLocation(filename,
1618 details=
"Electron microscopy density map, "
1619 "represented as a Gaussian Mixture Model (GMM)")
1622 l._allow_duplicates =
True
1623 d = ihm.dataset.EMDensityDataset(l)
1624 ret = {
'dataset':d,
'number_of_gaussians':
None}
1626 with open(filename)
as fh:
1628 if line.startswith(
'# data_fn: '):
1629 p = ihm.metadata.MRCParser()
1630 fn = line[11:].rstrip(
'\r\n')
1631 dataset = p.parse_file(os.path.join(
1632 os.path.dirname(filename), fn))[
'dataset']
1633 ret[
'dataset'].parents.append(dataset)
1634 elif line.startswith(
'# ncenters: '):
1635 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.
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.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
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.
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
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 ...