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
39 import ihm.representation
41 import ihm.cross_linkers
44 def _assign_id(obj, seen_objs, obj_by_id):
45 """Assign a unique ID to obj, and track all ids in obj_by_id."""
46 if obj
not in seen_objs:
47 if not hasattr(obj,
'id'):
49 obj.id = len(obj_by_id)
50 seen_objs[obj] = obj.id
52 obj.id = seen_objs[obj]
55 def _get_by_residue(p):
56 """Determine whether the given particle represents a specific residue
57 or a more coarse-grained object."""
61 class _ComponentMapper(object):
62 """Map a Particle to a component name"""
63 def __init__(self, prot):
66 self.name =
'cif-output'
67 self.o.dictionary_pdbs[self.name] = self.prot
68 self.o._init_dictchain(self.name, self.prot, multichar_chain=
True)
70 def __getitem__(self, p):
71 protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
75 class _AsymMapper(object):
76 """Map a Particle to an asym_unit"""
77 def __init__(self, simo, prot):
79 self._cm = _ComponentMapper(prot)
80 self._seen_ranges = {}
82 def __getitem__(self, p):
83 protname = self._cm[p]
84 return self.simo.asym_units[protname]
86 def get_feature(self, ps):
87 """Get an ihm.restraint.Feature that covers the given particles"""
95 rng = asym(rind, rind)
99 rng = asym(rinds[0], rinds[-1])
101 raise ValueError(
"Unsupported particle type %s" % str(p))
103 if len(rngs) > 0
and rngs[-1].asym == asym \
104 and rngs[-1].seq_id_range[1] == rng.seq_id_range[0] - 1:
105 rngs[-1].seq_id_range = (rngs[-1].seq_id_range[0],
112 if hrngs
in self._seen_ranges:
113 return self._seen_ranges[hrngs]
115 feat = ihm.restraint.ResidueFeature(rngs)
116 self._seen_ranges[hrngs] = feat
120 class _AllSoftware(object):
121 def __init__(self, system):
123 self.modeller_used = self.phyre2_used =
False
124 self.pmi = ihm.Software(
125 name=
"IMP PMI module",
126 version=IMP.pmi.__version__,
127 classification=
"integrative model building",
128 description=
"integrative model building",
129 location=
'https://integrativemodeling.org')
130 self.imp = ihm.Software(
131 name=
"Integrative Modeling Platform (IMP)",
132 version=IMP.__version__,
133 classification=
"integrative model building",
134 description=
"integrative model building",
135 location=
'https://integrativemodeling.org')
136 self.system.software.extend([self.pmi, self.imp])
138 def set_modeller_used(self, version, date):
139 if self.modeller_used:
141 self.modeller_used =
True
142 self.system.software.append(ihm.Software(
143 name=
'MODELLER', classification=
'comparative modeling',
144 description=
'Comparative modeling by satisfaction '
145 'of spatial restraints, build ' + date,
146 location=
'https://salilab.org/modeller/',
149 def set_phyre2_used(self):
152 self.phyre2_used =
True
153 self.system.software.append(ihm.Software(
154 name=
'Phyre2', classification=
'protein homology modeling',
155 description=
'Protein Homology/analogY Recognition '
158 location=
'http://www.sbg.bio.ic.ac.uk/~phyre2/'))
161 def _get_fragment_is_rigid(fragment):
162 """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."""
195 def __init__(self, state, component, start, end, count, hier, asym_unit):
196 super(_BeadsFragment, self).__init__(
197 asym_unit=asym_unit(start, end), rigid=
None, primitive=
'sphere',
199 self.state, self.component, self.hier = state, component, hier
201 rigid = property(
lambda self: _get_fragment_is_rigid(self),
202 lambda self, val:
None)
204 def combine(self, other):
206 if (type(other) == type(self)
and
207 other.asym_unit.seq_id_range[0]
208 == self.asym_unit.seq_id_range[1] + 1):
209 self.asym_unit.seq_id_range = (self.asym_unit.seq_id_range[0],
210 other.asym_unit.seq_id_range[1])
211 self.count += other.count
215 class _AllModelRepresentations(object):
216 def __init__(self, simo):
220 self.fragments = OrderedDict()
221 self._all_representations = {}
223 def copy_component(self, state, name, original, asym_unit):
224 """Copy all representation for `original` in `state` to `name`"""
227 newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
229 for rep
in self.fragments:
230 if original
in self.fragments[rep]:
231 if name
not in self.fragments[rep]:
232 self.fragments[rep][name] = OrderedDict()
233 self.fragments[rep][name][state] = [
234 copy_frag(f)
for f
in self.fragments[rep][original][state]]
237 first_state = list(self.fragments[rep][name].keys())[0]
238 if state
is first_state:
239 representation = self._all_representations[rep]
240 representation.extend(self.fragments[rep][name][state])
242 def add_fragment(self, state, representation, fragment):
243 """Add a model fragment."""
244 comp = fragment.component
245 id_rep = id(representation)
246 self._all_representations[id_rep] = representation
247 if id_rep
not in self.fragments:
248 self.fragments[id_rep] = OrderedDict()
249 if comp
not in self.fragments[id_rep]:
250 self.fragments[id_rep][comp] = OrderedDict()
251 if state
not in self.fragments[id_rep][comp]:
252 self.fragments[id_rep][comp][state] = []
253 fragments = self.fragments[id_rep][comp][state]
254 if len(fragments) == 0
or not fragments[-1].combine(fragment):
255 fragments.append(fragment)
258 first_state = list(self.fragments[id_rep][comp].keys())[0]
259 if state
is first_state:
260 representation.append(fragment)
263 class _AllDatasets(object):
264 """Track all datasets generated by PMI and add them to the ihm.System"""
265 def __init__(self, system):
268 self._datasets_by_state = {}
269 self._restraints_by_state = {}
271 def get_all_group(self, state):
272 """Get a DatasetGroup encompassing all datasets so far in this state"""
276 g = ihm.dataset.DatasetGroup(
277 self._datasets_by_state.get(state, [])
278 + [r.dataset
for r
in self._restraints_by_state.get(state, [])
282 def add(self, state, dataset):
283 """Add a new dataset."""
284 self._datasets.append(dataset)
285 if state
not in self._datasets_by_state:
286 self._datasets_by_state[state] = []
287 self._datasets_by_state[state].append(dataset)
290 self.system.orphan_datasets.append(dataset)
294 """Add the dataset for a restraint"""
295 if state
not in self._restraints_by_state:
296 self._restraints_by_state[state] = []
297 self._restraints_by_state[state].append(restraint)
300 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
301 """Restrain to a set of cross-links"""
304 _label_map = {
'wtDSS':
'DSS',
'scDSS':
'DSS',
'scEDC':
'EDC'}
305 _descriptor_map = {
'DSS': ihm.cross_linkers.dss,
306 'EDC': ihm.cross_linkers.edc}
308 def __init__(self, pmi_restraint):
309 self.pmi_restraint = pmi_restraint
312 linker = getattr(self.pmi_restraint,
'linker',
None)
313 label = self.pmi_restraint.label
315 super(_CrossLinkRestraint, self).__init__(
316 dataset=self.pmi_restraint.dataset,
317 linker=linker
or self._get_chem_descriptor(label))
320 def _get_chem_descriptor(cls, label):
322 label = cls._label_map.get(label, label)
323 if label
not in cls._descriptor_map:
327 d = ihm.ChemDescriptor(label)
328 cls._descriptor_map[label] = d
329 return cls._descriptor_map[label]
331 def _set_psi_sigma(self, model):
334 if model.m != self.pmi_restraint.model:
336 for resolution
in self.pmi_restraint.sigma_dictionary:
337 statname =
'ISDCrossLinkMS_Sigma_%s_%s' % (resolution, self.label)
338 if model.stats
and statname
in model.stats:
339 sigma = float(model.stats[statname])
340 p = self.pmi_restraint.sigma_dictionary[resolution][0]
341 old_values.append((p, p.get_scale()))
343 for psiindex
in self.pmi_restraint.psi_dictionary:
344 statname =
'ISDCrossLinkMS_Psi_%s_%s' % (psiindex, self.label)
345 if model.stats
and statname
in model.stats:
346 psi = float(model.stats[statname])
347 p = self.pmi_restraint.psi_dictionary[psiindex][0]
348 old_values.append((p, p.get_scale()))
352 return list(reversed(old_values))
354 def add_fits_from_model_statfile(self, model):
356 old_values = self._set_psi_sigma(model)
360 for xl
in self.cross_links:
362 xl.fits[model] = ihm.restraint.CrossLinkFit(
363 psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
365 for p, scale
in old_values:
369 def __set_dataset(self, val):
370 self.pmi_restraint.dataset = val
371 dataset = property(
lambda self: self.pmi_restraint.dataset,
375 def get_asym_mapper_for_state(simo, state, asym_map):
376 asym = asym_map.get(state,
None)
378 asym = _AsymMapper(simo, state.prot)
379 asym_map[state] = asym
383 class _PMICrossLink(object):
386 psi = property(
lambda self: self.psi_p.get_scale(),
387 lambda self, val:
None)
388 sigma1 = property(
lambda self: self.sigma1_p.get_scale(),
389 lambda self, val:
None)
390 sigma2 = property(
lambda self: self.sigma2_p.get_scale(),
391 lambda self, val:
None)
394 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
398 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
402 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
403 def __init__(self, state, pmi_restraint, image_number, resolution,
404 pixel_size, image_resolution, projection_number,
406 self.pmi_restraint, self.image_number = pmi_restraint, image_number
407 super(_EM2DRestraint, self).__init__(
408 dataset=pmi_restraint.datasets[image_number],
409 assembly=state.modeled_assembly,
410 segment=
False, number_raw_micrographs=micrographs_number,
411 pixel_size_width=pixel_size, pixel_size_height=pixel_size,
412 image_resolution=image_resolution,
413 number_of_projections=projection_number)
416 def __get_dataset(self):
417 return self.pmi_restraint.datasets[self.image_number]
419 def __set_dataset(self, val):
420 self.pmi_restraint.datasets[self.image_number] = val
422 dataset = property(__get_dataset, __set_dataset)
424 def add_fits_from_model_statfile(self, model):
425 ccc = self._get_cross_correlation(model)
426 transform = self._get_transformation(model)
427 rot = transform.get_rotation()
428 rm = [[e
for e
in rot.get_rotation_matrix_row(i)]
for i
in range(3)]
429 self.fits[model] = ihm.restraint.EM2DRestraintFit(
430 cross_correlation_coefficient=ccc,
432 tr_vector=transform.get_translation())
434 def _get_transformation(self, model):
435 """Get the transformation that places the model on the image"""
436 stats = model.em2d_stats
or model.stats
437 prefix =
'ElectronMicroscopy2D_%s_Image%d' % (self.pmi_restraint.label,
438 self.image_number + 1)
439 r = [float(stats[prefix +
'_Rotation%d' % i])
for i
in range(4)]
440 t = [float(stats[prefix +
'_Translation%d' % i])
444 inv = model.transform.get_inverse()
448 def _get_cross_correlation(self, model):
449 """Get the cross correlation coefficient between the model projection
451 stats = model.em2d_stats
or model.stats
452 return float(stats[
'ElectronMicroscopy2D_%s_Image%d_CCC'
453 % (self.pmi_restraint.label,
454 self.image_number + 1)])
457 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
459 def __init__(self, simo, state, pmi_restraint, target_ps, densities):
460 self.pmi_restraint = pmi_restraint
461 super(_EM3DRestraint, self).__init__(
462 dataset=pmi_restraint.dataset,
463 assembly=self._get_assembly(densities, simo, state),
464 fitting_method=
'Gaussian mixture models',
465 number_of_gaussians=len(target_ps))
468 def __set_dataset(self, val):
469 self.pmi_restraint.dataset = val
470 dataset = property(
lambda self: self.pmi_restraint.dataset,
473 def _get_assembly(self, densities, simo, state):
474 """Get the Assembly that this restraint acts on"""
475 cm = _ComponentMapper(state.prot)
478 components[cm[d]] =
None
479 a = simo._get_subassembly(
480 components, name=
"EM subassembly",
481 description=
"All components that fit the EM map")
484 def add_fits_from_model_statfile(self, model):
485 ccc = self._get_cross_correlation(model)
486 self.fits[model] = ihm.restraint.EM3DRestraintFit(
487 cross_correlation_coefficient=ccc)
489 def _get_cross_correlation(self, model):
490 """Get the cross correlation coefficient between the model
492 if model.stats
is not None:
493 return float(model.stats[
'GaussianEMRestraint_%s_CCC'
494 % self.pmi_restraint.label])
497 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
499 def __init__(self, simo, state, pmi_restraint, geometric_object,
500 feature, distance, sigma):
501 self.pmi_restraint = pmi_restraint
502 super(_GeometricRestraint, self).__init__(
503 dataset=pmi_restraint.dataset,
504 geometric_object=geometric_object, feature=feature,
505 distance=distance, harmonic_force_constant=1. / sigma,
509 def __set_dataset(self, val):
510 self.pmi_restraint.dataset = val
511 dataset = property(
lambda self: self.pmi_restraint.dataset,
515 class _ReplicaExchangeProtocolStep(ihm.protocol.Step):
516 def __init__(self, state, rex):
517 if rex.monte_carlo_sample_objects
is not None:
518 method =
'Replica exchange monte carlo'
520 method =
'Replica exchange molecular dynamics'
521 self.monte_carlo_temperature = rex.vars[
'monte_carlo_temperature']
522 self.replica_exchange_minimum_temperature = \
523 rex.vars[
'replica_exchange_minimum_temperature']
524 self.replica_exchange_maximum_temperature = \
525 rex.vars[
'replica_exchange_maximum_temperature']
526 super(_ReplicaExchangeProtocolStep, self).__init__(
527 assembly=state.modeled_assembly,
529 method=method, name=
'Sampling',
530 num_models_begin=
None,
531 num_models_end=rex.vars[
"number_of_frames"],
532 multi_scale=
True, multi_state=
False, ordered=
False)
535 class _ReplicaExchangeProtocolDumper(ihm.dumper.Dumper):
536 """Write IMP-specific information about replica exchange to mmCIF.
537 Note that IDs will have already been assigned by python-ihm's
538 standard modeling protocol dumper."""
539 def dump(self, system, writer):
540 with writer.loop(
"_imp_replica_exchange_protocol",
541 [
"protocol_id",
"step_id",
"monte_carlo_temperature",
542 "replica_exchange_minimum_temperature",
543 "replica_exchange_maximum_temperature"])
as lp:
544 for p
in system._all_protocols():
546 if isinstance(s, _ReplicaExchangeProtocolStep):
547 self._dump_step(p, s, lp)
549 def _dump_step(self, p, s, lp):
550 mintemp = s.replica_exchange_minimum_temperature
551 maxtemp = s.replica_exchange_maximum_temperature
552 lp.write(protocol_id=p._id, step_id=s._id,
553 monte_carlo_temperature=s.monte_carlo_temperature,
554 replica_exchange_minimum_temperature=mintemp,
555 replica_exchange_maximum_temperature=maxtemp)
558 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
559 category =
'_imp_replica_exchange_protocol'
561 """Read IMP-specific information about replica exchange from mmCIF."""
562 def __call__(self, protocol_id, step_id, monte_carlo_temperature,
563 replica_exchange_minimum_temperature,
564 replica_exchange_maximum_temperature):
565 p = self.sysr.protocols.get_by_id(protocol_id)
567 s = p.steps[int(step_id)-1]
569 s.__class__ = _ReplicaExchangeProtocolStep
570 s.monte_carlo_temperature = \
571 self.get_float(monte_carlo_temperature)
572 s.replica_exchange_minimum_temperature = \
573 self.get_float(replica_exchange_minimum_temperature)
574 s.replica_exchange_maximum_temperature = \
575 self.get_float(replica_exchange_maximum_temperature)
578 class _SimpleProtocolStep(ihm.protocol.Step):
579 def __init__(self, state, num_models_end, method):
580 super(_SimpleProtocolStep, self).__init__(
581 assembly=state.modeled_assembly,
583 method=method, name=
'Sampling',
584 num_models_begin=
None,
585 num_models_end=num_models_end,
586 multi_scale=
True, multi_state=
False, ordered=
False)
589 class _Chain(object):
590 """Represent a single chain in a Model"""
591 def __init__(self, pmi_chain_id, asym_unit):
592 self.pmi_chain_id, self.asym_unit = pmi_chain_id, asym_unit
596 def add(self, xyz, atom_type, residue_type, residue_index,
597 all_indexes, radius):
598 if atom_type
is None:
599 self.spheres.append((xyz, residue_type, residue_index,
600 all_indexes, radius))
602 self.atoms.append((xyz, atom_type, residue_type, residue_index,
603 all_indexes, radius))
604 orig_comp = property(
lambda self: self.comp)
607 class _TransformedChain(object):
608 """Represent a chain that is a transformed version of another"""
609 def __init__(self, orig_chain, asym_unit, transform):
610 self.orig_chain, self.asym_unit = orig_chain, asym_unit
611 self.transform = transform
613 def __get_spheres(self):
614 for (xyz, residue_type, residue_index, all_indexes,
615 radius)
in self.orig_chain.spheres:
616 yield (self.transform * xyz, residue_type, residue_index,
618 spheres = property(__get_spheres)
620 def __get_atoms(self):
621 for (xyz, atom_type, residue_type, residue_index, all_indexes,
622 radius)
in self.orig_chain.atoms:
623 yield (self.transform * xyz, atom_type, residue_type,
624 residue_index, all_indexes, radius)
625 atoms = property(__get_atoms)
627 entity = property(
lambda self: self.orig_chain.entity)
628 orig_comp = property(
lambda self: self.orig_chain.comp)
631 class _Excluder(object):
632 def __init__(self, component, simo):
633 self._seqranges = simo._exclude_coords.get(component, [])
635 def is_excluded(self, indexes):
636 """Return True iff the given sequence range is excluded."""
637 for seqrange
in self._seqranges:
638 if indexes[0] >= seqrange[0]
and indexes[-1] <= seqrange[1]:
642 class _Model(ihm.model.Model):
643 def __init__(self, prot, simo, protocol, assembly, representation):
644 super(_Model, self).__init__(assembly=assembly, protocol=protocol,
645 representation=representation)
646 self.simo = weakref.proxy(simo)
652 self.em2d_stats =
None
655 self._is_restrained =
True
658 self.m = prot.get_model()
659 o.dictionary_pdbs[name] = prot
660 o._init_dictchain(name, prot, multichar_chain=
True)
661 (particle_infos_for_pdb,
662 self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
664 self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
667 def all_chains(self, simo):
668 """Yield all chains, including transformed ones"""
670 for c
in self.chains:
672 chain_for_comp[c.comp] = c
673 for tc
in simo._transformed_components:
674 orig_chain = chain_for_comp.get(tc.original,
None)
676 asym = simo.asym_units[tc.name]
677 c = _TransformedChain(orig_chain, asym, tc.transform)
681 def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
682 entity_for_chain = {}
685 for protname, chain_id
in o.dictchain[name].items():
686 if protname
in simo.entities:
687 entity_for_chain[chain_id] = simo.entities[protname]
690 pn = protname.split(
'.')[0]
691 entity_for_chain[chain_id] = simo.entities[pn]
692 comp_for_chain[chain_id] = protname
696 correct_asym[chain_id] = simo.asym_units[protname]
703 for (xyz, atom_type, residue_type, chain_id, residue_index,
704 all_indexes, radius)
in particle_infos_for_pdb:
705 if chain
is None or chain.pmi_chain_id != chain_id:
706 chain = _Chain(chain_id, correct_asym[chain_id])
707 chain.entity = entity_for_chain[chain_id]
708 chain.comp = comp_for_chain[chain_id]
709 self.chains.append(chain)
710 excluder = _Excluder(chain.comp, simo)
711 if not excluder.is_excluded(all_indexes
if all_indexes
712 else [residue_index]):
713 chain.add(xyz, atom_type, residue_type, residue_index,
716 def parse_rmsf_file(self, fname, component):
717 self.rmsf[component] = rmsf = {}
718 with open(fname)
as fh:
720 resnum, blocknum, val = line.split()
721 rmsf[int(resnum)] = (int(blocknum), float(val))
723 def get_rmsf(self, component, indexes):
724 """Get the RMSF value for the given residue indexes."""
727 rmsf = self.rmsf[component]
728 blocknums = dict.fromkeys(rmsf[ind][0]
for ind
in indexes)
729 if len(blocknums) != 1:
730 raise ValueError(
"Residue indexes %s aren't all in the same block"
732 return rmsf[indexes[0]][1]
735 for chain
in self.all_chains(self.simo):
736 pmi_offset = chain.asym_unit.entity.pmi_offset
737 for atom
in chain.atoms:
738 (xyz, atom_type, residue_type, residue_index,
739 all_indexes, radius) = atom
740 pt = self.transform * xyz
741 yield ihm.model.Atom(
742 asym_unit=chain.asym_unit,
743 seq_id=residue_index - pmi_offset,
744 atom_id=atom_type.get_string(),
746 x=pt[0], y=pt[1], z=pt[2])
748 def get_spheres(self):
749 for chain
in self.all_chains(self.simo):
750 pmi_offset = chain.asym_unit.entity.pmi_offset
751 for sphere
in chain.spheres:
752 (xyz, residue_type, residue_index,
753 all_indexes, radius) = sphere
754 if all_indexes
is None:
755 all_indexes = (residue_index,)
756 pt = self.transform * xyz
757 yield ihm.model.Sphere(
758 asym_unit=chain.asym_unit,
759 seq_id_range=(all_indexes[0] - pmi_offset,
760 all_indexes[-1] - pmi_offset),
761 x=pt[0], y=pt[1], z=pt[2], radius=radius,
762 rmsf=self.get_rmsf(chain.orig_comp, all_indexes))
765 class _AllProtocols(object):
766 def __init__(self, simo):
767 self.simo = weakref.proxy(simo)
769 self.protocols = OrderedDict()
771 def add_protocol(self, state):
772 """Add a new Protocol"""
773 if state
not in self.protocols:
774 self.protocols[state] = []
775 p = ihm.protocol.Protocol()
776 self.simo.system.orphan_protocols.append(p)
777 self.protocols[state].append(p)
779 def add_step(self, step, state):
780 """Add a ProtocolStep to the last Protocol of the given State"""
781 if state
not in self.protocols:
782 self.add_protocol(state)
783 protocol = self.get_last_protocol(state)
784 if len(protocol.steps) == 0:
785 step.num_models_begin = 0
787 step.num_models_begin = protocol.steps[-1].num_models_end
788 protocol.steps.append(step)
789 step.id = len(protocol.steps)
791 step.dataset_group = self.simo.all_datasets.get_all_group(state)
793 def add_postproc(self, step, state):
794 """Add a postprocessing step to the last protocol"""
795 protocol = self.get_last_protocol(state)
796 if len(protocol.analyses) == 0:
797 protocol.analyses.append(ihm.analysis.Analysis())
798 protocol.analyses[-1].steps.append(step)
800 def get_last_protocol(self, state):
801 """Return the most recently-added _Protocol"""
802 return self.protocols[state][-1]
805 class _AllStartingModels(object):
806 def __init__(self, simo):
810 self.models = OrderedDict()
813 def add_pdb_fragment(self, fragment):
814 """Add a starting model PDB fragment."""
815 comp = fragment.component
816 state = fragment.state
817 if comp
not in self.models:
818 self.models[comp] = OrderedDict()
819 if state
not in self.models[comp]:
820 self.models[comp][state] = []
821 models = self.models[comp][state]
822 if len(models) == 0 \
823 or models[-1].fragments[0].pdbname != fragment.pdbname:
824 model = self._add_model(fragment)
828 models[-1].fragments.append(weakref.proxy(fragment))
830 pmi_offset = models[-1].asym_unit.entity.pmi_offset
831 sid_begin = min(fragment.start + fragment.offset - pmi_offset,
832 models[-1].asym_unit.seq_id_range[0])
833 sid_end = max(fragment.end + fragment.offset - pmi_offset,
834 models[-1].asym_unit.seq_id_range[1])
835 models[-1].asym_unit = fragment.asym_unit.asym(sid_begin, sid_end)
836 fragment.starting_model = models[-1]
838 def _add_model(self, f):
839 parser = ihm.metadata.PDBParser()
840 r = parser.parse_file(f.pdbname)
842 self.simo._add_dataset(r[
'dataset'])
844 templates = r[
'templates'].get(f.chain, [])
847 self.simo.system.locations.append(t.alignment_file)
849 self.simo._add_dataset(t.dataset)
850 source = r[
'entity_source'].get(f.chain)
852 f.asym_unit.entity.source = source
853 pmi_offset = f.asym_unit.entity.pmi_offset
855 asym_unit=f.asym_unit.asym.pmi_range(f.start, f.end),
856 dataset=r[
'dataset'], asym_id=f.chain,
857 templates=templates, offset=f.offset + pmi_offset,
858 metadata=r[
'metadata'],
859 software=r[
'software'][0]
if r[
'software']
else None,
860 script_file=r[
'script'])
861 m.fragments = [weakref.proxy(f)]
865 class _StartingModel(ihm.startmodel.StartingModel):
866 def get_seq_dif(self):
870 pmi_offset = self.asym_unit.entity.pmi_offset
871 mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
873 for f
in self.fragments:
876 residue_indexes=list(range(f.start - f.offset,
877 f.end - f.offset + 1)))
878 for a
in mh.get_ihm_atoms(sel.get_selected_particles(),
879 f.offset - pmi_offset):
881 self._seq_dif = mh._seq_dif
884 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
885 """Post processing using AnalysisReplicaExchange0 macro"""
887 def __init__(self, rex, num_models_begin):
890 for fname
in self.get_all_stat_files():
891 with open(fname)
as fh:
892 num_models_end += len(fh.readlines())
893 super(_ReplicaExchangeAnalysisPostProcess, self).__init__(
894 feature=
'RMSD', num_models_begin=num_models_begin,
895 num_models_end=num_models_end)
897 def get_stat_file(self, cluster_num):
898 return os.path.join(self.rex._outputdir,
"cluster.%d" % cluster_num,
901 def get_all_stat_files(self):
902 for i
in range(self.rex._number_of_clusters):
903 yield self.get_stat_file(i)
906 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
907 """Ensemble generated using AnalysisReplicaExchange0 macro"""
909 num_models_deposited =
None
911 def __init__(self, pp, cluster_num, model_group, num_deposit):
912 with open(pp.get_stat_file(cluster_num))
as fh:
913 num_models = len(fh.readlines())
914 super(_ReplicaExchangeAnalysisEnsemble, self).__init__(
915 num_models=num_models,
916 model_group=model_group, post_process=pp,
917 clustering_feature=pp.feature,
918 name=model_group.name)
919 self.cluster_num = cluster_num
920 self.num_models_deposited = num_deposit
922 def get_rmsf_file(self, component):
923 return os.path.join(self.post_process.rex._outputdir,
924 'cluster.%d' % self.cluster_num,
925 'rmsf.%s.dat' % component)
927 def load_rmsf(self, model, component):
928 fname = self.get_rmsf_file(component)
929 if os.path.exists(fname):
930 model.parse_rmsf_file(fname, component)
932 def get_localization_density_file(self, fname):
933 return os.path.join(self.post_process.rex._outputdir,
934 'cluster.%d' % self.cluster_num,
937 def load_localization_density(self, state, fname, select_tuple,
939 fullpath = self.get_localization_density_file(fname)
940 if os.path.exists(fullpath):
941 details =
"Localization density for %s %s" \
942 % (fname, self.model_group.name)
943 local_file = ihm.location.OutputFileLocation(fullpath,
945 for s
in select_tuple:
946 if isinstance(s, tuple)
and len(s) == 3:
947 asym = asym_units[s[2]].pmi_range(s[0], s[1])
950 den = ihm.model.LocalizationDensity(file=local_file,
952 self.densities.append(den)
954 def load_all_models(self, simo, state):
955 stat_fname = self.post_process.get_stat_file(self.cluster_num)
957 with open(stat_fname)
as fh:
958 stats = ast.literal_eval(fh.readline())
960 rmf_file = os.path.join(os.path.dirname(stat_fname),
961 "%d.rmf3" % model_num)
963 if os.path.exists(rmf_file):
964 rh = RMF.open_rmf_file_read_only(rmf_file)
965 system = state._pmi_object.system
971 if model_num >= self.num_models_deposited:
975 def _get_precision(self):
976 precfile = os.path.join(self.post_process.rex._outputdir,
977 "precision.%d.%d.out" % (self.cluster_num,
979 if not os.path.exists(precfile):
983 r'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
985 with open(precfile)
as fh:
989 return float(m.group(1))
991 precision = property(
lambda self: self._get_precision(),
992 lambda self, val:
None)
995 class _SimpleEnsemble(ihm.model.Ensemble):
996 """Simple manually-created ensemble"""
998 num_models_deposited =
None
1000 def __init__(self, pp, model_group, num_models, drmsd,
1001 num_models_deposited, ensemble_file):
1002 super(_SimpleEnsemble, self).__init__(
1003 model_group=model_group, post_process=pp, num_models=num_models,
1004 file=ensemble_file, precision=drmsd, name=model_group.name,
1005 clustering_feature=
'dRMSD')
1006 self.num_models_deposited = num_models_deposited
1008 def load_localization_density(self, state, component, asym, local_file):
1009 den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
1010 self.densities.append(den)
1013 class _CustomDNAAlphabet(object):
1014 """Custom DNA alphabet that maps A,C,G,T (rather than DA,DC,DG,DT
1015 as in python-ihm)"""
1016 _comps = dict([cc.code_canonical, cc]
1017 for cc
in ihm.DNAAlphabet._comps.values())
1020 class _EntityMapper(dict):
1021 """Handle mapping from IMP components (without copy number) to CIF
1022 entities. Multiple components may map to the same entity if they
1024 def __init__(self, system):
1025 super(_EntityMapper, self).__init__()
1026 self._sequence_dict = {}
1028 self.system = system
1030 def _get_alphabet(self, alphabet):
1031 """Map a PMI alphabet to an IHM alphabet"""
1035 alphabet_map = {
None: ihm.LPeptideAlphabet,
1036 IMP.pmi.alphabets.amino_acid: ihm.LPeptideAlphabet,
1037 IMP.pmi.alphabets.rna: ihm.RNAAlphabet,
1038 IMP.pmi.alphabets.dna: _CustomDNAAlphabet}
1039 if alphabet
in alphabet_map:
1040 return alphabet_map[alphabet]
1042 raise TypeError(
"Don't know how to handle %s" % alphabet)
1044 def add(self, component_name, sequence, offset, alphabet):
1045 def entity_seq(sequence):
1048 return [
'UNK' if s ==
'X' else s
for s
in sequence]
1051 if sequence
not in self._sequence_dict:
1054 d = component_name.split(
"@")[0].split(
".")[0]
1055 entity = Entity(entity_seq(sequence), description=d,
1057 alphabet=self._get_alphabet(alphabet))
1058 self.system.entities.append(entity)
1059 self._sequence_dict[sequence] = entity
1060 self[component_name] = self._sequence_dict[sequence]
1063 class _TransformedComponent(object):
1064 def __init__(self, name, original, transform):
1065 self.name, self.original, self.transform = name, original, transform
1068 class _SimpleRef(object):
1069 """Class with similar interface to weakref.ref, but keeps a strong ref"""
1070 def __init__(self, ref):
1077 class _State(ihm.model.State):
1078 """Representation of a single state in the system."""
1080 def __init__(self, pmi_object, po):
1085 self._pmi_object = weakref.proxy(pmi_object)
1086 if hasattr(pmi_object,
'state'):
1089 self._pmi_state = _SimpleRef(pmi_object.state)
1091 self._pmi_state = weakref.ref(pmi_object)
1093 old_name = self.name
1094 super(_State, self).__init__(experiment_type=
'Fraction of bulk')
1095 self.name = old_name
1099 self.modeled_assembly = ihm.Assembly(
1100 name=
"Modeled assembly",
1101 description=self.get_postfixed_name(
1102 "All components modeled by IMP"))
1103 po.system.orphan_assemblies.append(self.modeled_assembly)
1105 self.all_modeled_components = []
1108 return hash(self._pmi_state())
1110 def __eq__(self, other):
1111 return self._pmi_state() == other._pmi_state()
1113 def add_model_group(self, group):
1117 def get_prefixed_name(self, name):
1118 """Prefix the given name with the state name, if available."""
1120 return self.short_name +
' ' + name
1124 return name[0].upper() + name[1:]
if name
else ''
1126 def get_postfixed_name(self, name):
1127 """Postfix the given name with the state name, if available."""
1129 return "%s in state %s" % (name, self.short_name)
1133 short_name = property(
lambda self: self._pmi_state().short_name)
1134 long_name = property(
lambda self: self._pmi_state().long_name)
1136 def __get_name(self):
1137 return self._pmi_state().long_name
1139 def __set_name(self, val):
1140 self._pmi_state().long_name = val
1142 name = property(__get_name, __set_name)
1146 """A single entity in the system.
1147 This functions identically to the base ihm.Entity class, but it
1148 allows identifying residues by either the IHM numbering scheme
1149 (seq_id, which is always contiguous starting at 1) or by the PMI
1150 scheme (which is similar, but may not start at 1 if the sequence in
1151 the FASTA file has one or more N-terminal gaps"""
1152 def __init__(self, sequence, pmi_offset, *args, **keys):
1155 self.pmi_offset = pmi_offset
1156 super(Entity, self).__init__(sequence, *args, **keys)
1159 """Return a single IHM residue indexed using PMI numbering"""
1160 return self.residue(res_id - self.pmi_offset)
1163 """Return a range of IHM residues indexed using PMI numbering"""
1164 off = self.pmi_offset
1165 return self(res_id_begin - off, res_id_end - off)
1169 """A single asymmetric unit in the system.
1170 This functions identically to the base ihm.AsymUnit class, but it
1171 allows identifying residues by either the IHM numbering scheme
1172 (seq_id, which is always contiguous starting at 1) or by the PMI
1173 scheme (which is similar, but may not start at 1 if the sequence in
1174 the FASTA file has one or more N-terminal gaps"""
1176 def __init__(self, entity, *args, **keys):
1177 super(AsymUnit, self).__init__(
1178 entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1181 """Return a single IHM residue indexed using PMI numbering"""
1182 return self.residue(res_id - self.entity.pmi_offset)
1185 """Return a range of IHM residues indexed using PMI numbering"""
1186 off = self.entity.pmi_offset
1187 return self(res_id_begin - off, res_id_end - off)
1191 """Class to encode a modeling protocol as mmCIF.
1193 IMP has basic support for writing out files in mmCIF format, for
1194 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1195 After creating an instance of this class, attach it to an
1196 IMP.pmi.topology.System object. After this, any
1197 generated models and metadata are automatically collected in the
1198 `system` attribute, which is an
1199 [ihm.System](https://python-ihm.readthedocs.io/en/latest/main.html#ihm.System) object.
1200 Once the protocol is complete, call finalize() to make sure `system`
1201 contains everything, then use the
1202 [python-ihm API](https://python-ihm.readthedocs.io/en/latest/dumper.html#ihm.dumper.write)
1203 to write out files in mmCIF or BinaryCIF format.
1205 See also get_handlers(), get_dumpers().
1209 self.system = ihm.System()
1210 self._state_group = ihm.model.StateGroup()
1211 self.system.state_groups.append(self._state_group)
1213 self._state_ensemble_offset = 0
1214 self._main_script = os.path.abspath(sys.argv[0])
1217 loc = ihm.location.WorkflowFileLocation(
1218 path=self._main_script,
1219 details=
"The main integrative modeling script")
1220 self.system.locations.append(loc)
1223 self.__asym_states = {}
1224 self._working_directory = os.getcwd()
1226 "Default representation")
1227 self.entities = _EntityMapper(self.system)
1229 self.asym_units = {}
1230 self._all_components = {}
1231 self.all_modeled_components = []
1232 self._transformed_components = []
1233 self.sequence_dict = {}
1236 self._xy_plane = ihm.geometry.XYPlane()
1237 self._xz_plane = ihm.geometry.XZPlane()
1238 self._z_axis = ihm.geometry.ZAxis()
1239 self._center_origin = ihm.geometry.Center(0, 0, 0)
1240 self._identity_transform = ihm.geometry.Transformation.identity()
1243 self._exclude_coords = {}
1245 self.all_representations = _AllModelRepresentations(self)
1246 self.all_protocols = _AllProtocols(self)
1247 self.all_datasets = _AllDatasets(self.system)
1248 self.all_starting_models = _AllStartingModels(self)
1250 self.all_software = _AllSoftware(self.system)
1253 """Create a new Representation and return it. This can be
1254 passed to add_model(), add_bead_element() or add_pdb_element()."""
1255 r = ihm.representation.Representation(name=name)
1256 self.system.orphan_representations.append(r)
1260 """Don't record coordinates for the given domain.
1261 Coordinates for the given domain (specified by a component name
1262 and a 2-element tuple giving the start and end residue numbers)
1263 will be excluded from the mmCIF file. This can be used to exclude
1264 parts of the structure that weren't well resolved in modeling.
1265 Any bead or residue that lies wholly within this range will be
1266 excluded. Multiple ranges for a given component can be excluded
1267 by calling this method multiple times."""
1268 if component
not in self._exclude_coords:
1269 self._exclude_coords[component] = []
1270 self._exclude_coords[component].append(seqrange)
1272 def _is_excluded(self, component, start, end):
1273 """Return True iff this chunk of sequence should be excluded"""
1274 for seqrange
in self._exclude_coords.get(component, ()):
1275 if start >= seqrange[0]
and end <= seqrange[1]:
1278 def _add_state(self, state):
1279 """Create a new state and return a pointer to it."""
1280 self._state_ensemble_offset = len(self.system.ensembles)
1281 s = _State(state, self)
1282 self._state_group.append(s)
1283 self._last_state = s
1286 def _get_chain_for_component(self, name, output):
1287 """Get the chain ID for a component, if any."""
1289 if name
in self.asym_units:
1290 return self.asym_units[name]._id
1295 def _get_assembly_comps(self, assembly):
1296 """Get the names of the components in the given assembly"""
1300 comps[ca.details] =
None
1304 """Make a new component that's a transformed copy of another.
1305 All representation for the existing component is copied to the
1307 assembly_comps = self._get_assembly_comps(state.modeled_assembly)
1308 if name
in assembly_comps:
1309 raise ValueError(
"Component %s already exists" % name)
1310 elif original
not in assembly_comps:
1311 raise ValueError(
"Original component %s does not exist" % original)
1312 self.create_component(state, name,
True)
1313 self.add_component_sequence(state, name, self.sequence_dict[original])
1314 self._transformed_components.append(_TransformedComponent(
1315 name, original, transform))
1316 self.all_representations.copy_component(state, name, original,
1317 self.asym_units[name])
1319 def create_component(self, state, name, modeled, asym_name=None):
1320 if asym_name
is None:
1322 new_comp = name
not in self._all_components
1323 self._all_components[name] =
None
1325 state.all_modeled_components.append(name)
1326 if asym_name
not in self.asym_units:
1328 self.asym_units[asym_name] =
None
1330 self.all_modeled_components.append(name)
1332 def add_component_sequence(self, state, name, seq, asym_name=None,
1334 if asym_name
is None:
1337 def get_offset(seq):
1339 for i
in range(len(seq)):
1342 seq, offset = get_offset(seq)
1343 if name
in self.sequence_dict:
1344 if self.sequence_dict[name] != seq:
1345 raise ValueError(
"Sequence mismatch for component %s" % name)
1347 self.sequence_dict[name] = seq
1348 self.entities.add(name, seq, offset, alphabet)
1349 if asym_name
in self.asym_units:
1350 if self.asym_units[asym_name]
is None:
1352 entity = self.entities[name]
1353 asym =
AsymUnit(entity, details=asym_name)
1354 self.system.asym_units.append(asym)
1355 self.asym_units[asym_name] = asym
1356 state.modeled_assembly.append(self.asym_units[asym_name])
1358 def _add_restraint_model_fits(self):
1359 """Add fits to restraints for all known models"""
1360 for group, m
in self.system._all_models():
1361 if m._is_restrained:
1362 for r
in self.system.restraints:
1363 if hasattr(r,
'add_fits_from_model_statfile'):
1364 r.add_fits_from_model_statfile(m)
1367 """Do any final processing on the class hierarchy.
1368 After calling this method, the `system` member (an instance
1369 of `ihm.System`) completely reproduces the PMI modeling, and
1370 can be written out to an mmCIF file with `ihm.dumper.write`,
1371 and/or modified using the ihm API."""
1372 self._add_restraint_model_fits()
1374 def add_pdb_element(self, state, name, start, end, offset, pdbname,
1375 chain, hier, representation=
None):
1376 if self._is_excluded(name, start, end):
1378 if representation
is None:
1379 representation = self.default_representation
1380 asym = self.asym_units[name]
1381 p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1383 self.all_representations.add_fragment(state, representation, p)
1384 self.all_starting_models.add_pdb_fragment(p)
1386 def add_bead_element(self, state, name, start, end, num, hier,
1387 representation=
None):
1388 if self._is_excluded(name, start, end):
1390 if representation
is None:
1391 representation = self.default_representation
1392 asym = self.asym_units[name]
1393 pmi_offset = asym.entity.pmi_offset
1394 b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1396 self.all_representations.add_fragment(state, representation, b)
1398 def get_cross_link_group(self, pmi_restraint):
1399 r = _CrossLinkRestraint(pmi_restraint)
1400 self.system.restraints.append(r)
1401 self._add_restraint_dataset(r)
1404 def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1405 if c1
not in self._all_components
or c2
not in self._all_components:
1411 e1 = self.entities[c1]
1412 e2 = self.entities[c2]
1413 xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1414 residue2=e2.pmi_residue(r2))
1415 rsr.experimental_cross_links.append([xl])
1418 def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1421 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1422 d = ihm.restraint.UpperBoundDistanceRestraint(length)
1424 if _get_by_residue(p1)
and _get_by_residue(p2):
1425 cls = _ResidueCrossLink
1427 cls = _FeatureCrossLink
1428 xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1431 xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1432 rsr.cross_links.append(xl)
1434 def add_replica_exchange(self, state, rex):
1439 step = _ReplicaExchangeProtocolStep(state, rex)
1440 step.software = self.all_software.pmi
1441 self.all_protocols.add_step(step, state)
1443 def _add_simple_dynamics(self, num_models_end, method):
1445 state = self._last_state
1446 self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1449 def _add_protocol(self):
1451 state = self._last_state
1452 self.all_protocols.add_protocol(state)
1454 def _add_dataset(self, dataset):
1455 return self.all_datasets.add(self._last_state, dataset)
1457 def _add_restraint_dataset(self, restraint):
1458 return self.all_datasets.add_restraint(self._last_state, restraint)
1460 def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1462 state = self._last_state
1463 pp = ihm.analysis.ClusterStep(
'RMSD', num_models_begin, num_models_end)
1464 self.all_protocols.add_postproc(pp, state)
1467 def _add_no_postprocessing(self, num_models):
1469 state = self._last_state
1470 pp = ihm.analysis.EmptyStep()
1471 pp.num_models_begin = pp.num_models_end = num_models
1472 self.all_protocols.add_postproc(pp, state)
1475 def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1476 num_models_deposited, localization_densities,
1478 """Add an ensemble generated by ad hoc methods (not using PMI).
1479 This is currently only used by the Nup84 system."""
1481 state = self._last_state
1482 group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1483 state.add_model_group(group)
1485 self.system.locations.append(ensemble_file)
1486 e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1488 self.system.ensembles.append(e)
1489 for c
in state.all_modeled_components:
1490 den = localization_densities.get(c,
None)
1492 e.load_localization_density(state, c, self.asym_units[c], den)
1496 """Point a previously-created ensemble to an 'all-models' file.
1497 This could be a trajectory such as DCD, an RMF, or a multimodel
1499 self.system.locations.append(location)
1501 ind = i + self._state_ensemble_offset
1502 self.system.ensembles[ind].file = location
1504 def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1510 protocol = self.all_protocols.get_last_protocol(state)
1511 num_models = protocol.steps[-1].num_models_end
1512 pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1513 pp.software = self.all_software.pmi
1514 self.all_protocols.add_postproc(pp, state)
1515 for i
in range(rex._number_of_clusters):
1516 group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1517 'cluster %d' % (i + 1)))
1518 state.add_model_group(group)
1520 e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1521 self.system.ensembles.append(e)
1523 for fname, stuple
in sorted(density_custom_ranges.items()):
1524 e.load_localization_density(state, fname, stuple,
1526 for stats
in e.load_all_models(self, state):
1527 m = self.add_model(group)
1530 m.name =
'Best scoring model'
1533 for c
in state.all_modeled_components:
1536 def _get_subassembly(self, comps, name, description):
1537 """Get an Assembly consisting of the given components.
1538 `compdict` is a dictionary of the components to add, where keys
1539 are the component names and values are the sequence ranges (or
1540 None to use all residues in the component)."""
1542 for comp, seqrng
in comps.items():
1543 a = self.asym_units[comp]
1544 asyms.append(a
if seqrng
is None else a(*seqrng))
1546 a = ihm.Assembly(asyms, name=name, description=description)
1549 def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1551 """Add a basic FoXS fit. This is largely intended for use from the
1553 assembly = self._get_subassembly(
1555 name=
"SAXS subassembly",
1556 description=
"All components that fit SAXS data")
1557 r = ihm.restraint.SASRestraint(
1558 dataset, assembly, segment=
False,
1559 fitting_method=
'FoXS', fitting_atom_type=
'Heavy atoms',
1560 multi_state=
False, radius_of_gyration=rg, details=details)
1561 r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1562 self.system.restraints.append(r)
1563 self._add_restraint_dataset(r)
1565 def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1566 image_resolution, projection_number,
1567 micrographs_number):
1568 r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1569 image_resolution, projection_number,
1571 self.system.restraints.append(r)
1572 self._add_restraint_dataset(r)
1574 def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1576 r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1577 self.system.restraints.append(r)
1578 self._add_restraint_dataset(r)
1580 def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1581 sigma, pmi_restraint):
1582 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1583 sigma, pmi_restraint, self._xy_plane)
1585 def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1586 sigma, pmi_restraint):
1587 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1588 sigma, pmi_restraint, self._xz_plane)
1590 def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1591 sigma, pmi_restraint):
1592 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1593 sigma, pmi_restraint, self._z_axis)
1595 def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1596 sigma, pmi_restraint, geom):
1597 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1598 r = _GeometricRestraint(
1599 self, state, pmi_restraint, geom, asym.get_feature(ps),
1600 ihm.restraint.LowerUpperBoundDistanceRestraint(lower_bound,
1603 self.system.restraints.append(r)
1604 self._add_restraint_dataset(r)
1606 def _get_membrane(self, tor_R, tor_r, tor_th):
1607 """Get an object representing a half-torus membrane"""
1608 if not hasattr(self,
'_seen_membranes'):
1609 self._seen_membranes = {}
1612 membrane_id = tuple(int(x * 100.)
for x
in (tor_R, tor_r, tor_th))
1613 if membrane_id
not in self._seen_membranes:
1614 m = ihm.geometry.HalfTorus(
1615 center=self._center_origin,
1616 transformation=self._identity_transform,
1617 major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1618 inner=
True, name=
'Membrane')
1619 self._seen_membranes[membrane_id] = m
1620 return self._seen_membranes[membrane_id]
1622 def add_membrane_surface_location_restraint(
1623 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1624 self._add_membrane_restraint(
1625 state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint,
1626 ihm.restraint.UpperBoundDistanceRestraint(0.))
1628 def add_membrane_exclusion_restraint(
1629 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1630 self._add_membrane_restraint(
1631 state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint,
1632 ihm.restraint.LowerBoundDistanceRestraint(0.))
1634 def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1635 sigma, pmi_restraint, rsr):
1636 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1637 r = _GeometricRestraint(
1638 self, state, pmi_restraint,
1639 self._get_membrane(tor_R, tor_r, tor_th), asym.get_feature(ps),
1641 self.system.restraints.append(r)
1642 self._add_restraint_dataset(r)
1644 def add_model(self, group, assembly=None, representation=None):
1645 state = self._last_state
1646 if representation
is None:
1647 representation = self.default_representation
1648 protocol = self.all_protocols.get_last_protocol(state)
1649 m = _Model(state.prot, self, protocol,
1650 assembly
if assembly
else state.modeled_assembly,
1657 """Get custom python-ihm dumpers for writing PMI to from mmCIF.
1658 This returns a list of custom dumpers that can be passed as all or
1659 part of the `dumpers` argument to ihm.dumper.write(). They add
1660 PMI-specific information to mmCIF or BinaryCIF files written out
1662 return [_ReplicaExchangeProtocolDumper]
1666 """Get custom python-ihm handlers for reading PMI data from mmCIF.
1667 This returns a list of custom handlers that can be passed as all or
1668 part of the `handlers` argument to ihm.reader.read(). They read
1669 PMI-specific information from mmCIF or BinaryCIF files read in
1671 return [_ReplicaExchangeProtocolHandler]
1675 """Extract metadata from an EM density GMM file."""
1678 """Extract metadata from `filename`.
1679 @return a dict with key `dataset` pointing to the GMM file and
1680 `number_of_gaussians` to the number of GMMs (or None)"""
1681 loc = ihm.location.InputFileLocation(
1683 details=
"Electron microscopy density map, "
1684 "represented as a Gaussian Mixture Model (GMM)")
1687 loc._allow_duplicates =
True
1688 d = ihm.dataset.EMDensityDataset(loc)
1689 ret = {
'dataset': d,
'number_of_gaussians':
None}
1691 with open(filename)
as fh:
1693 if line.startswith(
'# data_fn: '):
1694 p = ihm.metadata.MRCParser()
1695 fn = line[11:].rstrip(
'\r\n')
1696 dataset = p.parse_file(os.path.join(
1697 os.path.dirname(filename), fn))[
'dataset']
1698 ret[
'dataset'].parents.append(dataset)
1699 elif line.startswith(
'# ncenters: '):
1700 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.
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.
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...
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 ...