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.
38 import ihm.representation
40 import ihm.cross_linkers
43 def _assign_id(obj, seen_objs, obj_by_id):
44 """Assign a unique ID to obj, and track all ids in obj_by_id."""
45 if obj
not in seen_objs:
46 if not hasattr(obj,
'id'):
48 obj.id = len(obj_by_id)
49 seen_objs[obj] = obj.id
51 obj.id = seen_objs[obj]
54 def _get_by_residue(p):
55 """Determine whether the given particle represents a specific residue
56 or a more coarse-grained object."""
60 class _ComponentMapper:
61 """Map a Particle to a component name"""
62 def __init__(self, prot):
65 self.name =
'cif-output'
66 self.o.dictionary_pdbs[self.name] = self.prot
67 self.o._init_dictchain(self.name, self.prot,
68 multichar_chain=
True, mmcif=
True)
70 def __getitem__(self, p):
71 protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
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
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')
138 if hasattr(self.imp,
'citation'):
139 javi =
'Vel\u00e1zquez-Muriel J'
140 self.imp.citation = ihm.Citation(
142 title=
'Putting the pieces together: integrative modeling '
143 'platform software for structure determination of '
144 'macromolecular assemblies',
145 journal=
'PLoS Biol', volume=10, page_range=
'e1001244',
147 authors=[
'Russel D',
'Lasker K',
'Webb B', javi,
'Tjioe E',
148 'Schneidman-Duhovny D',
'Peterson B',
'Sali A'],
149 doi=
'10.1371/journal.pbio.1001244')
150 self.pmi.citation = ihm.Citation(
152 title=
'Modeling Biological Complexes Using Integrative '
153 'Modeling Platform.',
154 journal=
'Methods Mol Biol', volume=2022, page_range=(353, 377),
156 authors=[
'Saltzberg D',
'Greenberg CH',
'Viswanath S',
157 'Chemmama I',
'Webb B',
'Pellarin R',
'Echeverria I',
159 doi=
'10.1007/978-1-4939-9608-7_15')
160 self.system.software.extend([self.pmi, self.imp])
162 def set_modeller_used(self, version, date):
163 if self.modeller_used:
165 self.modeller_used =
True
167 name=
'MODELLER', classification=
'comparative modeling',
168 description=
'Comparative modeling by satisfaction '
169 'of spatial restraints, build ' + date,
170 location=
'https://salilab.org/modeller/', version=version)
171 self.system.software.append(s)
172 if hasattr(s,
'citation'):
173 s.citation = ihm.Citation(
175 title=
'Comparative protein modelling by satisfaction of '
176 'spatial restraints.',
177 journal=
'J Mol Biol', volume=234, page_range=(779, 815),
178 year=1993, authors=[
'Sali A',
'Blundell TL'],
179 doi=
'10.1006/jmbi.1993.1626')
181 def set_phyre2_used(self):
184 self.phyre2_used =
True
186 name=
'Phyre2', classification=
'protein homology modeling',
187 description=
'Protein Homology/analogY Recognition Engine V 2.0',
188 version=
'2.0', location=
'http://www.sbg.bio.ic.ac.uk/~phyre2/')
189 if hasattr(s,
'citation'):
190 s.citation = ihm.Citation(
192 title=
'The Phyre2 web portal for protein modeling, '
193 'prediction and analysis.',
194 journal=
'Nat Protoc', volume=10, page_range=(845, 858),
195 authors=[
'Kelley LA',
'Mezulis S',
'Yates CM',
'Wass MN',
197 year=2015, doi=
'10.1038/nprot.2015.053')
198 self.system.software.append(s)
201 def _get_fragment_is_rigid(fragment):
202 """Determine whether a fragment is modeled rigidly"""
208 class _PDBFragment(ihm.representation.ResidueSegment):
209 """Record details about part of a PDB file used as input
211 def __init__(self, state, component, start, end, pdb_offset,
212 pdbname, chain, hier, asym_unit):
215 asym_unit=asym_unit.pmi_range(start, end),
216 rigid=
None, primitive=
'sphere')
217 self.component, self.start, self.end, self.offset, self.pdbname \
218 = component, start, end, pdb_offset, pdbname
219 self.state, self.chain, self.hier = state, chain, hier
224 rigid = property(
lambda self: _get_fragment_is_rigid(self),
225 lambda self, val:
None)
227 def combine(self, other):
231 class _BeadsFragment(ihm.representation.FeatureSegment):
232 """Record details about beads used to represent part of a component."""
235 def __init__(self, state, component, start, end, count, hier, asym_unit):
237 asym_unit=asym_unit(start, end), rigid=
None, primitive=
'sphere',
239 self.state, self.component, self.hier = state, component, hier
241 rigid = property(
lambda self: _get_fragment_is_rigid(self),
242 lambda self, val:
None)
244 def combine(self, other):
246 if (type(other) == type(self)
and
247 other.asym_unit.seq_id_range[0]
248 == self.asym_unit.seq_id_range[1] + 1):
249 self.asym_unit.seq_id_range = (self.asym_unit.seq_id_range[0],
250 other.asym_unit.seq_id_range[1])
251 self.count += other.count
255 class _AllModelRepresentations:
256 def __init__(self, simo):
260 self.fragments = OrderedDict()
261 self._all_representations = {}
263 def copy_component(self, state, name, original, asym_unit):
264 """Copy all representation for `original` in `state` to `name`"""
267 newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
269 for rep
in self.fragments:
270 if original
in self.fragments[rep]:
271 if name
not in self.fragments[rep]:
272 self.fragments[rep][name] = OrderedDict()
273 self.fragments[rep][name][state] = [
274 copy_frag(f)
for f
in self.fragments[rep][original][state]]
277 first_state = list(self.fragments[rep][name].keys())[0]
278 if state
is first_state:
279 representation = self._all_representations[rep]
280 representation.extend(self.fragments[rep][name][state])
282 def add_fragment(self, state, representation, fragment):
283 """Add a model fragment."""
284 comp = fragment.component
285 id_rep = id(representation)
286 self._all_representations[id_rep] = representation
287 if id_rep
not in self.fragments:
288 self.fragments[id_rep] = OrderedDict()
289 if comp
not in self.fragments[id_rep]:
290 self.fragments[id_rep][comp] = OrderedDict()
291 if state
not in self.fragments[id_rep][comp]:
292 self.fragments[id_rep][comp][state] = []
293 fragments = self.fragments[id_rep][comp][state]
294 if len(fragments) == 0
or not fragments[-1].combine(fragment):
295 fragments.append(fragment)
298 first_state = list(self.fragments[id_rep][comp].keys())[0]
299 if state
is first_state:
300 representation.append(fragment)
304 """Track all datasets generated by PMI and add them to the ihm.System"""
305 def __init__(self, system):
308 self._datasets_by_state = {}
309 self._restraints_by_state = {}
311 def get_all_group(self, state):
312 """Get a DatasetGroup encompassing all datasets so far in this state"""
316 g = ihm.dataset.DatasetGroup(
317 self._datasets_by_state.get(state, [])
318 + [r.dataset
for r
in self._restraints_by_state.get(state, [])
322 def add(self, state, dataset):
323 """Add a new dataset."""
324 self._datasets.append(dataset)
325 if state
not in self._datasets_by_state:
326 self._datasets_by_state[state] = []
327 self._datasets_by_state[state].append(dataset)
330 self.system.orphan_datasets.append(dataset)
334 """Add the dataset for a restraint"""
335 if state
not in self._restraints_by_state:
336 self._restraints_by_state[state] = []
337 self._restraints_by_state[state].append(restraint)
340 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
341 """Restrain to a set of cross-links"""
344 _label_map = {
'wtDSS':
'DSS',
'scDSS':
'DSS',
'scEDC':
'EDC'}
345 _descriptor_map = {
'DSS': ihm.cross_linkers.dss,
346 'EDC': ihm.cross_linkers.edc}
348 def __init__(self, pmi_restraint):
349 self.pmi_restraint = pmi_restraint
352 linker = getattr(self.pmi_restraint,
'linker',
None)
353 label = self.pmi_restraint.label
356 dataset=self.pmi_restraint.dataset,
357 linker=linker
or self._get_chem_descriptor(label))
360 def _get_chem_descriptor(cls, label):
362 label = cls._label_map.get(label, label)
363 if label
not in cls._descriptor_map:
367 d = ihm.ChemDescriptor(label)
368 cls._descriptor_map[label] = d
369 return cls._descriptor_map[label]
371 def _set_psi_sigma(self, model):
374 if model.m != self.pmi_restraint.model:
376 for resolution
in self.pmi_restraint.sigma_dictionary:
377 statname =
'ISDCrossLinkMS_Sigma_%s_%s' % (resolution, self.label)
378 if model.stats
and statname
in model.stats:
379 sigma = float(model.stats[statname])
380 p = self.pmi_restraint.sigma_dictionary[resolution][0]
381 old_values.append((p, p.get_scale()))
383 for psiindex
in self.pmi_restraint.psi_dictionary:
384 statname =
'ISDCrossLinkMS_Psi_%s_%s' % (psiindex, self.label)
385 if model.stats
and statname
in model.stats:
386 psi = float(model.stats[statname])
387 p = self.pmi_restraint.psi_dictionary[psiindex][0]
388 old_values.append((p, p.get_scale()))
392 return list(reversed(old_values))
394 def add_fits_from_model_statfile(self, model):
396 old_values = self._set_psi_sigma(model)
400 for xl
in self.cross_links:
402 xl.fits[model] = ihm.restraint.CrossLinkFit(
403 psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
405 for p, scale
in old_values:
409 def __set_dataset(self, val):
410 self.pmi_restraint.dataset = val
411 dataset = property(
lambda self: self.pmi_restraint.dataset,
415 def get_asym_mapper_for_state(simo, state, asym_map):
416 asym = asym_map.get(state,
None)
418 asym = _AsymMapper(simo, state.prot)
419 asym_map[state] = asym
426 psi = property(
lambda self: self.psi_p.get_scale(),
427 lambda self, val:
None)
428 sigma1 = property(
lambda self: self.sigma1_p.get_scale(),
429 lambda self, val:
None)
430 sigma2 = property(
lambda self: self.sigma2_p.get_scale(),
431 lambda self, val:
None)
434 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
438 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
442 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
443 def __init__(self, state, pmi_restraint, image_number, resolution,
444 pixel_size, image_resolution, projection_number,
446 self.pmi_restraint, self.image_number = pmi_restraint, image_number
448 dataset=pmi_restraint.datasets[image_number],
449 assembly=state.modeled_assembly,
450 segment=
False, number_raw_micrographs=micrographs_number,
451 pixel_size_width=pixel_size, pixel_size_height=pixel_size,
452 image_resolution=image_resolution,
453 number_of_projections=projection_number)
456 def __get_dataset(self):
457 return self.pmi_restraint.datasets[self.image_number]
459 def __set_dataset(self, val):
460 self.pmi_restraint.datasets[self.image_number] = val
462 dataset = property(__get_dataset, __set_dataset)
464 def add_fits_from_model_statfile(self, model):
465 ccc = self._get_cross_correlation(model)
466 transform = self._get_transformation(model)
467 rot = transform.get_rotation()
468 rm = [[e
for e
in rot.get_rotation_matrix_row(i)]
for i
in range(3)]
469 self.fits[model] = ihm.restraint.EM2DRestraintFit(
470 cross_correlation_coefficient=ccc,
472 tr_vector=transform.get_translation())
474 def _get_transformation(self, model):
475 """Get the transformation that places the model on the image"""
476 stats = model.em2d_stats
or model.stats
477 prefix =
'ElectronMicroscopy2D_%s_Image%d' % (self.pmi_restraint.label,
478 self.image_number + 1)
479 r = [float(stats[prefix +
'_Rotation%d' % i])
for i
in range(4)]
480 t = [float(stats[prefix +
'_Translation%d' % i])
484 inv = model.transform.get_inverse()
488 def _get_cross_correlation(self, model):
489 """Get the cross correlation coefficient between the model projection
491 stats = model.em2d_stats
or model.stats
492 return float(stats[
'ElectronMicroscopy2D_%s_Image%d_CCC'
493 % (self.pmi_restraint.label,
494 self.image_number + 1)])
497 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
499 def __init__(self, simo, state, pmi_restraint, target_ps, densities):
500 self.pmi_restraint = pmi_restraint
502 dataset=pmi_restraint.dataset,
503 assembly=self._get_assembly(densities, simo, state),
504 fitting_method=
'Gaussian mixture models',
505 number_of_gaussians=len(target_ps))
508 def __set_dataset(self, val):
509 self.pmi_restraint.dataset = val
510 dataset = property(
lambda self: self.pmi_restraint.dataset,
513 def _get_assembly(self, densities, simo, state):
514 """Get the Assembly that this restraint acts on"""
515 cm = _ComponentMapper(state.prot)
518 components[cm[d]] =
None
519 a = simo._get_subassembly(
520 components, name=
"EM subassembly",
521 description=
"All components that fit the EM map")
524 def add_fits_from_model_statfile(self, model):
525 ccc = self._get_cross_correlation(model)
526 self.fits[model] = ihm.restraint.EM3DRestraintFit(
527 cross_correlation_coefficient=ccc)
529 def _get_cross_correlation(self, model):
530 """Get the cross correlation coefficient between the model
532 if model.stats
is not None:
533 return float(model.stats[
'GaussianEMRestraint_%s_CCC'
534 % self.pmi_restraint.label])
537 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
539 def __init__(self, simo, state, pmi_restraint, geometric_object,
540 feature, distance, sigma):
541 self.pmi_restraint = pmi_restraint
543 dataset=pmi_restraint.dataset,
544 geometric_object=geometric_object, feature=feature,
545 distance=distance, harmonic_force_constant=1. / sigma,
549 def __set_dataset(self, val):
550 self.pmi_restraint.dataset = val
551 dataset = property(
lambda self: self.pmi_restraint.dataset,
555 class _ReplicaExchangeProtocolStep(ihm.protocol.Step):
556 def __init__(self, state, rex):
557 if rex.monte_carlo_sample_objects
is not None:
558 method =
'Replica exchange monte carlo'
560 method =
'Replica exchange molecular dynamics'
561 self.monte_carlo_temperature = rex.vars[
'monte_carlo_temperature']
562 self.replica_exchange_minimum_temperature = \
563 rex.vars[
'replica_exchange_minimum_temperature']
564 self.replica_exchange_maximum_temperature = \
565 rex.vars[
'replica_exchange_maximum_temperature']
567 assembly=state.modeled_assembly,
569 method=method, name=
'Sampling',
570 num_models_begin=
None,
571 num_models_end=rex.vars[
"number_of_frames"],
572 multi_scale=
True, multi_state=
False, ordered=
False, ensemble=
True)
575 class _ReplicaExchangeProtocolDumper(ihm.dumper.Dumper):
576 """Write IMP-specific information about replica exchange to mmCIF.
577 Note that IDs will have already been assigned by python-ihm's
578 standard modeling protocol dumper."""
579 def dump(self, system, writer):
580 with writer.loop(
"_imp_replica_exchange_protocol",
581 [
"protocol_id",
"step_id",
"monte_carlo_temperature",
582 "replica_exchange_minimum_temperature",
583 "replica_exchange_maximum_temperature"])
as lp:
584 for p
in system._all_protocols():
586 if isinstance(s, _ReplicaExchangeProtocolStep):
587 self._dump_step(p, s, lp)
589 def _dump_step(self, p, s, lp):
590 mintemp = s.replica_exchange_minimum_temperature
591 maxtemp = s.replica_exchange_maximum_temperature
592 lp.write(protocol_id=p._id, step_id=s._id,
593 monte_carlo_temperature=s.monte_carlo_temperature,
594 replica_exchange_minimum_temperature=mintemp,
595 replica_exchange_maximum_temperature=maxtemp)
598 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
599 category =
'_imp_replica_exchange_protocol'
601 """Read IMP-specific information about replica exchange from mmCIF."""
602 def __call__(self, protocol_id, step_id, monte_carlo_temperature,
603 replica_exchange_minimum_temperature,
604 replica_exchange_maximum_temperature):
605 p = self.sysr.protocols.get_by_id(protocol_id)
607 s = p.steps[int(step_id)-1]
609 s.__class__ = _ReplicaExchangeProtocolStep
610 s.monte_carlo_temperature = \
611 self.get_float(monte_carlo_temperature)
612 s.replica_exchange_minimum_temperature = \
613 self.get_float(replica_exchange_minimum_temperature)
614 s.replica_exchange_maximum_temperature = \
615 self.get_float(replica_exchange_maximum_temperature)
618 class _SimpleProtocolStep(ihm.protocol.Step):
619 def __init__(self, state, num_models_end, method):
621 assembly=state.modeled_assembly,
623 method=method, name=
'Sampling',
624 num_models_begin=
None,
625 num_models_end=num_models_end,
626 multi_scale=
True, multi_state=
False, ordered=
False,
631 """Represent a single chain in a Model"""
632 def __init__(self, pmi_chain_id, asym_unit):
633 self.pmi_chain_id, self.asym_unit = pmi_chain_id, asym_unit
637 def add(self, xyz, atom_type, residue_type, residue_index,
638 all_indexes, radius):
639 if atom_type
is None:
640 self.spheres.append((xyz, residue_type, residue_index,
641 all_indexes, radius))
643 self.atoms.append((xyz, atom_type, residue_type, residue_index,
644 all_indexes, radius))
645 orig_comp = property(
lambda self: self.comp)
648 class _TransformedChain:
649 """Represent a chain that is a transformed version of another"""
650 def __init__(self, orig_chain, asym_unit, transform):
651 self.orig_chain, self.asym_unit = orig_chain, asym_unit
652 self.transform = transform
654 def __get_spheres(self):
655 for (xyz, residue_type, residue_index, all_indexes,
656 radius)
in self.orig_chain.spheres:
657 yield (self.transform * xyz, residue_type, residue_index,
659 spheres = property(__get_spheres)
661 def __get_atoms(self):
662 for (xyz, atom_type, residue_type, residue_index, all_indexes,
663 radius)
in self.orig_chain.atoms:
664 yield (self.transform * xyz, atom_type, residue_type,
665 residue_index, all_indexes, radius)
666 atoms = property(__get_atoms)
668 entity = property(
lambda self: self.orig_chain.entity)
669 orig_comp = property(
lambda self: self.orig_chain.comp)
673 def __init__(self, component, simo):
674 self._seqranges = simo._exclude_coords.get(component, [])
676 def is_excluded(self, indexes):
677 """Return True iff the given sequence range is excluded."""
678 for seqrange
in self._seqranges:
679 if indexes[0] >= seqrange[0]
and indexes[-1] <= seqrange[1]:
683 class _Model(ihm.model.Model):
684 def __init__(self, prot, simo, protocol, assembly, representation):
685 super().__init__(assembly=assembly, protocol=protocol,
686 representation=representation)
687 self.simo = weakref.proxy(simo)
693 self.em2d_stats =
None
696 self._is_restrained =
True
699 self.m = prot.get_model()
700 o.dictionary_pdbs[name] = prot
701 o._init_dictchain(name, prot, multichar_chain=
True)
702 (particle_infos_for_pdb,
703 self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
705 self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
708 def all_chains(self, simo):
709 """Yield all chains, including transformed ones"""
711 for c
in self.chains:
713 chain_for_comp[c.comp] = c
714 for tc
in simo._transformed_components:
715 orig_chain = chain_for_comp.get(tc.original,
None)
717 asym = simo.asym_units[tc.name]
718 c = _TransformedChain(orig_chain, asym, tc.transform)
722 def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
723 entity_for_chain = {}
726 for protname, chain_id
in o.dictchain[name].items():
727 if protname
in simo.entities:
728 entity_for_chain[chain_id] = simo.entities[protname]
731 pn = protname.split(
'.')[0]
732 entity_for_chain[chain_id] = simo.entities[pn]
733 comp_for_chain[chain_id] = protname
737 correct_asym[chain_id] = simo.asym_units[protname]
744 for (xyz, atom_type, residue_type, chain_id, residue_index,
745 all_indexes, radius)
in particle_infos_for_pdb:
746 if chain
is None or chain.pmi_chain_id != chain_id:
747 chain = _Chain(chain_id, correct_asym[chain_id])
748 chain.entity = entity_for_chain[chain_id]
749 chain.comp = comp_for_chain[chain_id]
750 self.chains.append(chain)
751 excluder = _Excluder(chain.comp, simo)
752 if not excluder.is_excluded(all_indexes
if all_indexes
753 else [residue_index]):
754 chain.add(xyz, atom_type, residue_type, residue_index,
757 def parse_rmsf_file(self, fname, component):
758 self.rmsf[component] = rmsf = {}
759 with open(str(fname))
as fh:
761 resnum, blocknum, val = line.split()
762 rmsf[int(resnum)] = (int(blocknum), float(val))
764 def get_rmsf(self, component, indexes):
765 """Get the RMSF value for the given residue indexes."""
768 rmsf = self.rmsf[component]
769 blocknums = dict.fromkeys(rmsf[ind][0]
for ind
in indexes)
770 if len(blocknums) != 1:
771 raise ValueError(
"Residue indexes %s aren't all in the same block"
773 return rmsf[indexes[0]][1]
776 for chain
in self.all_chains(self.simo):
777 pmi_offset = chain.asym_unit.entity.pmi_offset
778 for atom
in chain.atoms:
779 (xyz, atom_type, residue_type, residue_index,
780 all_indexes, radius) = atom
781 pt = self.transform * xyz
782 yield ihm.model.Atom(
783 asym_unit=chain.asym_unit,
784 seq_id=residue_index - pmi_offset,
785 atom_id=atom_type.get_string(),
787 x=pt[0], y=pt[1], z=pt[2])
789 def get_spheres(self):
790 for chain
in self.all_chains(self.simo):
791 pmi_offset = chain.asym_unit.entity.pmi_offset
792 for sphere
in chain.spheres:
793 (xyz, residue_type, residue_index,
794 all_indexes, radius) = sphere
795 if all_indexes
is None:
796 all_indexes = (residue_index,)
797 pt = self.transform * xyz
798 yield ihm.model.Sphere(
799 asym_unit=chain.asym_unit,
800 seq_id_range=(all_indexes[0] - pmi_offset,
801 all_indexes[-1] - pmi_offset),
802 x=pt[0], y=pt[1], z=pt[2], radius=radius,
803 rmsf=self.get_rmsf(chain.orig_comp, all_indexes))
807 def __init__(self, simo):
808 self.simo = weakref.proxy(simo)
810 self.protocols = OrderedDict()
812 def add_protocol(self, state):
813 """Add a new Protocol"""
814 if state
not in self.protocols:
815 self.protocols[state] = []
816 p = ihm.protocol.Protocol()
817 self.simo.system.orphan_protocols.append(p)
818 self.protocols[state].append(p)
820 def add_step(self, step, state):
821 """Add a ProtocolStep to the last Protocol of the given State"""
822 if state
not in self.protocols:
823 self.add_protocol(state)
824 protocol = self.get_last_protocol(state)
825 if len(protocol.steps) == 0:
826 step.num_models_begin = 0
828 step.num_models_begin = protocol.steps[-1].num_models_end
829 protocol.steps.append(step)
830 step.id = len(protocol.steps)
832 step.dataset_group = self.simo.all_datasets.get_all_group(state)
834 def add_postproc(self, step, state):
835 """Add a postprocessing step to the last protocol"""
836 protocol = self.get_last_protocol(state)
837 if len(protocol.analyses) == 0:
838 protocol.analyses.append(ihm.analysis.Analysis())
839 protocol.analyses[-1].steps.append(step)
841 def get_last_protocol(self, state):
842 """Return the most recently-added _Protocol"""
843 return self.protocols[state][-1]
846 class _AllStartingModels:
847 def __init__(self, simo):
851 self.models = OrderedDict()
854 def add_pdb_fragment(self, fragment):
855 """Add a starting model PDB fragment."""
856 comp = fragment.component
857 state = fragment.state
858 if comp
not in self.models:
859 self.models[comp] = OrderedDict()
860 if state
not in self.models[comp]:
861 self.models[comp][state] = []
862 models = self.models[comp][state]
863 if len(models) == 0 \
864 or models[-1].fragments[0].pdbname != fragment.pdbname:
865 model = self._add_model(fragment)
869 models[-1].fragments.append(weakref.proxy(fragment))
873 pmi_offset = models[-1].asym_unit.entity.pmi_offset
874 sid_begin = min(fragment.start - pmi_offset,
875 models[-1].asym_unit.seq_id_range[0])
876 sid_end = max(fragment.end - pmi_offset,
877 models[-1].asym_unit.seq_id_range[1])
878 models[-1].asym_unit = fragment.asym_unit.asym(sid_begin, sid_end)
879 fragment.starting_model = models[-1]
881 def _add_model(self, f):
882 parser = ihm.metadata.PDBParser()
883 r = parser.parse_file(f.pdbname)
885 self.simo._add_dataset(r[
'dataset'])
887 templates = r[
'templates'].get(f.chain, [])
890 self.simo.system.locations.append(t.alignment_file)
892 self.simo._add_dataset(t.dataset)
893 source = r[
'entity_source'].get(f.chain)
895 f.asym_unit.entity.source = source
896 pmi_offset = f.asym_unit.entity.pmi_offset
898 asym_unit=f.asym_unit.asym.pmi_range(f.start, f.end),
899 dataset=r[
'dataset'], asym_id=f.chain,
900 templates=templates, offset=f.offset + pmi_offset,
901 metadata=r[
'metadata'],
902 software=r[
'software'][0]
if r[
'software']
else None,
903 script_file=r[
'script'])
904 m.fragments = [weakref.proxy(f)]
908 class _StartingModel(ihm.startmodel.StartingModel):
909 def get_seq_dif(self):
913 pmi_offset = self.asym_unit.entity.pmi_offset
914 mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
916 for f
in self.fragments:
919 residue_indexes=list(range(f.start - f.offset,
920 f.end - f.offset + 1)))
921 for a
in mh.get_ihm_atoms(sel.get_selected_particles(),
922 f.offset - pmi_offset):
924 self._seq_dif = mh._seq_dif
927 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
928 """Post processing using AnalysisReplicaExchange0 macro"""
930 def __init__(self, rex, num_models_begin):
933 for fname
in self.get_all_stat_files():
934 with open(str(fname))
as fh:
935 num_models_end += len(fh.readlines())
937 feature=
'RMSD', num_models_begin=num_models_begin,
938 num_models_end=num_models_end)
940 def get_stat_file(self, cluster_num):
941 return self.rex._outputdir / (
"cluster.%d" % cluster_num) /
'stat.out'
943 def get_all_stat_files(self):
944 for i
in range(self.rex._number_of_clusters):
945 yield self.get_stat_file(i)
948 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
949 """Ensemble generated using AnalysisReplicaExchange0 macro"""
951 num_models_deposited =
None
953 def __init__(self, pp, cluster_num, model_group, num_deposit):
954 with open(str(pp.get_stat_file(cluster_num)))
as fh:
955 num_models = len(fh.readlines())
957 num_models=num_models,
958 model_group=model_group, post_process=pp,
959 clustering_feature=pp.feature,
960 name=model_group.name)
961 self.cluster_num = cluster_num
962 self.num_models_deposited = num_deposit
964 def get_rmsf_file(self, component):
965 return (self.post_process.rex._outputdir
966 / (
'cluster.%d' % self.cluster_num)
967 / (
'rmsf.%s.dat' % component))
969 def load_rmsf(self, model, component):
970 fname = self.get_rmsf_file(component)
972 model.parse_rmsf_file(fname, component)
974 def get_localization_density_file(self, fname):
975 return (self.post_process.rex._outputdir
976 / (
'cluster.%d' % self.cluster_num)
977 / (
'%s.mrc' % fname))
979 def load_localization_density(self, state, fname, select_tuple,
981 fullpath = self.get_localization_density_file(fname)
982 if fullpath.exists():
983 details =
"Localization density for %s %s" \
984 % (fname, self.model_group.name)
985 local_file = ihm.location.OutputFileLocation(str(fullpath),
987 for s
in select_tuple:
988 if isinstance(s, tuple)
and len(s) == 3:
989 asym = asym_units[s[2]].pmi_range(s[0], s[1])
992 den = ihm.model.LocalizationDensity(file=local_file,
994 self.densities.append(den)
996 def load_all_models(self, simo, state):
997 stat_fname = self.post_process.get_stat_file(self.cluster_num)
999 with open(str(stat_fname))
as fh:
1000 stats = ast.literal_eval(fh.readline())
1002 rmf_file = stat_fname.parent / (
"%d.rmf3" % model_num)
1004 if rmf_file.exists():
1005 rh = RMF.open_rmf_file_read_only(str(rmf_file))
1006 system = state._pmi_object.system
1012 if model_num >= self.num_models_deposited:
1016 def _get_precision(self):
1017 precfile = (self.post_process.rex._outputdir /
1018 (
"precision.%d.%d.out" % (self.cluster_num,
1020 if not precfile.exists():
1024 r'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
1026 with open(str(precfile))
as fh:
1030 return float(m.group(1))
1032 precision = property(
lambda self: self._get_precision(),
1033 lambda self, val:
None)
1036 class _SimpleEnsemble(ihm.model.Ensemble):
1037 """Simple manually-created ensemble"""
1039 num_models_deposited =
None
1041 def __init__(self, pp, model_group, num_models, drmsd,
1042 num_models_deposited, ensemble_file):
1044 model_group=model_group, post_process=pp, num_models=num_models,
1045 file=ensemble_file, precision=drmsd, name=model_group.name,
1046 clustering_feature=
'dRMSD')
1047 self.num_models_deposited = num_models_deposited
1049 def load_localization_density(self, state, component, asym, local_file):
1050 den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
1051 self.densities.append(den)
1054 class _CustomDNAAlphabet:
1055 """Custom DNA alphabet that maps A,C,G,T (rather than DA,DC,DG,DT
1056 as in python-ihm)"""
1057 _comps = dict([cc.code_canonical, cc]
1058 for cc
in ihm.DNAAlphabet._comps.values())
1061 class _EntityMapper(dict):
1062 """Handle mapping from IMP components (without copy number) to CIF
1063 entities. Multiple components may map to the same entity if they
1065 def __init__(self, system):
1067 self._sequence_dict = {}
1069 self.system = system
1071 def _get_alphabet(self, alphabet):
1072 """Map a PMI alphabet to an IHM alphabet"""
1076 alphabet_map = {
None: ihm.LPeptideAlphabet,
1077 IMP.pmi.alphabets.amino_acid: ihm.LPeptideAlphabet,
1078 IMP.pmi.alphabets.rna: ihm.RNAAlphabet,
1079 IMP.pmi.alphabets.dna: _CustomDNAAlphabet}
1080 if alphabet
in alphabet_map:
1081 return alphabet_map[alphabet]
1083 raise TypeError(
"Don't know how to handle %s" % alphabet)
1085 def add(self, component_name, sequence, offset, alphabet):
1086 def entity_seq(sequence):
1089 return [
'UNK' if s ==
'X' else s
for s
in sequence]
1092 if sequence
not in self._sequence_dict:
1095 d = component_name.split(
"@")[0].split(
".")[0]
1096 entity = Entity(entity_seq(sequence), description=d,
1098 alphabet=self._get_alphabet(alphabet))
1099 self.system.entities.append(entity)
1100 self._sequence_dict[sequence] = entity
1101 self[component_name] = self._sequence_dict[sequence]
1104 class _TransformedComponent:
1105 def __init__(self, name, original, transform):
1106 self.name, self.original, self.transform = name, original, transform
1110 """Class with similar interface to weakref.ref, but keeps a strong ref"""
1111 def __init__(self, ref):
1118 class _State(ihm.model.State):
1119 """Representation of a single state in the system."""
1121 def __init__(self, pmi_object, po):
1126 self._pmi_object = weakref.proxy(pmi_object)
1127 if hasattr(pmi_object,
'state'):
1130 self._pmi_state = _SimpleRef(pmi_object.state)
1132 self._pmi_state = weakref.ref(pmi_object)
1134 old_name = self.name
1135 super().__init__(experiment_type=
'Fraction of bulk')
1136 self.name = old_name
1140 self.modeled_assembly = ihm.Assembly(
1141 name=
"Modeled assembly",
1142 description=self.get_postfixed_name(
1143 "All components modeled by IMP"))
1144 po.system.orphan_assemblies.append(self.modeled_assembly)
1146 self.all_modeled_components = []
1149 return hash(self._pmi_state())
1151 def __eq__(self, other):
1152 return self._pmi_state() == other._pmi_state()
1154 def add_model_group(self, group):
1158 def get_prefixed_name(self, name):
1159 """Prefix the given name with the state name, if available."""
1161 return self.short_name +
' ' + name
1165 return name[0].upper() + name[1:]
if name
else ''
1167 def get_postfixed_name(self, name):
1168 """Postfix the given name with the state name, if available."""
1170 return "%s in state %s" % (name, self.short_name)
1174 short_name = property(
lambda self: self._pmi_state().short_name)
1175 long_name = property(
lambda self: self._pmi_state().long_name)
1177 def __get_name(self):
1178 return self._pmi_state().long_name
1180 def __set_name(self, val):
1181 self._pmi_state().long_name = val
1183 name = property(__get_name, __set_name)
1187 """A single entity in the system.
1188 This functions identically to the base ihm.Entity class, but it
1189 allows identifying residues by either the IHM numbering scheme
1190 (seq_id, which is always contiguous starting at 1) or by the PMI
1191 scheme (which is similar, but may not start at 1 if the sequence in
1192 the FASTA file has one or more N-terminal gaps"""
1193 def __init__(self, sequence, pmi_offset, *args, **keys):
1196 self.pmi_offset = pmi_offset
1197 super().__init__(sequence, *args, **keys)
1200 """Return a single IHM residue indexed using PMI numbering"""
1201 return self.residue(res_id - self.pmi_offset)
1204 """Return a range of IHM residues indexed using PMI numbering"""
1205 off = self.pmi_offset
1206 return self(res_id_begin - off, res_id_end - off)
1210 """A single asymmetric unit in the system.
1211 This functions identically to the base ihm.AsymUnit class, but it
1212 allows identifying residues by either the IHM numbering scheme
1213 (seq_id, which is always contiguous starting at 1) or by the PMI
1214 scheme (which is similar, but may not start at 1 if the sequence in
1215 the FASTA file has one or more N-terminal gaps"""
1217 def __init__(self, entity, *args, **keys):
1219 entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1222 """Return a single IHM residue indexed using PMI numbering"""
1223 return self.residue(res_id - self.entity.pmi_offset)
1226 """Return a range of IHM residues indexed using PMI numbering"""
1227 off = self.entity.pmi_offset
1228 return self(res_id_begin - off, res_id_end - off)
1232 """Class to encode a modeling protocol as mmCIF.
1234 IMP has basic support for writing out files in mmCIF format, for
1235 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1236 After creating an instance of this class, attach it to an
1237 IMP.pmi.topology.System object. After this, any
1238 generated models and metadata are automatically collected in the
1239 `system` attribute, which is an
1240 [ihm.System](https://python-ihm.readthedocs.io/en/latest/main.html#ihm.System) object.
1241 Once the protocol is complete, call finalize() to make sure `system`
1242 contains everything, then use the
1243 [python-ihm API](https://python-ihm.readthedocs.io/en/latest/dumper.html#ihm.dumper.write)
1244 to write out files in mmCIF or BinaryCIF format.
1246 See also get_handlers(), get_dumpers().
1250 self.system = ihm.System()
1251 self._state_group = ihm.model.StateGroup()
1252 self.system.state_groups.append(self._state_group)
1254 self._state_ensemble_offset = 0
1255 self._main_script = os.path.abspath(sys.argv[0])
1258 loc = ihm.location.WorkflowFileLocation(
1259 path=self._main_script,
1260 details=
"The main integrative modeling script")
1261 self.system.locations.append(loc)
1264 self.__asym_states = {}
1265 self._working_directory = os.getcwd()
1267 "Default representation")
1268 self.entities = _EntityMapper(self.system)
1270 self.asym_units = {}
1271 self._all_components = {}
1272 self.all_modeled_components = []
1273 self._transformed_components = []
1274 self.sequence_dict = {}
1277 self._xy_plane = ihm.geometry.XYPlane()
1278 self._xz_plane = ihm.geometry.XZPlane()
1279 self._z_axis = ihm.geometry.ZAxis()
1280 self._center_origin = ihm.geometry.Center(0, 0, 0)
1281 self._identity_transform = ihm.geometry.Transformation.identity()
1284 self._exclude_coords = {}
1286 self.all_representations = _AllModelRepresentations(self)
1287 self.all_protocols = _AllProtocols(self)
1288 self.all_datasets = _AllDatasets(self.system)
1289 self.all_starting_models = _AllStartingModels(self)
1291 self.all_software = _AllSoftware(self.system)
1294 """Create a new Representation and return it. This can be
1295 passed to add_model(), add_bead_element() or add_pdb_element()."""
1296 r = ihm.representation.Representation(name=name)
1297 self.system.orphan_representations.append(r)
1301 """Don't record coordinates for the given domain.
1302 Coordinates for the given domain (specified by a component name
1303 and a 2-element tuple giving the start and end residue numbers)
1304 will be excluded from the mmCIF file. This can be used to exclude
1305 parts of the structure that weren't well resolved in modeling.
1306 Any bead or residue that lies wholly within this range will be
1307 excluded. Multiple ranges for a given component can be excluded
1308 by calling this method multiple times."""
1309 if component
not in self._exclude_coords:
1310 self._exclude_coords[component] = []
1311 self._exclude_coords[component].append(seqrange)
1313 def _is_excluded(self, component, start, end):
1314 """Return True iff this chunk of sequence should be excluded"""
1315 for seqrange
in self._exclude_coords.get(component, ()):
1316 if start >= seqrange[0]
and end <= seqrange[1]:
1319 def _add_state(self, state):
1320 """Create a new state and return a pointer to it."""
1321 self._state_ensemble_offset = len(self.system.ensembles)
1322 s = _State(state, self)
1323 self._state_group.append(s)
1324 self._last_state = s
1327 def _get_chain_for_component(self, name, output):
1328 """Get the chain ID for a component, if any."""
1330 if name
in self.asym_units:
1331 return self.asym_units[name]._id
1336 def _get_assembly_comps(self, assembly):
1337 """Get the names of the components in the given assembly"""
1341 comps[ca.details] =
None
1345 """Make a new component that's a transformed copy of another.
1346 All representation for the existing component is copied to the
1348 assembly_comps = self._get_assembly_comps(state.modeled_assembly)
1349 if name
in assembly_comps:
1350 raise ValueError(
"Component %s already exists" % name)
1351 elif original
not in assembly_comps:
1352 raise ValueError(
"Original component %s does not exist" % original)
1353 self.create_component(state, name,
True)
1354 self.add_component_sequence(state, name, self.sequence_dict[original])
1355 self._transformed_components.append(_TransformedComponent(
1356 name, original, transform))
1357 self.all_representations.copy_component(state, name, original,
1358 self.asym_units[name])
1360 def create_component(self, state, name, modeled, asym_name=None):
1361 if asym_name
is None:
1363 new_comp = name
not in self._all_components
1364 self._all_components[name] =
None
1366 state.all_modeled_components.append(name)
1367 if asym_name
not in self.asym_units:
1369 self.asym_units[asym_name] =
None
1371 self.all_modeled_components.append(name)
1373 def add_component_sequence(self, state, name, seq, asym_name=None,
1375 if asym_name
is None:
1378 def get_offset(seq):
1380 for i
in range(len(seq)):
1383 seq, offset = get_offset(seq)
1384 if name
in self.sequence_dict:
1385 if self.sequence_dict[name] != seq:
1386 raise ValueError(
"Sequence mismatch for component %s" % name)
1388 self.sequence_dict[name] = seq
1389 self.entities.add(name, seq, offset, alphabet)
1390 if asym_name
in self.asym_units:
1391 if self.asym_units[asym_name]
is None:
1393 entity = self.entities[name]
1394 asym =
AsymUnit(entity, details=asym_name)
1395 self.system.asym_units.append(asym)
1396 self.asym_units[asym_name] = asym
1397 state.modeled_assembly.append(self.asym_units[asym_name])
1399 def _add_restraint_model_fits(self):
1400 """Add fits to restraints for all known models"""
1401 for group, m
in self.system._all_models():
1402 if m._is_restrained:
1403 for r
in self.system.restraints:
1404 if hasattr(r,
'add_fits_from_model_statfile'):
1405 r.add_fits_from_model_statfile(m)
1408 """Do any final processing on the class hierarchy.
1409 After calling this method, the `system` member (an instance
1410 of `ihm.System`) completely reproduces the PMI modeling, and
1411 can be written out to an mmCIF file with `ihm.dumper.write`,
1412 and/or modified using the ihm API."""
1413 self._add_restraint_model_fits()
1415 def add_pdb_element(self, state, name, start, end, offset, pdbname,
1416 chain, hier, representation=
None):
1417 if self._is_excluded(name, start, end):
1419 if representation
is None:
1420 representation = self.default_representation
1421 asym = self.asym_units[name]
1422 p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1424 self.all_representations.add_fragment(state, representation, p)
1425 self.all_starting_models.add_pdb_fragment(p)
1427 def add_bead_element(self, state, name, start, end, num, hier,
1428 representation=
None):
1429 if self._is_excluded(name, start, end):
1431 if representation
is None:
1432 representation = self.default_representation
1433 asym = self.asym_units[name]
1434 pmi_offset = asym.entity.pmi_offset
1435 b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1437 self.all_representations.add_fragment(state, representation, b)
1439 def get_cross_link_group(self, pmi_restraint):
1440 r = _CrossLinkRestraint(pmi_restraint)
1441 self.system.restraints.append(r)
1442 self._add_restraint_dataset(r)
1445 def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1446 if c1
not in self._all_components
or c2
not in self._all_components:
1452 e1 = self.entities[c1]
1453 e2 = self.entities[c2]
1454 xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1455 residue2=e2.pmi_residue(r2))
1456 rsr.experimental_cross_links.append([xl])
1459 def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1462 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1463 d = ihm.restraint.UpperBoundDistanceRestraint(length)
1465 if _get_by_residue(p1)
and _get_by_residue(p2):
1466 cls = _ResidueCrossLink
1468 cls = _FeatureCrossLink
1469 xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1472 xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1473 rsr.cross_links.append(xl)
1475 def add_replica_exchange(self, state, rex):
1480 step = _ReplicaExchangeProtocolStep(state, rex)
1481 step.software = self.all_software.pmi
1482 self.all_protocols.add_step(step, state)
1484 def _add_simple_dynamics(self, num_models_end, method):
1486 state = self._last_state
1487 self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1490 def _add_protocol(self):
1492 state = self._last_state
1493 self.all_protocols.add_protocol(state)
1495 def _add_dataset(self, dataset):
1496 return self.all_datasets.add(self._last_state, dataset)
1498 def _add_restraint_dataset(self, restraint):
1499 return self.all_datasets.add_restraint(self._last_state, restraint)
1501 def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1503 state = self._last_state
1504 pp = ihm.analysis.ClusterStep(
'RMSD', num_models_begin, num_models_end)
1505 self.all_protocols.add_postproc(pp, state)
1508 def _add_no_postprocessing(self, num_models):
1510 state = self._last_state
1511 pp = ihm.analysis.EmptyStep()
1512 pp.num_models_begin = pp.num_models_end = num_models
1513 self.all_protocols.add_postproc(pp, state)
1516 def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1517 num_models_deposited, localization_densities,
1519 """Add an ensemble generated by ad hoc methods (not using PMI).
1520 This is currently only used by the Nup84 system."""
1522 state = self._last_state
1523 group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1524 state.add_model_group(group)
1526 self.system.locations.append(ensemble_file)
1527 e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1529 self.system.ensembles.append(e)
1530 for c
in state.all_modeled_components:
1531 den = localization_densities.get(c,
None)
1533 e.load_localization_density(state, c, self.asym_units[c], den)
1537 """Point a previously-created ensemble to an 'all-models' file.
1538 This could be a trajectory such as DCD, an RMF, or a multimodel
1540 self.system.locations.append(location)
1542 ind = i + self._state_ensemble_offset
1543 self.system.ensembles[ind].file = location
1545 def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1551 protocol = self.all_protocols.get_last_protocol(state)
1552 num_models = protocol.steps[-1].num_models_end
1553 pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1554 pp.software = self.all_software.pmi
1555 self.all_protocols.add_postproc(pp, state)
1556 for i
in range(rex._number_of_clusters):
1557 group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1558 'cluster %d' % (i + 1)))
1559 state.add_model_group(group)
1561 e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1562 self.system.ensembles.append(e)
1564 for fname, stuple
in sorted(density_custom_ranges.items()):
1565 e.load_localization_density(state, fname, stuple,
1567 for stats
in e.load_all_models(self, state):
1568 m = self.add_model(group)
1571 m.name =
'Best scoring model'
1574 for c
in state.all_modeled_components:
1577 def _get_subassembly(self, comps, name, description):
1578 """Get an Assembly consisting of the given components.
1579 `compdict` is a dictionary of the components to add, where keys
1580 are the component names and values are the sequence ranges (or
1581 None to use all residues in the component)."""
1583 for comp, seqrng
in comps.items():
1584 a = self.asym_units[comp]
1585 asyms.append(a
if seqrng
is None else a(*seqrng))
1587 a = ihm.Assembly(asyms, name=name, description=description)
1590 def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1592 """Add a basic FoXS fit. This is largely intended for use from the
1594 assembly = self._get_subassembly(
1596 name=
"SAXS subassembly",
1597 description=
"All components that fit SAXS data")
1598 r = ihm.restraint.SASRestraint(
1599 dataset, assembly, segment=
False,
1600 fitting_method=
'FoXS', fitting_atom_type=
'Heavy atoms',
1601 multi_state=
False, radius_of_gyration=rg, details=details)
1602 r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1603 self.system.restraints.append(r)
1604 self._add_restraint_dataset(r)
1606 def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1607 image_resolution, projection_number,
1608 micrographs_number):
1609 r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1610 image_resolution, projection_number,
1612 self.system.restraints.append(r)
1613 self._add_restraint_dataset(r)
1615 def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1617 r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1618 self.system.restraints.append(r)
1619 self._add_restraint_dataset(r)
1621 def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1622 sigma, pmi_restraint):
1623 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1624 sigma, pmi_restraint, self._xy_plane)
1626 def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1627 sigma, pmi_restraint):
1628 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1629 sigma, pmi_restraint, self._xz_plane)
1631 def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1632 sigma, pmi_restraint):
1633 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1634 sigma, pmi_restraint, self._z_axis)
1636 def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1637 sigma, pmi_restraint, geom):
1638 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1639 r = _GeometricRestraint(
1640 self, state, pmi_restraint, geom, asym.get_feature(ps),
1641 ihm.restraint.LowerUpperBoundDistanceRestraint(lower_bound,
1644 self.system.restraints.append(r)
1645 self._add_restraint_dataset(r)
1647 def _get_membrane(self, tor_R, tor_r, tor_th):
1648 """Get an object representing a half-torus membrane"""
1649 if not hasattr(self,
'_seen_membranes'):
1650 self._seen_membranes = {}
1653 membrane_id = tuple(int(x * 100.)
for x
in (tor_R, tor_r, tor_th))
1654 if membrane_id
not in self._seen_membranes:
1655 m = ihm.geometry.HalfTorus(
1656 center=self._center_origin,
1657 transformation=self._identity_transform,
1658 major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1659 inner=
True, name=
'Membrane')
1660 self._seen_membranes[membrane_id] = m
1661 return self._seen_membranes[membrane_id]
1663 def add_membrane_surface_location_restraint(
1664 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1665 self._add_membrane_restraint(
1666 state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint,
1667 ihm.restraint.UpperBoundDistanceRestraint(0.))
1669 def add_membrane_exclusion_restraint(
1670 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1671 self._add_membrane_restraint(
1672 state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint,
1673 ihm.restraint.LowerBoundDistanceRestraint(0.))
1675 def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1676 sigma, pmi_restraint, rsr):
1677 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1678 r = _GeometricRestraint(
1679 self, state, pmi_restraint,
1680 self._get_membrane(tor_R, tor_r, tor_th), asym.get_feature(ps),
1682 self.system.restraints.append(r)
1683 self._add_restraint_dataset(r)
1685 def add_model(self, group, assembly=None, representation=None):
1686 state = self._last_state
1687 if representation
is None:
1688 representation = self.default_representation
1689 protocol = self.all_protocols.get_last_protocol(state)
1690 m = _Model(state.prot, self, protocol,
1691 assembly
if assembly
else state.modeled_assembly,
1698 """Get custom python-ihm dumpers for writing PMI to from mmCIF.
1699 This returns a list of custom dumpers that can be passed as all or
1700 part of the `dumpers` argument to ihm.dumper.write(). They add
1701 PMI-specific information to mmCIF or BinaryCIF files written out
1703 return [_ReplicaExchangeProtocolDumper]
1707 """Get custom python-ihm handlers for reading PMI data from mmCIF.
1708 This returns a list of custom handlers that can be passed as all or
1709 part of the `handlers` argument to ihm.reader.read(). They read
1710 PMI-specific information from mmCIF or BinaryCIF files read in
1712 return [_ReplicaExchangeProtocolHandler]
1716 """Extract metadata from an EM density GMM file."""
1719 """Extract metadata from `filename`.
1720 @return a dict with key `dataset` pointing to the GMM file and
1721 `number_of_gaussians` to the number of GMMs (or None)"""
1722 loc = ihm.location.InputFileLocation(
1724 details=
"Electron microscopy density map, "
1725 "represented as a Gaussian Mixture Model (GMM)")
1728 loc._allow_duplicates =
True
1729 d = ihm.dataset.EMDensityDataset(loc)
1730 ret = {
'dataset': d,
'number_of_gaussians':
None}
1732 with open(filename)
as fh:
1734 if line.startswith(
'# data_fn: '):
1735 p = ihm.metadata.MRCParser()
1736 fn = line[11:].rstrip(
'\r\n')
1737 dataset = p.parse_file(os.path.join(
1738 os.path.dirname(filename), fn))[
'dataset']
1739 ret[
'dataset'].parents.append(dataset)
1740 elif line.startswith(
'# ncenters: '):
1741 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 ...