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
42 import ihm.representation
44 import ihm.cross_linkers
46 def _assign_id(obj, seen_objs, obj_by_id):
47 """Assign a unique ID to obj, and track all ids in obj_by_id."""
48 if obj
not in seen_objs:
49 if not hasattr(obj,
'id'):
51 obj.id = len(obj_by_id)
52 seen_objs[obj] = obj.id
54 obj.id = seen_objs[obj]
57 def _get_by_residue(p):
58 """Determine whether the given particle represents a specific residue
59 or a more coarse-grained object."""
63 class _ComponentMapper(object):
64 """Map a Particle to a component name"""
65 def __init__(self, prot):
68 self.name =
'cif-output'
69 self.o.dictionary_pdbs[self.name] = self.prot
70 self.o._init_dictchain(self.name, self.prot, multichar_chain=
True)
72 def __getitem__(self, p):
73 protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
77 class _AsymMapper(object):
78 """Map a Particle to an asym_unit"""
79 def __init__(self, simo, prot):
81 self._cm = _ComponentMapper(prot)
82 self._seen_ranges = {}
84 def __getitem__(self, p):
85 protname = self._cm[p]
86 return self.simo.asym_units[protname]
88 def get_feature(self, ps):
89 """Get an ihm.restraint.Feature that covers the given particles"""
97 rng = asym(rind, rind)
101 rng = asym(rinds[0], rinds[-1])
103 raise ValueError(
"Unsupported particle type %s" % str(p))
105 if len(rngs) > 0
and rngs[-1].asym == asym \
106 and rngs[-1].seq_id_range[1] == rng.seq_id_range[0] - 1:
107 rngs[-1].seq_id_range = (rngs[-1].seq_id_range[0],
114 if hrngs
in self._seen_ranges:
115 return self._seen_ranges[hrngs]
117 feat = ihm.restraint.ResidueFeature(rngs)
118 self._seen_ranges[hrngs] = feat
122 class _AllSoftware(object):
123 def __init__(self, system):
125 self.modeller_used = self.phyre2_used =
False
126 self.pmi = ihm.Software(
127 name=
"IMP PMI module",
128 version=IMP.pmi.__version__,
129 classification=
"integrative model building",
130 description=
"integrative model building",
131 location=
'https://integrativemodeling.org')
132 self.imp = ihm.Software(
133 name=
"Integrative Modeling Platform (IMP)",
134 version=IMP.__version__,
135 classification=
"integrative model building",
136 description=
"integrative model building",
137 location=
'https://integrativemodeling.org')
138 self.system.software.extend([self.pmi, self.imp])
140 def set_modeller_used(self, version, date):
141 if self.modeller_used:
143 self.modeller_used =
True
144 self.system.software.append(ihm.Software(
145 name=
'MODELLER', classification=
'comparative modeling',
146 description=
'Comparative modeling by satisfaction '
147 'of spatial restraints, build ' + date,
148 location=
'https://salilab.org/modeller/',
151 def set_phyre2_used(self):
154 self.phyre2_used =
True
155 self.system.software.append(ihm.Software(
156 name=
'Phyre2', classification=
'protein homology modeling',
157 description=
'Protein Homology/analogY Recognition '
160 location=
'http://www.sbg.bio.ic.ac.uk/~phyre2/'))
162 def _get_fragment_is_rigid(fragment):
163 """Determine whether a fragment is modeled rigidly"""
168 class _PDBFragment(ihm.representation.ResidueSegment):
169 """Record details about part of a PDB file used as input
171 def __init__(self, state, component, start, end, pdb_offset,
172 pdbname, chain, hier, asym_unit):
174 super(_PDBFragment, self).__init__(
175 asym_unit=asym_unit.pmi_range(start, end),
176 rigid=
None, primitive=
'sphere')
177 self.component, self.start, self.end, self.offset, self.pdbname \
178 = component, start, end, pdb_offset, pdbname
179 self.state, self.chain, self.hier = state, chain, hier
184 rigid = property(
lambda self: _get_fragment_is_rigid(self),
185 lambda self, val:
None)
187 def combine(self, other):
191 class _BeadsFragment(ihm.representation.FeatureSegment):
192 """Record details about beads used to represent part of a component."""
194 def __init__(self, state, component, start, end, count, hier, asym_unit):
195 super(_BeadsFragment, self).__init__(asym_unit=asym_unit(start, end),
196 rigid=
None, primitive=
'sphere', count=count)
197 self.state, self.component, self.hier \
198 = state, component, hier
200 rigid = property(
lambda self: _get_fragment_is_rigid(self),
201 lambda self, val:
None)
203 def combine(self, other):
205 if type(other) == type(self)
and \
206 other.asym_unit.seq_id_range[0] == self.asym_unit.seq_id_range[1] + 1:
207 self.asym_unit.seq_id_range = (self.asym_unit.seq_id_range[0],
208 other.asym_unit.seq_id_range[1])
209 self.count += other.count
213 class _AllModelRepresentations(object):
214 def __init__(self, simo):
218 self.fragments = OrderedDict()
219 self._all_representations = {}
221 def copy_component(self, state, name, original, asym_unit):
222 """Copy all representation for `original` in `state` to `name`"""
225 newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
227 for rep
in self.fragments:
228 if original
in self.fragments[rep]:
229 if name
not in self.fragments[rep]:
230 self.fragments[rep][name] = OrderedDict()
231 self.fragments[rep][name][state] = [copy_frag(f)
232 for f
in self.fragments[rep][original][state]]
235 first_state = list(self.fragments[rep][name].keys())[0]
236 if state
is first_state:
237 representation = self._all_representations[rep]
238 representation.extend(self.fragments[rep][name][state])
240 def add_fragment(self, state, representation, fragment):
241 """Add a model fragment."""
242 comp = fragment.component
243 id_rep = id(representation)
244 self._all_representations[id_rep] = representation
245 if id_rep
not in self.fragments:
246 self.fragments[id_rep] = OrderedDict()
247 if comp
not in self.fragments[id_rep]:
248 self.fragments[id_rep][comp] = OrderedDict()
249 if state
not in self.fragments[id_rep][comp]:
250 self.fragments[id_rep][comp][state] = []
251 fragments = self.fragments[id_rep][comp][state]
252 if len(fragments) == 0
or not fragments[-1].combine(fragment):
253 fragments.append(fragment)
256 first_state = list(self.fragments[id_rep][comp].keys())[0]
257 if state
is first_state:
258 representation.append(fragment)
261 class _AllDatasets(object):
262 """Track all datasets generated by PMI and add them to the ihm.System"""
263 def __init__(self, system):
266 self._datasets_by_state = {}
267 self._restraints_by_state = {}
269 def get_all_group(self, state):
270 """Get a DatasetGroup encompassing all datasets so far in this state"""
274 g = ihm.dataset.DatasetGroup(self._datasets_by_state.get(state, [])
275 + [r.dataset
for r
in self._restraints_by_state.get(state, [])
279 def add(self, state, dataset):
280 """Add a new dataset."""
281 self._datasets.append(dataset)
282 if state
not in self._datasets_by_state:
283 self._datasets_by_state[state] = []
284 self._datasets_by_state[state].append(dataset)
287 self.system.orphan_datasets.append(dataset)
291 """Add the dataset for a restraint"""
292 if state
not in self._restraints_by_state:
293 self._restraints_by_state[state] = []
294 self._restraints_by_state[state].append(restraint)
297 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
298 """Restrain to a set of cross-links"""
301 _label_map = {
'wtDSS':
'DSS',
'scDSS':
'DSS',
'scEDC':
'EDC'}
302 _descriptor_map = {
'DSS': ihm.cross_linkers.dss,
303 'EDC': ihm.cross_linkers.edc}
305 def __init__(self, pmi_restraint):
306 self.pmi_restraint = pmi_restraint
309 linker = getattr(self.pmi_restraint,
'linker',
None)
310 label = self.pmi_restraint.label
312 super(_CrossLinkRestraint, self).__init__(
313 dataset=self.pmi_restraint.dataset,
314 linker=linker
or self._get_chem_descriptor(label))
317 def _get_chem_descriptor(cls, label):
319 label = cls._label_map.get(label, label)
320 if label
not in cls._descriptor_map:
324 d = ihm.ChemDescriptor(label)
325 cls._descriptor_map[label] = d
326 return cls._descriptor_map[label]
328 def _set_psi_sigma(self, model):
331 if model.m != self.pmi_restraint.model:
333 for resolution
in self.pmi_restraint.sigma_dictionary:
334 statname =
'ISDCrossLinkMS_Sigma_%s_%s' % (resolution, self.label)
335 if model.stats
and statname
in model.stats:
336 sigma = float(model.stats[statname])
337 p = self.pmi_restraint.sigma_dictionary[resolution][0]
338 old_values.append((p, p.get_scale()))
340 for psiindex
in self.pmi_restraint.psi_dictionary:
341 statname =
'ISDCrossLinkMS_Psi_%s_%s' % (psiindex, self.label)
342 if model.stats
and statname
in model.stats:
343 psi = float(model.stats[statname])
344 p = self.pmi_restraint.psi_dictionary[psiindex][0]
345 old_values.append((p, p.get_scale()))
349 return list(reversed(old_values))
351 def add_fits_from_model_statfile(self, model):
353 old_values = self._set_psi_sigma(model)
357 for xl
in self.cross_links:
359 xl.fits[model] = ihm.restraint.CrossLinkFit(
360 psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
362 for p, scale
in old_values:
366 def __set_dataset(self, val):
367 self.pmi_restraint.dataset = val
368 dataset = property(
lambda self: self.pmi_restraint.dataset,
372 def get_asym_mapper_for_state(simo, state, asym_map):
373 asym = asym_map.get(state,
None)
375 asym = _AsymMapper(simo, state.prot)
376 asym_map[state] = asym
379 class _PMICrossLink(object):
382 psi = property(
lambda self: self.psi_p.get_scale(),
383 lambda self, val:
None)
384 sigma1 = property(
lambda self: self.sigma1_p.get_scale(),
385 lambda self, val:
None)
386 sigma2 = property(
lambda self: self.sigma2_p.get_scale(),
387 lambda self, val:
None)
390 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
394 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
398 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
399 def __init__(self, state, pmi_restraint, image_number, resolution,
400 pixel_size, image_resolution, projection_number,
402 self.pmi_restraint, self.image_number = pmi_restraint, image_number
403 super(_EM2DRestraint, self).__init__(
404 dataset=pmi_restraint.datasets[image_number],
405 assembly=state.modeled_assembly,
406 segment=
False, number_raw_micrographs=micrographs_number,
407 pixel_size_width=pixel_size, pixel_size_height=pixel_size,
408 image_resolution=image_resolution,
409 number_of_projections=projection_number)
412 def __get_dataset(self):
413 return self.pmi_restraint.datasets[self.image_number]
414 def __set_dataset(self, val):
415 self.pmi_restraint.datasets[self.image_number] = val
416 dataset = property(__get_dataset, __set_dataset)
418 def add_fits_from_model_statfile(self, model):
419 ccc = self._get_cross_correlation(model)
420 transform = self._get_transformation(model)
421 rot = transform.get_rotation()
422 rm = [[e
for e
in rot.get_rotation_matrix_row(i)]
for i
in range(3)]
423 self.fits[model] = ihm.restraint.EM2DRestraintFit(
424 cross_correlation_coefficient=ccc,
426 tr_vector=transform.get_translation())
428 def _get_transformation(self, model):
429 """Get the transformation that places the model on the image"""
430 stats = model.em2d_stats
or model.stats
431 prefix =
'ElectronMicroscopy2D_%s_Image%d' % (self.pmi_restraint.label,
432 self.image_number + 1)
433 r = [float(stats[prefix +
'_Rotation%d' % i])
for i
in range(4)]
434 t = [float(stats[prefix +
'_Translation%d' % i])
438 inv = model.transform.get_inverse()
442 def _get_cross_correlation(self, model):
443 """Get the cross correlation coefficient between the model projection
445 stats = model.em2d_stats
or model.stats
446 return float(stats[
'ElectronMicroscopy2D_%s_Image%d_CCC'
447 % (self.pmi_restraint.label,
448 self.image_number + 1)])
451 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
453 def __init__(self, simo, state, pmi_restraint, target_ps, densities):
454 self.pmi_restraint = pmi_restraint
455 super(_EM3DRestraint, self).__init__(
456 dataset=pmi_restraint.dataset,
457 assembly=self._get_assembly(densities, simo, state),
458 fitting_method=
'Gaussian mixture models',
459 number_of_gaussians=len(target_ps))
462 def __set_dataset(self, val):
463 self.pmi_restraint.dataset = val
464 dataset = property(
lambda self: self.pmi_restraint.dataset,
467 def _get_assembly(self, densities, simo, state):
468 """Get the Assembly that this restraint acts on"""
469 cm = _ComponentMapper(state.prot)
472 components[cm[d]] =
None
473 a = simo._get_subassembly(components,
474 name=
"EM subassembly",
475 description=
"All components that fit the EM map")
478 def add_fits_from_model_statfile(self, model):
479 ccc = self._get_cross_correlation(model)
480 self.fits[model] = ihm.restraint.EM3DRestraintFit(
481 cross_correlation_coefficient=ccc)
483 def _get_cross_correlation(self, model):
484 """Get the cross correlation coefficient between the model
486 if model.stats
is not None:
487 return float(model.stats[
'GaussianEMRestraint_%s_CCC'
488 % self.pmi_restraint.label])
491 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
493 def __init__(self, simo, state, pmi_restraint, geometric_object,
494 feature, distance, sigma):
495 self.pmi_restraint = pmi_restraint
496 super(_GeometricRestraint, self).__init__(
497 dataset=pmi_restraint.dataset,
498 geometric_object=geometric_object, feature=feature,
499 distance=distance, harmonic_force_constant = 1. / sigma,
503 def __set_dataset(self, val):
504 self.pmi_restraint.dataset = val
505 dataset = property(
lambda self: self.pmi_restraint.dataset,
509 class _ReplicaExchangeProtocolStep(ihm.protocol.Step):
510 def __init__(self, state, rex):
511 if rex.monte_carlo_sample_objects
is not None:
512 method =
'Replica exchange monte carlo'
514 method =
'Replica exchange molecular dynamics'
515 self.monte_carlo_temperature = rex.vars[
'monte_carlo_temperature']
516 self.replica_exchange_minimum_temperature = \
517 rex.vars[
'replica_exchange_minimum_temperature']
518 self.replica_exchange_maximum_temperature = \
519 rex.vars[
'replica_exchange_maximum_temperature']
520 super(_ReplicaExchangeProtocolStep, self).__init__(
521 assembly=state.modeled_assembly,
523 method=method, name=
'Sampling',
524 num_models_begin=
None,
525 num_models_end=rex.vars[
"number_of_frames"],
526 multi_scale=
True, multi_state=
False, ordered=
False)
529 class _ReplicaExchangeProtocolDumper(ihm.dumper.Dumper):
530 """Write IMP-specific information about replica exchange to mmCIF.
531 Note that IDs will have already been assigned by python-ihm's
532 standard modeling protocol dumper."""
533 def dump(self, system, writer):
534 with writer.loop(
"_imp_replica_exchange_protocol",
535 [
"protocol_id",
"step_id",
"monte_carlo_temperature",
536 "replica_exchange_minimum_temperature",
537 "replica_exchange_maximum_temperature"])
as l:
538 for p
in system._all_protocols():
540 if isinstance(s, _ReplicaExchangeProtocolStep):
541 self._dump_step(p, s, l)
543 def _dump_step(self, p, s, l):
544 l.write(protocol_id=p._id, step_id=s._id,
545 monte_carlo_temperature=s.monte_carlo_temperature,
546 replica_exchange_minimum_temperature= \
547 s.replica_exchange_minimum_temperature,
548 replica_exchange_maximum_temperature= \
549 s.replica_exchange_maximum_temperature)
552 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
553 category =
'_imp_replica_exchange_protocol'
555 """Read IMP-specific information about replica exchange from mmCIF."""
556 def __call__(self, protocol_id, step_id, monte_carlo_temperature,
557 replica_exchange_minimum_temperature,
558 replica_exchange_maximum_temperature):
559 p = self.sysr.protocols.get_by_id(protocol_id)
561 s = p.steps[int(step_id)-1]
563 s.__class__ = _ReplicaExchangeProtocolStep
564 s.monte_carlo_temperature = \
565 self.get_float(monte_carlo_temperature)
566 s.replica_exchange_minimum_temperature = \
567 self.get_float(replica_exchange_minimum_temperature)
568 s.replica_exchange_maximum_temperature = \
569 self.get_float(replica_exchange_maximum_temperature)
572 class _SimpleProtocolStep(ihm.protocol.Step):
573 def __init__(self, state, num_models_end, method):
574 super(_SimpleProtocolStep, self).__init__(
575 assembly=state.modeled_assembly,
577 method=method, name=
'Sampling',
578 num_models_begin=
None,
579 num_models_end=num_models_end,
580 multi_scale=
True, multi_state=
False, ordered=
False)
583 class _Chain(object):
584 """Represent a single chain in a Model"""
585 def __init__(self, pmi_chain_id, asym_unit):
586 self.pmi_chain_id, self.asym_unit = pmi_chain_id, asym_unit
590 def add(self, xyz, atom_type, residue_type, residue_index,
591 all_indexes, radius):
592 if atom_type
is None:
593 self.spheres.append((xyz, residue_type, residue_index,
594 all_indexes, radius))
596 self.atoms.append((xyz, atom_type, residue_type, residue_index,
597 all_indexes, radius))
598 orig_comp = property(
lambda self: self.comp)
600 class _TransformedChain(object):
601 """Represent a chain that is a transformed version of another"""
602 def __init__(self, orig_chain, asym_unit, transform):
603 self.orig_chain, self.asym_unit = orig_chain, asym_unit
604 self.transform = transform
606 def __get_spheres(self):
607 for (xyz, residue_type, residue_index, all_indexes,
608 radius)
in self.orig_chain.spheres:
609 yield (self.transform * xyz, residue_type, residue_index,
611 spheres = property(__get_spheres)
613 def __get_atoms(self):
614 for (xyz, atom_type, residue_type, residue_index, all_indexes,
615 radius)
in self.orig_chain.atoms:
616 yield (self.transform * xyz, atom_type, residue_type, residue_index,
618 atoms = property(__get_atoms)
620 entity = property(
lambda self: self.orig_chain.entity)
621 orig_comp = property(
lambda self: self.orig_chain.comp)
623 class _Excluder(object):
624 def __init__(self, component, simo):
625 self._seqranges = simo._exclude_coords.get(component, [])
627 def is_excluded(self, indexes):
628 """Return True iff the given sequence range is excluded."""
629 for seqrange
in self._seqranges:
630 if indexes[0] >= seqrange[0]
and indexes[-1] <= seqrange[1]:
634 class _Model(ihm.model.Model):
635 def __init__(self, prot, simo, protocol, assembly, representation):
636 super(_Model, self).__init__(assembly=assembly, protocol=protocol,
637 representation=representation)
638 self.simo = weakref.proxy(simo)
644 self.em2d_stats =
None
647 self._is_restrained =
True
650 self.m = prot.get_model()
651 o.dictionary_pdbs[name] = prot
652 o._init_dictchain(name, prot, multichar_chain=
True)
653 (particle_infos_for_pdb,
654 self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
656 self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
659 def all_chains(self, simo):
660 """Yield all chains, including transformed ones"""
662 for c
in self.chains:
664 chain_for_comp[c.comp] = c
665 for tc
in simo._transformed_components:
666 orig_chain = chain_for_comp.get(tc.original,
None)
668 asym = simo.asym_units[tc.name]
669 c = _TransformedChain(orig_chain, asym, tc.transform)
673 def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
674 entity_for_chain = {}
677 for protname, chain_id
in o.dictchain[name].items():
678 if protname
in simo.entities:
679 entity_for_chain[chain_id] = simo.entities[protname]
682 pn = protname.split(
'.')[0]
683 entity_for_chain[chain_id] = simo.entities[pn]
684 comp_for_chain[chain_id] = protname
688 correct_asym[chain_id] = simo.asym_units[protname]
695 for (xyz, atom_type, residue_type, chain_id, residue_index,
696 all_indexes, radius)
in particle_infos_for_pdb:
697 if chain
is None or chain.pmi_chain_id != chain_id:
698 chain = _Chain(chain_id, correct_asym[chain_id])
699 chain.entity = entity_for_chain[chain_id]
700 chain.comp = comp_for_chain[chain_id]
701 self.chains.append(chain)
702 excluder = _Excluder(chain.comp, simo)
703 if not excluder.is_excluded(all_indexes
if all_indexes
704 else [residue_index]):
705 chain.add(xyz, atom_type, residue_type, residue_index,
708 def parse_rmsf_file(self, fname, component):
709 self.rmsf[component] = rmsf = {}
710 with open(fname)
as fh:
712 resnum, blocknum, val = line.split()
713 rmsf[int(resnum)] = (int(blocknum), float(val))
715 def get_rmsf(self, component, indexes):
716 """Get the RMSF value for the given residue indexes."""
719 rmsf = self.rmsf[component]
720 blocknums = dict.fromkeys(rmsf[ind][0]
for ind
in indexes)
721 if len(blocknums) != 1:
722 raise ValueError(
"Residue indexes %s aren't all in the same block"
724 return rmsf[indexes[0]][1]
727 for chain
in self.all_chains(self.simo):
728 pmi_offset = chain.asym_unit.entity.pmi_offset
729 for atom
in chain.atoms:
730 (xyz, atom_type, residue_type, residue_index,
731 all_indexes, radius) = atom
732 pt = self.transform * xyz
733 yield ihm.model.Atom(asym_unit=chain.asym_unit,
734 seq_id=residue_index - pmi_offset,
735 atom_id=atom_type.get_string(),
737 x=pt[0], y=pt[1], z=pt[2])
739 def get_spheres(self):
740 for chain
in self.all_chains(self.simo):
741 pmi_offset = chain.asym_unit.entity.pmi_offset
742 for sphere
in chain.spheres:
743 (xyz, residue_type, residue_index,
744 all_indexes, radius) = sphere
745 if all_indexes
is None:
746 all_indexes = (residue_index,)
747 pt = self.transform * xyz
748 yield ihm.model.Sphere(asym_unit=chain.asym_unit,
749 seq_id_range=(all_indexes[0] - pmi_offset,
750 all_indexes[-1] - pmi_offset),
751 x=pt[0], y=pt[1], z=pt[2], radius=radius,
752 rmsf=self.get_rmsf(chain.orig_comp, all_indexes))
755 class _AllProtocols(object):
756 def __init__(self, simo):
757 self.simo = weakref.proxy(simo)
759 self.protocols = OrderedDict()
761 def add_protocol(self, state):
762 """Add a new Protocol"""
763 if state
not in self.protocols:
764 self.protocols[state] = []
765 p = ihm.protocol.Protocol()
766 self.simo.system.orphan_protocols.append(p)
767 self.protocols[state].append(p)
769 def add_step(self, step, state):
770 """Add a ProtocolStep to the last Protocol of the given State"""
771 if state
not in self.protocols:
772 self.add_protocol(state)
773 protocol = self.get_last_protocol(state)
774 if len(protocol.steps) == 0:
775 step.num_models_begin = 0
777 step.num_models_begin = protocol.steps[-1].num_models_end
778 protocol.steps.append(step)
779 step.id = len(protocol.steps)
781 step.dataset_group = self.simo.all_datasets.get_all_group(state)
783 def add_postproc(self, step, state):
784 """Add a postprocessing step to the last protocol"""
785 protocol = self.get_last_protocol(state)
786 if len(protocol.analyses) == 0:
787 protocol.analyses.append(ihm.analysis.Analysis())
788 protocol.analyses[-1].steps.append(step)
790 def get_last_protocol(self, state):
791 """Return the most recently-added _Protocol"""
792 return self.protocols[state][-1]
795 class _AllStartingModels(object):
796 def __init__(self, simo):
800 self.models = OrderedDict()
803 def add_pdb_fragment(self, fragment):
804 """Add a starting model PDB fragment."""
805 comp = fragment.component
806 state = fragment.state
807 if comp
not in self.models:
808 self.models[comp] = OrderedDict()
809 if state
not in self.models[comp]:
810 self.models[comp][state] = []
811 models = self.models[comp][state]
812 if len(models) == 0 \
813 or models[-1].fragments[0].pdbname != fragment.pdbname:
814 model = self._add_model(fragment)
818 models[-1].fragments.append(weakref.proxy(fragment))
820 pmi_offset = models[-1].asym_unit.entity.pmi_offset
821 sid_begin = min(fragment.start + fragment.offset - pmi_offset,
822 models[-1].asym_unit.seq_id_range[0])
823 sid_end = max(fragment.end + fragment.offset - pmi_offset,
824 models[-1].asym_unit.seq_id_range[1])
825 models[-1].asym_unit = fragment.asym_unit.asym(sid_begin, sid_end)
826 fragment.starting_model = models[-1]
828 def _add_model(self, f):
829 parser = ihm.metadata.PDBParser()
830 r = parser.parse_file(f.pdbname)
832 self.simo._add_dataset(r[
'dataset'])
834 templates = r[
'templates'].get(f.chain, [])
837 self.simo.system.locations.append(t.alignment_file)
839 self.simo._add_dataset(t.dataset)
840 source = r[
'entity_source'].get(f.chain)
842 f.asym_unit.entity.source = source
843 pmi_offset = f.asym_unit.entity.pmi_offset
845 asym_unit=f.asym_unit.asym.pmi_range(f.start, f.end),
846 dataset=r[
'dataset'], asym_id=f.chain,
847 templates=templates, offset=f.offset + pmi_offset,
848 metadata=r[
'metadata'],
849 software=r[
'software'][0]
if r[
'software']
else None,
850 script_file=r[
'script'])
851 m.fragments = [weakref.proxy(f)]
855 class _StartingModel(ihm.startmodel.StartingModel):
856 def get_seq_dif(self):
860 pmi_offset = self.asym_unit.entity.pmi_offset
861 mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
863 for f
in self.fragments:
865 residue_indexes=list(range(f.start - f.offset,
866 f.end - f.offset + 1)))
867 for a
in mh.get_ihm_atoms(sel.get_selected_particles(),
868 f.offset - pmi_offset):
870 self._seq_dif = mh._seq_dif
873 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
874 """Post processing using AnalysisReplicaExchange0 macro"""
876 def __init__(self, rex, num_models_begin):
879 for fname
in self.get_all_stat_files():
880 with open(fname)
as fh:
881 num_models_end += len(fh.readlines())
882 super(_ReplicaExchangeAnalysisPostProcess, self).__init__(
883 feature=
'RMSD', num_models_begin=num_models_begin,
884 num_models_end=num_models_end)
886 def get_stat_file(self, cluster_num):
887 return os.path.join(self.rex._outputdir,
"cluster.%d" % cluster_num,
890 def get_all_stat_files(self):
891 for i
in range(self.rex._number_of_clusters):
892 yield self.get_stat_file(i)
895 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
896 """Ensemble generated using AnalysisReplicaExchange0 macro"""
898 num_models_deposited =
None
900 def __init__(self, pp, cluster_num, model_group, num_deposit):
901 with open(pp.get_stat_file(cluster_num))
as fh:
902 num_models = len(fh.readlines())
903 super(_ReplicaExchangeAnalysisEnsemble, self).__init__(
904 num_models=num_models,
905 model_group=model_group, post_process=pp,
906 clustering_feature=pp.feature,
907 name=model_group.name)
908 self.cluster_num = cluster_num
909 self.num_models_deposited = num_deposit
911 def get_rmsf_file(self, component):
912 return os.path.join(self.post_process.rex._outputdir,
913 'cluster.%d' % self.cluster_num,
914 'rmsf.%s.dat' % component)
916 def load_rmsf(self, model, component):
917 fname = self.get_rmsf_file(component)
918 if os.path.exists(fname):
919 model.parse_rmsf_file(fname, component)
921 def get_localization_density_file(self, fname):
922 return os.path.join(self.post_process.rex._outputdir,
923 'cluster.%d' % self.cluster_num,
926 def load_localization_density(self, state, fname, select_tuple, asym_units):
927 fullpath = self.get_localization_density_file(fname)
928 if os.path.exists(fullpath):
929 details =
"Localization density for %s %s" \
930 % (fname, self.model_group.name)
931 local_file = ihm.location.OutputFileLocation(fullpath,
933 for s
in select_tuple:
934 if isinstance(s, tuple)
and len(s) == 3:
935 asym = asym_units[s[2]].pmi_range(s[0], s[1])
938 den = ihm.model.LocalizationDensity(file=local_file,
940 self.densities.append(den)
942 def load_all_models(self, simo, state):
943 stat_fname = self.post_process.get_stat_file(self.cluster_num)
945 with open(stat_fname)
as fh:
946 stats = ast.literal_eval(fh.readline())
948 rmf_file = os.path.join(os.path.dirname(stat_fname),
949 "%d.rmf3" % model_num)
951 if os.path.exists(rmf_file):
952 rh = RMF.open_rmf_file_read_only(rmf_file)
953 system = state._pmi_object.system
959 if model_num >= self.num_models_deposited:
963 def _get_precision(self):
964 precfile = os.path.join(self.post_process.rex._outputdir,
965 "precision.%d.%d.out" % (self.cluster_num,
967 if not os.path.exists(precfile):
970 r = re.compile(
'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
972 with open(precfile)
as fh:
976 return float(m.group(1))
978 precision = property(
lambda self: self._get_precision(),
979 lambda self, val:
None)
981 class _SimpleEnsemble(ihm.model.Ensemble):
982 """Simple manually-created ensemble"""
984 num_models_deposited =
None
986 def __init__(self, pp, model_group, num_models, drmsd,
987 num_models_deposited, ensemble_file):
988 super(_SimpleEnsemble, self).__init__(
989 model_group=model_group, post_process=pp, num_models=num_models,
990 file=ensemble_file, precision=drmsd, name=model_group.name,
991 clustering_feature=
'dRMSD')
992 self.num_models_deposited = num_models_deposited
994 def load_localization_density(self, state, component, asym, local_file):
995 den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
996 self.densities.append(den)
999 class _CustomDNAAlphabet(object):
1000 """Custom DNA alphabet that maps A,C,G,T (rather than DA,DC,DG,DT
1001 as in python-ihm)"""
1002 _comps = dict([cc.code_canonical, cc]
1003 for cc
in ihm.DNAAlphabet._comps.values())
1006 class _EntityMapper(dict):
1007 """Handle mapping from IMP components (without copy number) to CIF entities.
1008 Multiple components may map to the same entity if they share sequence."""
1009 def __init__(self, system):
1010 super(_EntityMapper, self).__init__()
1011 self._sequence_dict = {}
1013 self.system = system
1015 def _get_alphabet(self, alphabet):
1016 """Map a PMI alphabet to an IHM alphabet"""
1020 alphabet_map = {
None: ihm.LPeptideAlphabet,
1021 IMP.pmi.alphabets.amino_acid: ihm.LPeptideAlphabet,
1022 IMP.pmi.alphabets.rna: ihm.RNAAlphabet,
1023 IMP.pmi.alphabets.dna: _CustomDNAAlphabet}
1024 if alphabet
in alphabet_map:
1025 return alphabet_map[alphabet]
1027 raise TypeError(
"Don't know how to handle %s" % alphabet)
1029 def add(self, component_name, sequence, offset, alphabet):
1030 def entity_seq(sequence):
1033 return [
'UNK' if s ==
'X' else s
for s
in sequence]
1036 if sequence
not in self._sequence_dict:
1039 d = component_name.split(
"@")[0].split(
".")[0]
1040 entity = Entity(entity_seq(sequence), description=d,
1042 alphabet=self._get_alphabet(alphabet))
1043 self.system.entities.append(entity)
1044 self._sequence_dict[sequence] = entity
1045 self[component_name] = self._sequence_dict[sequence]
1048 class _TransformedComponent(object):
1049 def __init__(self, name, original, transform):
1050 self.name, self.original, self.transform = name, original, transform
1053 class _SimpleRef(object):
1054 """Class with similar interface to weakref.ref, but keeps a strong ref"""
1055 def __init__(self, ref):
1061 class _State(ihm.model.State):
1062 """Representation of a single state in the system."""
1064 def __init__(self, pmi_object, po):
1069 self._pmi_object = weakref.proxy(pmi_object)
1070 if hasattr(pmi_object,
'state'):
1073 self._pmi_state = _SimpleRef(pmi_object.state)
1075 self._pmi_state = weakref.ref(pmi_object)
1077 old_name = self.name
1078 super(_State, self).__init__(experiment_type=
'Fraction of bulk')
1079 self.name = old_name
1083 self.modeled_assembly = ihm.Assembly(
1084 name=
"Modeled assembly",
1085 description=self.get_postfixed_name(
1086 "All components modeled by IMP"))
1087 po.system.orphan_assemblies.append(self.modeled_assembly)
1089 self.all_modeled_components = []
1092 return hash(self._pmi_state())
1093 def __eq__(self, other):
1094 return self._pmi_state() == other._pmi_state()
1096 def add_model_group(self, group):
1100 def get_prefixed_name(self, name):
1101 """Prefix the given name with the state name, if available."""
1103 return self.short_name +
' ' + name
1107 return name[0].upper() + name[1:]
if name
else ''
1109 def get_postfixed_name(self, name):
1110 """Postfix the given name with the state name, if available."""
1112 return "%s in state %s" % (name, self.short_name)
1116 short_name = property(
lambda self: self._pmi_state().short_name)
1117 long_name = property(
lambda self: self._pmi_state().long_name)
1118 def __get_name(self):
1119 return self._pmi_state().long_name
1120 def __set_name(self, val):
1121 self._pmi_state().long_name = val
1122 name = property(__get_name, __set_name)
1126 """A single entity in the system.
1127 This functions identically to the base ihm.Entity class, but it
1128 allows identifying residues by either the IHM numbering scheme
1129 (seq_id, which is always contiguous starting at 1) or by the PMI
1130 scheme (which is similar, but may not start at 1 if the sequence in
1131 the FASTA file has one or more N-terminal gaps"""
1132 def __init__(self, sequence, pmi_offset, *args, **keys):
1135 self.pmi_offset = pmi_offset
1136 super(Entity, self).__init__(sequence, *args, **keys)
1139 """Return a single IHM residue indexed using PMI numbering"""
1140 return self.residue(res_id - self.pmi_offset)
1143 """Return a range of IHM residues indexed using PMI numbering"""
1144 off = self.pmi_offset
1145 return self(res_id_begin - off, res_id_end - off)
1149 """A single asymmetric unit in the system.
1150 This functions identically to the base ihm.AsymUnit class, but it
1151 allows identifying residues by either the IHM numbering scheme
1152 (seq_id, which is always contiguous starting at 1) or by the PMI
1153 scheme (which is similar, but may not start at 1 if the sequence in
1154 the FASTA file has one or more N-terminal gaps"""
1156 def __init__(self, entity, *args, **keys):
1157 super(AsymUnit, self).__init__(
1158 entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1161 """Return a single IHM residue indexed using PMI numbering"""
1162 return self.residue(res_id - self.entity.pmi_offset)
1165 """Return a range of IHM residues indexed using PMI numbering"""
1166 off = self.entity.pmi_offset
1167 return self(res_id_begin - off, res_id_end - off)
1171 """Class to encode a modeling protocol as mmCIF.
1173 IMP has basic support for writing out files in mmCIF format, for
1174 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1175 After creating an instance of this class, attach it to an
1176 IMP.pmi.topology.System object. After this, any
1177 generated models and metadata are automatically collected in the
1178 `system` attribute, which is an
1179 [ihm.System](https://python-ihm.readthedocs.io/en/latest/main.html#ihm.System) object.
1180 Once the protocol is complete, call finalize() to make sure `system`
1181 contains everything, then use the
1182 [python-ihm API](https://python-ihm.readthedocs.io/en/latest/dumper.html#ihm.dumper.write)
1183 to write out files in mmCIF or BinaryCIF format.
1185 See also get_handlers(), get_dumpers().
1187 def __init__(self, fh=None):
1189 self.system = ihm.System()
1190 self._state_group = ihm.model.StateGroup()
1191 self.system.state_groups.append(self._state_group)
1196 "Passing a file handle to ProtocolOutput is deprecated. "
1197 "Create it with no arguments, then use the python-ihm API to "
1198 "output files (after calling finalize())")
1199 self._state_ensemble_offset = 0
1200 self._main_script = os.path.abspath(sys.argv[0])
1203 loc = ihm.location.WorkflowFileLocation(path=self._main_script,
1204 details=
"The main integrative modeling script")
1205 self.system.locations.append(loc)
1208 self.__asym_states = {}
1209 self._working_directory = os.getcwd()
1211 "Default representation")
1212 self.entities = _EntityMapper(self.system)
1214 self.asym_units = {}
1215 self._all_components = {}
1216 self.all_modeled_components = []
1217 self._transformed_components = []
1218 self.sequence_dict = {}
1221 self._xy_plane = ihm.geometry.XYPlane()
1222 self._xz_plane = ihm.geometry.XZPlane()
1223 self._z_axis = ihm.geometry.ZAxis()
1224 self._center_origin = ihm.geometry.Center(0,0,0)
1225 self._identity_transform = ihm.geometry.Transformation.identity()
1228 self._exclude_coords = {}
1230 self.all_representations = _AllModelRepresentations(self)
1231 self.all_protocols = _AllProtocols(self)
1232 self.all_datasets = _AllDatasets(self.system)
1233 self.all_starting_models = _AllStartingModels(self)
1235 self.all_software = _AllSoftware(self.system)
1238 """Create a new Representation and return it. This can be
1239 passed to add_model(), add_bead_element() or add_pdb_element()."""
1240 r = ihm.representation.Representation(name=name)
1241 self.system.orphan_representations.append(r)
1245 """Don't record coordinates for the given domain.
1246 Coordinates for the given domain (specified by a component name
1247 and a 2-element tuple giving the start and end residue numbers)
1248 will be excluded from the mmCIF file. This can be used to exclude
1249 parts of the structure that weren't well resolved in modeling.
1250 Any bead or residue that lies wholly within this range will be
1251 excluded. Multiple ranges for a given component can be excluded
1252 by calling this method multiple times."""
1253 if component
not in self._exclude_coords:
1254 self._exclude_coords[component] = []
1255 self._exclude_coords[component].append(seqrange)
1257 def _is_excluded(self, component, start, end):
1258 """Return True iff this chunk of sequence should be excluded"""
1259 for seqrange
in self._exclude_coords.get(component, ()):
1260 if start >= seqrange[0]
and end <= seqrange[1]:
1263 def _add_state(self, state):
1264 """Create a new state and return a pointer to it."""
1265 self._state_ensemble_offset = len(self.system.ensembles)
1266 s = _State(state, self)
1267 self._state_group.append(s)
1268 self._last_state = s
1271 def _get_chain_for_component(self, name, output):
1272 """Get the chain ID for a component, if any."""
1274 if name
in self.asym_units:
1275 return self.asym_units[name]._id
1280 def _get_assembly_comps(self, assembly):
1281 """Get the names of the components in the given assembly"""
1285 comps[ca.details] =
None
1289 """Make a new component that's a transformed copy of another.
1290 All representation for the existing component is copied to the
1292 assembly_comps = self._get_assembly_comps(state.modeled_assembly)
1293 if name
in assembly_comps:
1294 raise ValueError(
"Component %s already exists" % name)
1295 elif original
not in assembly_comps:
1296 raise ValueError(
"Original component %s does not exist" % original)
1297 self.create_component(state, name,
True)
1298 self.add_component_sequence(state, name, self.sequence_dict[original])
1299 self._transformed_components.append(_TransformedComponent(
1300 name, original, transform))
1301 self.all_representations.copy_component(state, name, original,
1302 self.asym_units[name])
1304 def create_component(self, state, name, modeled, asym_name=None):
1305 if asym_name
is None:
1307 new_comp = name
not in self._all_components
1308 self._all_components[name] =
None
1310 state.all_modeled_components.append(name)
1311 if asym_name
not in self.asym_units:
1313 self.asym_units[asym_name] =
None
1315 self.all_modeled_components.append(name)
1317 def add_component_sequence(self, state, name, seq, asym_name=None,
1319 if asym_name
is None:
1321 def get_offset(seq):
1323 for i
in range(len(seq)):
1326 seq, offset = get_offset(seq)
1327 if name
in self.sequence_dict:
1328 if self.sequence_dict[name] != seq:
1329 raise ValueError(
"Sequence mismatch for component %s" % name)
1331 self.sequence_dict[name] = seq
1332 self.entities.add(name, seq, offset, alphabet)
1333 if asym_name
in self.asym_units:
1334 if self.asym_units[asym_name]
is None:
1336 entity = self.entities[name]
1337 asym =
AsymUnit(entity, details=asym_name)
1338 self.system.asym_units.append(asym)
1339 self.asym_units[asym_name] = asym
1340 state.modeled_assembly.append(self.asym_units[asym_name])
1342 def _add_restraint_model_fits(self):
1343 """Add fits to restraints for all known models"""
1344 for group, m
in self.system._all_models():
1345 if m._is_restrained:
1346 for r
in self.system.restraints:
1347 if hasattr(r,
'add_fits_from_model_statfile'):
1348 r.add_fits_from_model_statfile(m)
1351 "Use finalize() then the python-ihm API instead")
1352 def flush(self, format='mmCIF'):
1353 """Write out all information to the file.
1354 Information can be written in any format supported by
1355 the ihm library (typically this is 'mmCIF' or 'BCIF')."""
1357 ihm.dumper.write(self.fh, [self.system], format, dumpers=
get_dumpers())
1360 """Do any final processing on the class hierarchy.
1361 After calling this method, the `system` member (an instance
1362 of `ihm.System`) completely reproduces the PMI modeling, and
1363 can be written out to an mmCIF file with `ihm.dumper.write`,
1364 and/or modified using the ihm API."""
1365 self._add_restraint_model_fits()
1367 def add_pdb_element(self, state, name, start, end, offset, pdbname,
1368 chain, hier, representation=
None):
1369 if self._is_excluded(name, start, end):
1371 if representation
is None:
1372 representation = self.default_representation
1373 asym = self.asym_units[name]
1374 p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1376 self.all_representations.add_fragment(state, representation, p)
1377 self.all_starting_models.add_pdb_fragment(p)
1379 def add_bead_element(self, state, name, start, end, num, hier,
1380 representation=
None):
1381 if self._is_excluded(name, start, end):
1383 if representation
is None:
1384 representation = self.default_representation
1385 asym = self.asym_units[name]
1386 pmi_offset = asym.entity.pmi_offset
1387 b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1389 self.all_representations.add_fragment(state, representation, b)
1391 def get_cross_link_group(self, pmi_restraint):
1392 r = _CrossLinkRestraint(pmi_restraint)
1393 self.system.restraints.append(r)
1394 self._add_restraint_dataset(r)
1397 def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1398 if c1
not in self._all_components
or c2
not in self._all_components:
1404 e1 = self.entities[c1]
1405 e2 = self.entities[c2]
1406 xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1407 residue2=e2.pmi_residue(r2))
1408 rsr.experimental_cross_links.append([xl])
1411 def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1414 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1415 d = ihm.restraint.UpperBoundDistanceRestraint(length)
1417 if _get_by_residue(p1)
and _get_by_residue(p2):
1418 cls = _ResidueCrossLink
1420 cls = _FeatureCrossLink
1421 xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1424 xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1425 rsr.cross_links.append(xl)
1427 def add_replica_exchange(self, state, rex):
1432 step = _ReplicaExchangeProtocolStep(state, rex)
1433 step.software = self.all_software.pmi
1434 self.all_protocols.add_step(step, state)
1436 def _add_simple_dynamics(self, num_models_end, method):
1438 state = self._last_state
1439 self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1442 def _add_protocol(self):
1444 state = self._last_state
1445 self.all_protocols.add_protocol(state)
1447 def _add_dataset(self, dataset):
1448 return self.all_datasets.add(self._last_state, dataset)
1450 def _add_restraint_dataset(self, restraint):
1451 return self.all_datasets.add_restraint(self._last_state, restraint)
1453 def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1455 state = self._last_state
1456 pp = ihm.analysis.ClusterStep(
'RMSD', num_models_begin, num_models_end)
1457 self.all_protocols.add_postproc(pp, state)
1460 def _add_no_postprocessing(self, num_models):
1462 state = self._last_state
1463 pp = ihm.analysis.EmptyStep()
1464 pp.num_models_begin = pp.num_models_end = num_models
1465 self.all_protocols.add_postproc(pp, state)
1468 def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1469 num_models_deposited, localization_densities,
1471 """Add an ensemble generated by ad hoc methods (not using PMI).
1472 This is currently only used by the Nup84 system."""
1474 state = self._last_state
1475 group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1476 state.add_model_group(group)
1478 self.system.locations.append(ensemble_file)
1479 e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1481 self.system.ensembles.append(e)
1482 for c
in state.all_modeled_components:
1483 den = localization_densities.get(c,
None)
1485 e.load_localization_density(state, c, self.asym_units[c], den)
1489 """Point a previously-created ensemble to an 'all-models' file.
1490 This could be a trajectory such as DCD, an RMF, or a multimodel
1492 self.system.locations.append(location)
1494 ind = i + self._state_ensemble_offset
1495 self.system.ensembles[ind].file = location
1497 def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1503 protocol = self.all_protocols.get_last_protocol(state)
1504 num_models = protocol.steps[-1].num_models_end
1505 pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1506 pp.software = self.all_software.pmi
1507 self.all_protocols.add_postproc(pp, state)
1508 for i
in range(rex._number_of_clusters):
1509 group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1510 'cluster %d' % (i + 1)))
1511 state.add_model_group(group)
1513 e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1514 self.system.ensembles.append(e)
1516 for fname, stuple
in sorted(density_custom_ranges.items()):
1517 e.load_localization_density(state, fname, stuple,
1519 for stats
in e.load_all_models(self, state):
1520 m = self.add_model(group)
1523 m.name =
'Best scoring model'
1526 for c
in state.all_modeled_components:
1529 def _get_subassembly(self, comps, name, description):
1530 """Get an Assembly consisting of the given components.
1531 `compdict` is a dictionary of the components to add, where keys
1532 are the component names and values are the sequence ranges (or
1533 None to use all residues in the component)."""
1535 for comp, seqrng
in comps.items():
1536 a = self.asym_units[comp]
1537 asyms.append(a
if seqrng
is None else a(*seqrng))
1539 a = ihm.Assembly(asyms, name=name, description=description)
1542 def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1544 """Add a basic FoXS fit. This is largely intended for use from the
1546 assembly = self._get_subassembly({comp:seqrange},
1547 name=
"SAXS subassembly",
1548 description=
"All components that fit SAXS data")
1549 r = ihm.restraint.SASRestraint(dataset, assembly, segment=
False,
1550 fitting_method=
'FoXS', fitting_atom_type=
'Heavy atoms',
1551 multi_state=
False, radius_of_gyration=rg, details=details)
1552 r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1553 self.system.restraints.append(r)
1554 self._add_restraint_dataset(r)
1556 def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1557 image_resolution, projection_number,
1558 micrographs_number):
1559 r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1560 image_resolution, projection_number,
1562 self.system.restraints.append(r)
1563 self._add_restraint_dataset(r)
1565 def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1567 r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1568 self.system.restraints.append(r)
1569 self._add_restraint_dataset(r)
1571 def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1572 sigma, pmi_restraint):
1573 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1574 sigma, pmi_restraint, self._xy_plane)
1576 def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1577 sigma, pmi_restraint):
1578 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1579 sigma, pmi_restraint, self._xz_plane)
1581 def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1582 sigma, pmi_restraint):
1583 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1584 sigma, pmi_restraint, self._z_axis)
1586 def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1587 sigma, pmi_restraint, geom):
1588 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1589 r = _GeometricRestraint(self, state, pmi_restraint, geom,
1590 asym.get_feature(ps),
1591 ihm.restraint.LowerUpperBoundDistanceRestraint(
1592 lower_bound, upper_bound),
1594 self.system.restraints.append(r)
1595 self._add_restraint_dataset(r)
1597 def _get_membrane(self, tor_R, tor_r, tor_th):
1598 """Get an object representing a half-torus membrane"""
1599 if not hasattr(self,
'_seen_membranes'):
1600 self._seen_membranes = {}
1603 membrane_id = tuple(int(x * 100.)
for x
in (tor_R, tor_r, tor_th))
1604 if membrane_id
not in self._seen_membranes:
1605 m = ihm.geometry.HalfTorus(center=self._center_origin,
1606 transformation=self._identity_transform,
1607 major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1608 inner=
True, name=
'Membrane')
1609 self._seen_membranes[membrane_id] = m
1610 return self._seen_membranes[membrane_id]
1612 def add_membrane_surface_location_restraint(
1613 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1614 self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1615 pmi_restraint, ihm.restraint.UpperBoundDistanceRestraint(0.))
1617 def add_membrane_exclusion_restraint(
1618 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1619 self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1620 pmi_restraint, ihm.restraint.LowerBoundDistanceRestraint(0.))
1622 def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1623 sigma, pmi_restraint, rsr):
1624 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1625 r = _GeometricRestraint(self, state, pmi_restraint,
1626 self._get_membrane(tor_R, tor_r, tor_th),
1627 asym.get_feature(ps), rsr, sigma)
1628 self.system.restraints.append(r)
1629 self._add_restraint_dataset(r)
1631 def add_model(self, group, assembly=None, representation=None):
1632 state = self._last_state
1633 if representation
is None:
1634 representation = self.default_representation
1635 protocol = self.all_protocols.get_last_protocol(state)
1636 m = _Model(state.prot, self, protocol,
1637 assembly
if assembly
else state.modeled_assembly,
1644 """Get custom python-ihm dumpers for writing PMI to from mmCIF.
1645 This returns a list of custom dumpers that can be passed as all or
1646 part of the `dumpers` argument to ihm.dumper.write(). They add
1647 PMI-specific information to mmCIF or BinaryCIF files written out
1649 return [_ReplicaExchangeProtocolDumper]
1653 """Get custom python-ihm handlers for reading PMI data from mmCIF.
1654 This returns a list of custom handlers that can be passed as all or
1655 part of the `handlers` argument to ihm.reader.read(). They read
1656 PMI-specific information from mmCIF or BinaryCIF files read in
1658 return [_ReplicaExchangeProtocolHandler]
1662 "Use ihm.reader.read() instead with handlers=get_handlers()")
1663 def read(fh, model_class=ihm.model.Model, format='mmCIF', handlers=[]):
1664 """Read data from the mmCIF file handle `fh`.
1665 This is a simple wrapper around `ihm.reader.read()`, which also
1666 reads PMI-specific information from the mmCIF or BinaryCIF file."""
1667 return ihm.reader.read(fh, model_class=model_class, format=format,
1672 """Extract metadata from an EM density GMM file."""
1675 """Extract metadata from `filename`.
1676 @return a dict with key `dataset` pointing to the GMM file and
1677 `number_of_gaussians` to the number of GMMs (or None)"""
1678 l = ihm.location.InputFileLocation(filename,
1679 details=
"Electron microscopy density map, "
1680 "represented as a Gaussian Mixture Model (GMM)")
1683 l._allow_duplicates =
True
1684 d = ihm.dataset.EMDensityDataset(l)
1685 ret = {
'dataset':d,
'number_of_gaussians':
None}
1687 with open(filename)
as fh:
1689 if line.startswith(
'# data_fn: '):
1690 p = ihm.metadata.MRCParser()
1691 fn = line[11:].rstrip(
'\r\n')
1692 dataset = p.parse_file(os.path.join(
1693 os.path.dirname(filename), fn))[
'dataset']
1694 ret[
'dataset'].parents.append(dataset)
1695 elif line.startswith(
'# ncenters: '):
1696 ret[
'number_of_gaussians'] = int(line[12:])
Select non water and non hydrogen atoms.
def get_handlers
Get custom python-ihm handlers for reading PMI data from mmCIF.
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.
void handle_use_deprecated(std::string message)
def deprecated_function
Python decorator to mark a function as deprecated.
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 deprecated_method
Python decorator to mark a method as deprecated.
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.
def get_dumpers
Get custom python-ihm dumpers for writing PMI to from mmCIF.
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.
Mapping between FASTA one-letter codes and residue types.
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 ...