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
30 from pathlib
import Path
32 from IMP._compat_pathlib
import Path
43 import ihm.representation
45 import ihm.cross_linkers
48 def _assign_id(obj, seen_objs, obj_by_id):
49 """Assign a unique ID to obj, and track all ids in obj_by_id."""
50 if obj
not in seen_objs:
51 if not hasattr(obj,
'id'):
53 obj.id = len(obj_by_id)
54 seen_objs[obj] = obj.id
56 obj.id = seen_objs[obj]
59 def _get_by_residue(p):
60 """Determine whether the given particle represents a specific residue
61 or a more coarse-grained object."""
65 class _ComponentMapper(object):
66 """Map a Particle to a component name"""
67 def __init__(self, prot):
70 self.name =
'cif-output'
71 self.o.dictionary_pdbs[self.name] = self.prot
72 self.o._init_dictchain(self.name, self.prot,
73 multichar_chain=
True, mmcif=
True)
75 def __getitem__(self, p):
76 protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
80 class _AsymMapper(object):
81 """Map a Particle to an asym_unit"""
82 def __init__(self, simo, prot):
84 self._cm = _ComponentMapper(prot)
85 self._seen_ranges = {}
87 def __getitem__(self, p):
88 protname = self._cm[p]
89 return self.simo.asym_units[protname]
91 def get_feature(self, ps):
92 """Get an ihm.restraint.Feature that covers the given particles"""
100 rng = asym(rind, rind)
104 rng = asym(rinds[0], rinds[-1])
106 raise ValueError(
"Unsupported particle type %s" % str(p))
108 if len(rngs) > 0
and rngs[-1].asym == asym \
109 and rngs[-1].seq_id_range[1] == rng.seq_id_range[0] - 1:
110 rngs[-1].seq_id_range = (rngs[-1].seq_id_range[0],
117 if hrngs
in self._seen_ranges:
118 return self._seen_ranges[hrngs]
120 feat = ihm.restraint.ResidueFeature(rngs)
121 self._seen_ranges[hrngs] = feat
125 class _AllSoftware(object):
126 def __init__(self, system):
128 self.modeller_used = self.phyre2_used =
False
129 self.pmi = ihm.Software(
130 name=
"IMP PMI module",
131 version=IMP.pmi.__version__,
132 classification=
"integrative model building",
133 description=
"integrative model building",
134 location=
'https://integrativemodeling.org')
135 self.imp = ihm.Software(
136 name=
"Integrative Modeling Platform (IMP)",
137 version=IMP.__version__,
138 classification=
"integrative model building",
139 description=
"integrative model building",
140 location=
'https://integrativemodeling.org')
143 if hasattr(self.imp,
'citation'):
144 if sys.version_info[0] > 2:
146 javi =
'Vel\u00e1zquez-Muriel J'
148 javi =
'Velazquez-Muriel J'
149 self.imp.citation = ihm.Citation(
151 title=
'Putting the pieces together: integrative modeling '
152 'platform software for structure determination of '
153 'macromolecular assemblies',
154 journal=
'PLoS Biol', volume=10, page_range=
'e1001244',
156 authors=[
'Russel D',
'Lasker K',
'Webb B', javi,
'Tjioe E',
157 'Schneidman-Duhovny D',
'Peterson B',
'Sali A'],
158 doi=
'10.1371/journal.pbio.1001244')
159 self.pmi.citation = ihm.Citation(
161 title=
'Modeling Biological Complexes Using Integrative '
162 'Modeling Platform.',
163 journal=
'Methods Mol Biol', volume=2022, page_range=(353, 377),
165 authors=[
'Saltzberg D',
'Greenberg CH',
'Viswanath S',
166 'Chemmama I',
'Webb B',
'Pellarin R',
'Echeverria I',
168 doi=
'10.1007/978-1-4939-9608-7_15')
169 self.system.software.extend([self.pmi, self.imp])
171 def set_modeller_used(self, version, date):
172 if self.modeller_used:
174 self.modeller_used =
True
176 name=
'MODELLER', classification=
'comparative modeling',
177 description=
'Comparative modeling by satisfaction '
178 'of spatial restraints, build ' + date,
179 location=
'https://salilab.org/modeller/', version=version)
180 self.system.software.append(s)
181 if hasattr(s,
'citation'):
182 s.citation = ihm.Citation(
184 title=
'Comparative protein modelling by satisfaction of '
185 'spatial restraints.',
186 journal=
'J Mol Biol', volume=234, page_range=(779, 815),
187 year=1993, authors=[
'Sali A',
'Blundell TL'],
188 doi=
'10.1006/jmbi.1993.1626')
190 def set_phyre2_used(self):
193 self.phyre2_used =
True
195 name=
'Phyre2', classification=
'protein homology modeling',
196 description=
'Protein Homology/analogY Recognition Engine V 2.0',
197 version=
'2.0', location=
'http://www.sbg.bio.ic.ac.uk/~phyre2/')
198 if hasattr(s,
'citation'):
199 s.citation = ihm.Citation(
201 title=
'The Phyre2 web portal for protein modeling, '
202 'prediction and analysis.',
203 journal=
'Nat Protoc', volume=10, page_range=(845, 858),
204 authors=[
'Kelley LA',
'Mezulis S',
'Yates CM',
'Wass MN',
206 year=2015, doi=
'10.1038/nprot.2015.053')
207 self.system.software.append(s)
210 def _get_fragment_is_rigid(fragment):
211 """Determine whether a fragment is modeled rigidly"""
217 class _PDBFragment(ihm.representation.ResidueSegment):
218 """Record details about part of a PDB file used as input
220 def __init__(self, state, component, start, end, pdb_offset,
221 pdbname, chain, hier, asym_unit):
223 super(_PDBFragment, self).__init__(
224 asym_unit=asym_unit.pmi_range(start, end),
225 rigid=
None, primitive=
'sphere')
226 self.component, self.start, self.end, self.offset, self.pdbname \
227 = component, start, end, pdb_offset, pdbname
228 self.state, self.chain, self.hier = state, chain, hier
233 rigid = property(
lambda self: _get_fragment_is_rigid(self),
234 lambda self, val:
None)
236 def combine(self, other):
240 class _BeadsFragment(ihm.representation.FeatureSegment):
241 """Record details about beads used to represent part of a component."""
244 def __init__(self, state, component, start, end, count, hier, asym_unit):
245 super(_BeadsFragment, self).__init__(
246 asym_unit=asym_unit(start, end), rigid=
None, primitive=
'sphere',
248 self.state, self.component, self.hier = state, component, hier
250 rigid = property(
lambda self: _get_fragment_is_rigid(self),
251 lambda self, val:
None)
253 def combine(self, other):
255 if (type(other) == type(self)
and
256 other.asym_unit.seq_id_range[0]
257 == self.asym_unit.seq_id_range[1] + 1):
258 self.asym_unit.seq_id_range = (self.asym_unit.seq_id_range[0],
259 other.asym_unit.seq_id_range[1])
260 self.count += other.count
264 class _AllModelRepresentations(object):
265 def __init__(self, simo):
269 self.fragments = OrderedDict()
270 self._all_representations = {}
272 def copy_component(self, state, name, original, asym_unit):
273 """Copy all representation for `original` in `state` to `name`"""
276 newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
278 for rep
in self.fragments:
279 if original
in self.fragments[rep]:
280 if name
not in self.fragments[rep]:
281 self.fragments[rep][name] = OrderedDict()
282 self.fragments[rep][name][state] = [
283 copy_frag(f)
for f
in self.fragments[rep][original][state]]
286 first_state = list(self.fragments[rep][name].keys())[0]
287 if state
is first_state:
288 representation = self._all_representations[rep]
289 representation.extend(self.fragments[rep][name][state])
291 def add_fragment(self, state, representation, fragment):
292 """Add a model fragment."""
293 comp = fragment.component
294 id_rep = id(representation)
295 self._all_representations[id_rep] = representation
296 if id_rep
not in self.fragments:
297 self.fragments[id_rep] = OrderedDict()
298 if comp
not in self.fragments[id_rep]:
299 self.fragments[id_rep][comp] = OrderedDict()
300 if state
not in self.fragments[id_rep][comp]:
301 self.fragments[id_rep][comp][state] = []
302 fragments = self.fragments[id_rep][comp][state]
303 if len(fragments) == 0
or not fragments[-1].combine(fragment):
304 fragments.append(fragment)
307 first_state = list(self.fragments[id_rep][comp].keys())[0]
308 if state
is first_state:
309 representation.append(fragment)
312 class _AllDatasets(object):
313 """Track all datasets generated by PMI and add them to the ihm.System"""
314 def __init__(self, system):
317 self._datasets_by_state = {}
318 self._restraints_by_state = {}
320 def get_all_group(self, state):
321 """Get a DatasetGroup encompassing all datasets so far in this state"""
325 g = ihm.dataset.DatasetGroup(
326 self._datasets_by_state.get(state, [])
327 + [r.dataset
for r
in self._restraints_by_state.get(state, [])
331 def add(self, state, dataset):
332 """Add a new dataset."""
333 self._datasets.append(dataset)
334 if state
not in self._datasets_by_state:
335 self._datasets_by_state[state] = []
336 self._datasets_by_state[state].append(dataset)
339 self.system.orphan_datasets.append(dataset)
343 """Add the dataset for a restraint"""
344 if state
not in self._restraints_by_state:
345 self._restraints_by_state[state] = []
346 self._restraints_by_state[state].append(restraint)
349 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
350 """Restrain to a set of cross-links"""
353 _label_map = {
'wtDSS':
'DSS',
'scDSS':
'DSS',
'scEDC':
'EDC'}
354 _descriptor_map = {
'DSS': ihm.cross_linkers.dss,
355 'EDC': ihm.cross_linkers.edc}
357 def __init__(self, pmi_restraint):
358 self.pmi_restraint = pmi_restraint
361 linker = getattr(self.pmi_restraint,
'linker',
None)
362 label = self.pmi_restraint.label
364 super(_CrossLinkRestraint, self).__init__(
365 dataset=self.pmi_restraint.dataset,
366 linker=linker
or self._get_chem_descriptor(label))
369 def _get_chem_descriptor(cls, label):
371 label = cls._label_map.get(label, label)
372 if label
not in cls._descriptor_map:
376 d = ihm.ChemDescriptor(label)
377 cls._descriptor_map[label] = d
378 return cls._descriptor_map[label]
380 def _set_psi_sigma(self, model):
383 if model.m != self.pmi_restraint.model:
385 for resolution
in self.pmi_restraint.sigma_dictionary:
386 statname =
'ISDCrossLinkMS_Sigma_%s_%s' % (resolution, self.label)
387 if model.stats
and statname
in model.stats:
388 sigma = float(model.stats[statname])
389 p = self.pmi_restraint.sigma_dictionary[resolution][0]
390 old_values.append((p, p.get_scale()))
392 for psiindex
in self.pmi_restraint.psi_dictionary:
393 statname =
'ISDCrossLinkMS_Psi_%s_%s' % (psiindex, self.label)
394 if model.stats
and statname
in model.stats:
395 psi = float(model.stats[statname])
396 p = self.pmi_restraint.psi_dictionary[psiindex][0]
397 old_values.append((p, p.get_scale()))
401 return list(reversed(old_values))
403 def add_fits_from_model_statfile(self, model):
405 old_values = self._set_psi_sigma(model)
409 for xl
in self.cross_links:
411 xl.fits[model] = ihm.restraint.CrossLinkFit(
412 psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
414 for p, scale
in old_values:
418 def __set_dataset(self, val):
419 self.pmi_restraint.dataset = val
420 dataset = property(
lambda self: self.pmi_restraint.dataset,
424 def get_asym_mapper_for_state(simo, state, asym_map):
425 asym = asym_map.get(state,
None)
427 asym = _AsymMapper(simo, state.prot)
428 asym_map[state] = asym
432 class _PMICrossLink(object):
435 psi = property(
lambda self: self.psi_p.get_scale(),
436 lambda self, val:
None)
437 sigma1 = property(
lambda self: self.sigma1_p.get_scale(),
438 lambda self, val:
None)
439 sigma2 = property(
lambda self: self.sigma2_p.get_scale(),
440 lambda self, val:
None)
443 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
447 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
451 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
452 def __init__(self, state, pmi_restraint, image_number, resolution,
453 pixel_size, image_resolution, projection_number,
455 self.pmi_restraint, self.image_number = pmi_restraint, image_number
456 super(_EM2DRestraint, self).__init__(
457 dataset=pmi_restraint.datasets[image_number],
458 assembly=state.modeled_assembly,
459 segment=
False, number_raw_micrographs=micrographs_number,
460 pixel_size_width=pixel_size, pixel_size_height=pixel_size,
461 image_resolution=image_resolution,
462 number_of_projections=projection_number)
465 def __get_dataset(self):
466 return self.pmi_restraint.datasets[self.image_number]
468 def __set_dataset(self, val):
469 self.pmi_restraint.datasets[self.image_number] = val
471 dataset = property(__get_dataset, __set_dataset)
473 def add_fits_from_model_statfile(self, model):
474 ccc = self._get_cross_correlation(model)
475 transform = self._get_transformation(model)
476 rot = transform.get_rotation()
477 rm = [[e
for e
in rot.get_rotation_matrix_row(i)]
for i
in range(3)]
478 self.fits[model] = ihm.restraint.EM2DRestraintFit(
479 cross_correlation_coefficient=ccc,
481 tr_vector=transform.get_translation())
483 def _get_transformation(self, model):
484 """Get the transformation that places the model on the image"""
485 stats = model.em2d_stats
or model.stats
486 prefix =
'ElectronMicroscopy2D_%s_Image%d' % (self.pmi_restraint.label,
487 self.image_number + 1)
488 r = [float(stats[prefix +
'_Rotation%d' % i])
for i
in range(4)]
489 t = [float(stats[prefix +
'_Translation%d' % i])
493 inv = model.transform.get_inverse()
497 def _get_cross_correlation(self, model):
498 """Get the cross correlation coefficient between the model projection
500 stats = model.em2d_stats
or model.stats
501 return float(stats[
'ElectronMicroscopy2D_%s_Image%d_CCC'
502 % (self.pmi_restraint.label,
503 self.image_number + 1)])
506 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
508 def __init__(self, simo, state, pmi_restraint, target_ps, densities):
509 self.pmi_restraint = pmi_restraint
510 super(_EM3DRestraint, self).__init__(
511 dataset=pmi_restraint.dataset,
512 assembly=self._get_assembly(densities, simo, state),
513 fitting_method=
'Gaussian mixture models',
514 number_of_gaussians=len(target_ps))
517 def __set_dataset(self, val):
518 self.pmi_restraint.dataset = val
519 dataset = property(
lambda self: self.pmi_restraint.dataset,
522 def _get_assembly(self, densities, simo, state):
523 """Get the Assembly that this restraint acts on"""
524 cm = _ComponentMapper(state.prot)
527 components[cm[d]] =
None
528 a = simo._get_subassembly(
529 components, name=
"EM subassembly",
530 description=
"All components that fit the EM map")
533 def add_fits_from_model_statfile(self, model):
534 ccc = self._get_cross_correlation(model)
535 self.fits[model] = ihm.restraint.EM3DRestraintFit(
536 cross_correlation_coefficient=ccc)
538 def _get_cross_correlation(self, model):
539 """Get the cross correlation coefficient between the model
541 if model.stats
is not None:
542 return float(model.stats[
'GaussianEMRestraint_%s_CCC'
543 % self.pmi_restraint.label])
546 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
548 def __init__(self, simo, state, pmi_restraint, geometric_object,
549 feature, distance, sigma):
550 self.pmi_restraint = pmi_restraint
551 super(_GeometricRestraint, self).__init__(
552 dataset=pmi_restraint.dataset,
553 geometric_object=geometric_object, feature=feature,
554 distance=distance, harmonic_force_constant=1. / sigma,
558 def __set_dataset(self, val):
559 self.pmi_restraint.dataset = val
560 dataset = property(
lambda self: self.pmi_restraint.dataset,
564 class _ReplicaExchangeProtocolStep(ihm.protocol.Step):
565 def __init__(self, state, rex):
566 if rex.monte_carlo_sample_objects
is not None:
567 method =
'Replica exchange monte carlo'
569 method =
'Replica exchange molecular dynamics'
570 self.monte_carlo_temperature = rex.vars[
'monte_carlo_temperature']
571 self.replica_exchange_minimum_temperature = \
572 rex.vars[
'replica_exchange_minimum_temperature']
573 self.replica_exchange_maximum_temperature = \
574 rex.vars[
'replica_exchange_maximum_temperature']
575 super(_ReplicaExchangeProtocolStep, self).__init__(
576 assembly=state.modeled_assembly,
578 method=method, name=
'Sampling',
579 num_models_begin=
None,
580 num_models_end=rex.vars[
"number_of_frames"],
581 multi_scale=
True, multi_state=
False, ordered=
False, ensemble=
True)
584 class _ReplicaExchangeProtocolDumper(ihm.dumper.Dumper):
585 """Write IMP-specific information about replica exchange to mmCIF.
586 Note that IDs will have already been assigned by python-ihm's
587 standard modeling protocol dumper."""
588 def dump(self, system, writer):
589 with writer.loop(
"_imp_replica_exchange_protocol",
590 [
"protocol_id",
"step_id",
"monte_carlo_temperature",
591 "replica_exchange_minimum_temperature",
592 "replica_exchange_maximum_temperature"])
as lp:
593 for p
in system._all_protocols():
595 if isinstance(s, _ReplicaExchangeProtocolStep):
596 self._dump_step(p, s, lp)
598 def _dump_step(self, p, s, lp):
599 mintemp = s.replica_exchange_minimum_temperature
600 maxtemp = s.replica_exchange_maximum_temperature
601 lp.write(protocol_id=p._id, step_id=s._id,
602 monte_carlo_temperature=s.monte_carlo_temperature,
603 replica_exchange_minimum_temperature=mintemp,
604 replica_exchange_maximum_temperature=maxtemp)
607 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
608 category =
'_imp_replica_exchange_protocol'
610 """Read IMP-specific information about replica exchange from mmCIF."""
611 def __call__(self, protocol_id, step_id, monte_carlo_temperature,
612 replica_exchange_minimum_temperature,
613 replica_exchange_maximum_temperature):
614 p = self.sysr.protocols.get_by_id(protocol_id)
616 s = p.steps[int(step_id)-1]
618 s.__class__ = _ReplicaExchangeProtocolStep
619 s.monte_carlo_temperature = \
620 self.get_float(monte_carlo_temperature)
621 s.replica_exchange_minimum_temperature = \
622 self.get_float(replica_exchange_minimum_temperature)
623 s.replica_exchange_maximum_temperature = \
624 self.get_float(replica_exchange_maximum_temperature)
627 class _SimpleProtocolStep(ihm.protocol.Step):
628 def __init__(self, state, num_models_end, method):
629 super(_SimpleProtocolStep, self).__init__(
630 assembly=state.modeled_assembly,
632 method=method, name=
'Sampling',
633 num_models_begin=
None,
634 num_models_end=num_models_end,
635 multi_scale=
True, multi_state=
False, ordered=
False,
639 class _Chain(object):
640 """Represent a single chain in a Model"""
641 def __init__(self, pmi_chain_id, asym_unit):
642 self.pmi_chain_id, self.asym_unit = pmi_chain_id, asym_unit
646 def add(self, xyz, atom_type, residue_type, residue_index,
647 all_indexes, radius):
648 if atom_type
is None:
649 self.spheres.append((xyz, residue_type, residue_index,
650 all_indexes, radius))
652 self.atoms.append((xyz, atom_type, residue_type, residue_index,
653 all_indexes, radius))
654 orig_comp = property(
lambda self: self.comp)
657 class _TransformedChain(object):
658 """Represent a chain that is a transformed version of another"""
659 def __init__(self, orig_chain, asym_unit, transform):
660 self.orig_chain, self.asym_unit = orig_chain, asym_unit
661 self.transform = transform
663 def __get_spheres(self):
664 for (xyz, residue_type, residue_index, all_indexes,
665 radius)
in self.orig_chain.spheres:
666 yield (self.transform * xyz, residue_type, residue_index,
668 spheres = property(__get_spheres)
670 def __get_atoms(self):
671 for (xyz, atom_type, residue_type, residue_index, all_indexes,
672 radius)
in self.orig_chain.atoms:
673 yield (self.transform * xyz, atom_type, residue_type,
674 residue_index, all_indexes, radius)
675 atoms = property(__get_atoms)
677 entity = property(
lambda self: self.orig_chain.entity)
678 orig_comp = property(
lambda self: self.orig_chain.comp)
681 class _Excluder(object):
682 def __init__(self, component, simo):
683 self._seqranges = simo._exclude_coords.get(component, [])
685 def is_excluded(self, indexes):
686 """Return True iff the given sequence range is excluded."""
687 for seqrange
in self._seqranges:
688 if indexes[0] >= seqrange[0]
and indexes[-1] <= seqrange[1]:
692 class _Model(ihm.model.Model):
693 def __init__(self, prot, simo, protocol, assembly, representation):
694 super(_Model, self).__init__(assembly=assembly, protocol=protocol,
695 representation=representation)
696 self.simo = weakref.proxy(simo)
702 self.em2d_stats =
None
705 self._is_restrained =
True
708 self.m = prot.get_model()
709 o.dictionary_pdbs[name] = prot
710 o._init_dictchain(name, prot, multichar_chain=
True)
711 (particle_infos_for_pdb,
712 self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
714 self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
717 def all_chains(self, simo):
718 """Yield all chains, including transformed ones"""
720 for c
in self.chains:
722 chain_for_comp[c.comp] = c
723 for tc
in simo._transformed_components:
724 orig_chain = chain_for_comp.get(tc.original,
None)
726 asym = simo.asym_units[tc.name]
727 c = _TransformedChain(orig_chain, asym, tc.transform)
731 def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
732 entity_for_chain = {}
735 for protname, chain_id
in o.dictchain[name].items():
736 if protname
in simo.entities:
737 entity_for_chain[chain_id] = simo.entities[protname]
740 pn = protname.split(
'.')[0]
741 entity_for_chain[chain_id] = simo.entities[pn]
742 comp_for_chain[chain_id] = protname
746 correct_asym[chain_id] = simo.asym_units[protname]
753 for (xyz, atom_type, residue_type, chain_id, residue_index,
754 all_indexes, radius)
in particle_infos_for_pdb:
755 if chain
is None or chain.pmi_chain_id != chain_id:
756 chain = _Chain(chain_id, correct_asym[chain_id])
757 chain.entity = entity_for_chain[chain_id]
758 chain.comp = comp_for_chain[chain_id]
759 self.chains.append(chain)
760 excluder = _Excluder(chain.comp, simo)
761 if not excluder.is_excluded(all_indexes
if all_indexes
762 else [residue_index]):
763 chain.add(xyz, atom_type, residue_type, residue_index,
766 def parse_rmsf_file(self, fname, component):
767 self.rmsf[component] = rmsf = {}
768 with open(str(fname))
as fh:
770 resnum, blocknum, val = line.split()
771 rmsf[int(resnum)] = (int(blocknum), float(val))
773 def get_rmsf(self, component, indexes):
774 """Get the RMSF value for the given residue indexes."""
777 rmsf = self.rmsf[component]
778 blocknums = dict.fromkeys(rmsf[ind][0]
for ind
in indexes)
779 if len(blocknums) != 1:
780 raise ValueError(
"Residue indexes %s aren't all in the same block"
782 return rmsf[indexes[0]][1]
785 for chain
in self.all_chains(self.simo):
786 pmi_offset = chain.asym_unit.entity.pmi_offset
787 for atom
in chain.atoms:
788 (xyz, atom_type, residue_type, residue_index,
789 all_indexes, radius) = atom
790 pt = self.transform * xyz
791 yield ihm.model.Atom(
792 asym_unit=chain.asym_unit,
793 seq_id=residue_index - pmi_offset,
794 atom_id=atom_type.get_string(),
796 x=pt[0], y=pt[1], z=pt[2])
798 def get_spheres(self):
799 for chain
in self.all_chains(self.simo):
800 pmi_offset = chain.asym_unit.entity.pmi_offset
801 for sphere
in chain.spheres:
802 (xyz, residue_type, residue_index,
803 all_indexes, radius) = sphere
804 if all_indexes
is None:
805 all_indexes = (residue_index,)
806 pt = self.transform * xyz
807 yield ihm.model.Sphere(
808 asym_unit=chain.asym_unit,
809 seq_id_range=(all_indexes[0] - pmi_offset,
810 all_indexes[-1] - pmi_offset),
811 x=pt[0], y=pt[1], z=pt[2], radius=radius,
812 rmsf=self.get_rmsf(chain.orig_comp, all_indexes))
815 class _AllProtocols(object):
816 def __init__(self, simo):
817 self.simo = weakref.proxy(simo)
819 self.protocols = OrderedDict()
821 def add_protocol(self, state):
822 """Add a new Protocol"""
823 if state
not in self.protocols:
824 self.protocols[state] = []
825 p = ihm.protocol.Protocol()
826 self.simo.system.orphan_protocols.append(p)
827 self.protocols[state].append(p)
829 def add_step(self, step, state):
830 """Add a ProtocolStep to the last Protocol of the given State"""
831 if state
not in self.protocols:
832 self.add_protocol(state)
833 protocol = self.get_last_protocol(state)
834 if len(protocol.steps) == 0:
835 step.num_models_begin = 0
837 step.num_models_begin = protocol.steps[-1].num_models_end
838 protocol.steps.append(step)
839 step.id = len(protocol.steps)
841 step.dataset_group = self.simo.all_datasets.get_all_group(state)
843 def add_postproc(self, step, state):
844 """Add a postprocessing step to the last protocol"""
845 protocol = self.get_last_protocol(state)
846 if len(protocol.analyses) == 0:
847 protocol.analyses.append(ihm.analysis.Analysis())
848 protocol.analyses[-1].steps.append(step)
850 def get_last_protocol(self, state):
851 """Return the most recently-added _Protocol"""
852 return self.protocols[state][-1]
855 class _AllStartingModels(object):
856 def __init__(self, simo):
860 self.models = OrderedDict()
863 def add_pdb_fragment(self, fragment):
864 """Add a starting model PDB fragment."""
865 comp = fragment.component
866 state = fragment.state
867 if comp
not in self.models:
868 self.models[comp] = OrderedDict()
869 if state
not in self.models[comp]:
870 self.models[comp][state] = []
871 models = self.models[comp][state]
872 if len(models) == 0 \
873 or models[-1].fragments[0].pdbname != fragment.pdbname:
874 model = self._add_model(fragment)
878 models[-1].fragments.append(weakref.proxy(fragment))
882 pmi_offset = models[-1].asym_unit.entity.pmi_offset
883 sid_begin = min(fragment.start - pmi_offset,
884 models[-1].asym_unit.seq_id_range[0])
885 sid_end = max(fragment.end - pmi_offset,
886 models[-1].asym_unit.seq_id_range[1])
887 models[-1].asym_unit = fragment.asym_unit.asym(sid_begin, sid_end)
888 fragment.starting_model = models[-1]
890 def _add_model(self, f):
891 parser = ihm.metadata.PDBParser()
892 r = parser.parse_file(f.pdbname)
894 self.simo._add_dataset(r[
'dataset'])
896 templates = r[
'templates'].get(f.chain, [])
899 self.simo.system.locations.append(t.alignment_file)
901 self.simo._add_dataset(t.dataset)
902 source = r[
'entity_source'].get(f.chain)
904 f.asym_unit.entity.source = source
905 pmi_offset = f.asym_unit.entity.pmi_offset
907 asym_unit=f.asym_unit.asym.pmi_range(f.start, f.end),
908 dataset=r[
'dataset'], asym_id=f.chain,
909 templates=templates, offset=f.offset + pmi_offset,
910 metadata=r[
'metadata'],
911 software=r[
'software'][0]
if r[
'software']
else None,
912 script_file=r[
'script'])
913 m.fragments = [weakref.proxy(f)]
917 class _StartingModel(ihm.startmodel.StartingModel):
918 def get_seq_dif(self):
922 pmi_offset = self.asym_unit.entity.pmi_offset
923 mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
925 for f
in self.fragments:
928 residue_indexes=list(range(f.start - f.offset,
929 f.end - f.offset + 1)))
930 for a
in mh.get_ihm_atoms(sel.get_selected_particles(),
931 f.offset - pmi_offset):
933 self._seq_dif = mh._seq_dif
936 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
937 """Post processing using AnalysisReplicaExchange0 macro"""
939 def __init__(self, rex, num_models_begin):
942 for fname
in self.get_all_stat_files():
943 with open(str(fname))
as fh:
944 num_models_end += len(fh.readlines())
945 super(_ReplicaExchangeAnalysisPostProcess, self).__init__(
946 feature=
'RMSD', num_models_begin=num_models_begin,
947 num_models_end=num_models_end)
949 def get_stat_file(self, cluster_num):
950 return self.rex._outputdir / (
"cluster.%d" % cluster_num) /
'stat.out'
952 def get_all_stat_files(self):
953 for i
in range(self.rex._number_of_clusters):
954 yield self.get_stat_file(i)
957 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
958 """Ensemble generated using AnalysisReplicaExchange0 macro"""
960 num_models_deposited =
None
962 def __init__(self, pp, cluster_num, model_group, num_deposit):
963 with open(str(pp.get_stat_file(cluster_num)))
as fh:
964 num_models = len(fh.readlines())
965 super(_ReplicaExchangeAnalysisEnsemble, self).__init__(
966 num_models=num_models,
967 model_group=model_group, post_process=pp,
968 clustering_feature=pp.feature,
969 name=model_group.name)
970 self.cluster_num = cluster_num
971 self.num_models_deposited = num_deposit
973 def get_rmsf_file(self, component):
974 return (self.post_process.rex._outputdir
975 / (
'cluster.%d' % self.cluster_num)
976 / (
'rmsf.%s.dat' % component))
978 def load_rmsf(self, model, component):
979 fname = self.get_rmsf_file(component)
981 model.parse_rmsf_file(fname, component)
983 def get_localization_density_file(self, fname):
984 return (self.post_process.rex._outputdir
985 / (
'cluster.%d' % self.cluster_num)
986 / (
'%s.mrc' % fname))
988 def load_localization_density(self, state, fname, select_tuple,
990 fullpath = self.get_localization_density_file(fname)
991 if fullpath.exists():
992 details =
"Localization density for %s %s" \
993 % (fname, self.model_group.name)
994 local_file = ihm.location.OutputFileLocation(str(fullpath),
996 for s
in select_tuple:
997 if isinstance(s, tuple)
and len(s) == 3:
998 asym = asym_units[s[2]].pmi_range(s[0], s[1])
1000 asym = asym_units[s]
1001 den = ihm.model.LocalizationDensity(file=local_file,
1003 self.densities.append(den)
1005 def load_all_models(self, simo, state):
1006 stat_fname = self.post_process.get_stat_file(self.cluster_num)
1008 with open(str(stat_fname))
as fh:
1009 stats = ast.literal_eval(fh.readline())
1011 rmf_file = stat_fname.parent / (
"%d.rmf3" % model_num)
1013 if rmf_file.exists():
1014 rh = RMF.open_rmf_file_read_only(str(rmf_file))
1015 system = state._pmi_object.system
1021 if model_num >= self.num_models_deposited:
1025 def _get_precision(self):
1026 precfile = (self.post_process.rex._outputdir /
1027 (
"precision.%d.%d.out" % (self.cluster_num,
1029 if not precfile.exists():
1033 r'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
1035 with open(str(precfile))
as fh:
1039 return float(m.group(1))
1041 precision = property(
lambda self: self._get_precision(),
1042 lambda self, val:
None)
1045 class _SimpleEnsemble(ihm.model.Ensemble):
1046 """Simple manually-created ensemble"""
1048 num_models_deposited =
None
1050 def __init__(self, pp, model_group, num_models, drmsd,
1051 num_models_deposited, ensemble_file):
1052 super(_SimpleEnsemble, self).__init__(
1053 model_group=model_group, post_process=pp, num_models=num_models,
1054 file=ensemble_file, precision=drmsd, name=model_group.name,
1055 clustering_feature=
'dRMSD')
1056 self.num_models_deposited = num_models_deposited
1058 def load_localization_density(self, state, component, asym, local_file):
1059 den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
1060 self.densities.append(den)
1063 class _CustomDNAAlphabet(object):
1064 """Custom DNA alphabet that maps A,C,G,T (rather than DA,DC,DG,DT
1065 as in python-ihm)"""
1066 _comps = dict([cc.code_canonical, cc]
1067 for cc
in ihm.DNAAlphabet._comps.values())
1070 class _EntityMapper(dict):
1071 """Handle mapping from IMP components (without copy number) to CIF
1072 entities. Multiple components may map to the same entity if they
1074 def __init__(self, system):
1075 super(_EntityMapper, self).__init__()
1076 self._sequence_dict = {}
1078 self.system = system
1080 def _get_alphabet(self, alphabet):
1081 """Map a PMI alphabet to an IHM alphabet"""
1085 alphabet_map = {
None: ihm.LPeptideAlphabet,
1086 IMP.pmi.alphabets.amino_acid: ihm.LPeptideAlphabet,
1087 IMP.pmi.alphabets.rna: ihm.RNAAlphabet,
1088 IMP.pmi.alphabets.dna: _CustomDNAAlphabet}
1089 if alphabet
in alphabet_map:
1090 return alphabet_map[alphabet]
1092 raise TypeError(
"Don't know how to handle %s" % alphabet)
1094 def add(self, component_name, sequence, offset, alphabet):
1095 def entity_seq(sequence):
1098 return [
'UNK' if s ==
'X' else s
for s
in sequence]
1101 if sequence
not in self._sequence_dict:
1104 d = component_name.split(
"@")[0].split(
".")[0]
1105 entity = Entity(entity_seq(sequence), description=d,
1107 alphabet=self._get_alphabet(alphabet))
1108 self.system.entities.append(entity)
1109 self._sequence_dict[sequence] = entity
1110 self[component_name] = self._sequence_dict[sequence]
1113 class _TransformedComponent(object):
1114 def __init__(self, name, original, transform):
1115 self.name, self.original, self.transform = name, original, transform
1118 class _SimpleRef(object):
1119 """Class with similar interface to weakref.ref, but keeps a strong ref"""
1120 def __init__(self, ref):
1127 class _State(ihm.model.State):
1128 """Representation of a single state in the system."""
1130 def __init__(self, pmi_object, po):
1135 self._pmi_object = weakref.proxy(pmi_object)
1136 if hasattr(pmi_object,
'state'):
1139 self._pmi_state = _SimpleRef(pmi_object.state)
1141 self._pmi_state = weakref.ref(pmi_object)
1143 old_name = self.name
1144 super(_State, self).__init__(experiment_type=
'Fraction of bulk')
1145 self.name = old_name
1149 self.modeled_assembly = ihm.Assembly(
1150 name=
"Modeled assembly",
1151 description=self.get_postfixed_name(
1152 "All components modeled by IMP"))
1153 po.system.orphan_assemblies.append(self.modeled_assembly)
1155 self.all_modeled_components = []
1158 return hash(self._pmi_state())
1160 def __eq__(self, other):
1161 return self._pmi_state() == other._pmi_state()
1163 def add_model_group(self, group):
1167 def get_prefixed_name(self, name):
1168 """Prefix the given name with the state name, if available."""
1170 return self.short_name +
' ' + name
1174 return name[0].upper() + name[1:]
if name
else ''
1176 def get_postfixed_name(self, name):
1177 """Postfix the given name with the state name, if available."""
1179 return "%s in state %s" % (name, self.short_name)
1183 short_name = property(
lambda self: self._pmi_state().short_name)
1184 long_name = property(
lambda self: self._pmi_state().long_name)
1186 def __get_name(self):
1187 return self._pmi_state().long_name
1189 def __set_name(self, val):
1190 self._pmi_state().long_name = val
1192 name = property(__get_name, __set_name)
1196 """A single entity in the system.
1197 This functions identically to the base ihm.Entity class, but it
1198 allows identifying residues by either the IHM numbering scheme
1199 (seq_id, which is always contiguous starting at 1) or by the PMI
1200 scheme (which is similar, but may not start at 1 if the sequence in
1201 the FASTA file has one or more N-terminal gaps"""
1202 def __init__(self, sequence, pmi_offset, *args, **keys):
1205 self.pmi_offset = pmi_offset
1206 super(Entity, self).__init__(sequence, *args, **keys)
1209 """Return a single IHM residue indexed using PMI numbering"""
1210 return self.residue(res_id - self.pmi_offset)
1213 """Return a range of IHM residues indexed using PMI numbering"""
1214 off = self.pmi_offset
1215 return self(res_id_begin - off, res_id_end - off)
1219 """A single asymmetric unit in the system.
1220 This functions identically to the base ihm.AsymUnit class, but it
1221 allows identifying residues by either the IHM numbering scheme
1222 (seq_id, which is always contiguous starting at 1) or by the PMI
1223 scheme (which is similar, but may not start at 1 if the sequence in
1224 the FASTA file has one or more N-terminal gaps"""
1226 def __init__(self, entity, *args, **keys):
1227 super(AsymUnit, self).__init__(
1228 entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1231 """Return a single IHM residue indexed using PMI numbering"""
1232 return self.residue(res_id - self.entity.pmi_offset)
1235 """Return a range of IHM residues indexed using PMI numbering"""
1236 off = self.entity.pmi_offset
1237 return self(res_id_begin - off, res_id_end - off)
1241 """Class to encode a modeling protocol as mmCIF.
1243 IMP has basic support for writing out files in mmCIF format, for
1244 deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1245 After creating an instance of this class, attach it to an
1246 IMP.pmi.topology.System object. After this, any
1247 generated models and metadata are automatically collected in the
1248 `system` attribute, which is an
1249 [ihm.System](https://python-ihm.readthedocs.io/en/latest/main.html#ihm.System) object.
1250 Once the protocol is complete, call finalize() to make sure `system`
1251 contains everything, then use the
1252 [python-ihm API](https://python-ihm.readthedocs.io/en/latest/dumper.html#ihm.dumper.write)
1253 to write out files in mmCIF or BinaryCIF format.
1255 See also get_handlers(), get_dumpers().
1259 self.system = ihm.System()
1260 self._state_group = ihm.model.StateGroup()
1261 self.system.state_groups.append(self._state_group)
1263 self._state_ensemble_offset = 0
1264 self._main_script = os.path.abspath(sys.argv[0])
1267 loc = ihm.location.WorkflowFileLocation(
1268 path=self._main_script,
1269 details=
"The main integrative modeling script")
1270 self.system.locations.append(loc)
1273 self.__asym_states = {}
1274 self._working_directory = os.getcwd()
1276 "Default representation")
1277 self.entities = _EntityMapper(self.system)
1279 self.asym_units = {}
1280 self._all_components = {}
1281 self.all_modeled_components = []
1282 self._transformed_components = []
1283 self.sequence_dict = {}
1286 self._xy_plane = ihm.geometry.XYPlane()
1287 self._xz_plane = ihm.geometry.XZPlane()
1288 self._z_axis = ihm.geometry.ZAxis()
1289 self._center_origin = ihm.geometry.Center(0, 0, 0)
1290 self._identity_transform = ihm.geometry.Transformation.identity()
1293 self._exclude_coords = {}
1295 self.all_representations = _AllModelRepresentations(self)
1296 self.all_protocols = _AllProtocols(self)
1297 self.all_datasets = _AllDatasets(self.system)
1298 self.all_starting_models = _AllStartingModels(self)
1300 self.all_software = _AllSoftware(self.system)
1303 """Create a new Representation and return it. This can be
1304 passed to add_model(), add_bead_element() or add_pdb_element()."""
1305 r = ihm.representation.Representation(name=name)
1306 self.system.orphan_representations.append(r)
1310 """Don't record coordinates for the given domain.
1311 Coordinates for the given domain (specified by a component name
1312 and a 2-element tuple giving the start and end residue numbers)
1313 will be excluded from the mmCIF file. This can be used to exclude
1314 parts of the structure that weren't well resolved in modeling.
1315 Any bead or residue that lies wholly within this range will be
1316 excluded. Multiple ranges for a given component can be excluded
1317 by calling this method multiple times."""
1318 if component
not in self._exclude_coords:
1319 self._exclude_coords[component] = []
1320 self._exclude_coords[component].append(seqrange)
1322 def _is_excluded(self, component, start, end):
1323 """Return True iff this chunk of sequence should be excluded"""
1324 for seqrange
in self._exclude_coords.get(component, ()):
1325 if start >= seqrange[0]
and end <= seqrange[1]:
1328 def _add_state(self, state):
1329 """Create a new state and return a pointer to it."""
1330 self._state_ensemble_offset = len(self.system.ensembles)
1331 s = _State(state, self)
1332 self._state_group.append(s)
1333 self._last_state = s
1336 def _get_chain_for_component(self, name, output):
1337 """Get the chain ID for a component, if any."""
1339 if name
in self.asym_units:
1340 return self.asym_units[name]._id
1345 def _get_assembly_comps(self, assembly):
1346 """Get the names of the components in the given assembly"""
1350 comps[ca.details] =
None
1354 """Make a new component that's a transformed copy of another.
1355 All representation for the existing component is copied to the
1357 assembly_comps = self._get_assembly_comps(state.modeled_assembly)
1358 if name
in assembly_comps:
1359 raise ValueError(
"Component %s already exists" % name)
1360 elif original
not in assembly_comps:
1361 raise ValueError(
"Original component %s does not exist" % original)
1362 self.create_component(state, name,
True)
1363 self.add_component_sequence(state, name, self.sequence_dict[original])
1364 self._transformed_components.append(_TransformedComponent(
1365 name, original, transform))
1366 self.all_representations.copy_component(state, name, original,
1367 self.asym_units[name])
1369 def create_component(self, state, name, modeled, asym_name=None):
1370 if asym_name
is None:
1372 new_comp = name
not in self._all_components
1373 self._all_components[name] =
None
1375 state.all_modeled_components.append(name)
1376 if asym_name
not in self.asym_units:
1378 self.asym_units[asym_name] =
None
1380 self.all_modeled_components.append(name)
1382 def add_component_sequence(self, state, name, seq, asym_name=None,
1384 if asym_name
is None:
1387 def get_offset(seq):
1389 for i
in range(len(seq)):
1392 seq, offset = get_offset(seq)
1393 if name
in self.sequence_dict:
1394 if self.sequence_dict[name] != seq:
1395 raise ValueError(
"Sequence mismatch for component %s" % name)
1397 self.sequence_dict[name] = seq
1398 self.entities.add(name, seq, offset, alphabet)
1399 if asym_name
in self.asym_units:
1400 if self.asym_units[asym_name]
is None:
1402 entity = self.entities[name]
1403 asym =
AsymUnit(entity, details=asym_name)
1404 self.system.asym_units.append(asym)
1405 self.asym_units[asym_name] = asym
1406 state.modeled_assembly.append(self.asym_units[asym_name])
1408 def _add_restraint_model_fits(self):
1409 """Add fits to restraints for all known models"""
1410 for group, m
in self.system._all_models():
1411 if m._is_restrained:
1412 for r
in self.system.restraints:
1413 if hasattr(r,
'add_fits_from_model_statfile'):
1414 r.add_fits_from_model_statfile(m)
1417 """Do any final processing on the class hierarchy.
1418 After calling this method, the `system` member (an instance
1419 of `ihm.System`) completely reproduces the PMI modeling, and
1420 can be written out to an mmCIF file with `ihm.dumper.write`,
1421 and/or modified using the ihm API."""
1422 self._add_restraint_model_fits()
1424 def add_pdb_element(self, state, name, start, end, offset, pdbname,
1425 chain, hier, representation=
None):
1426 if self._is_excluded(name, start, end):
1428 if representation
is None:
1429 representation = self.default_representation
1430 asym = self.asym_units[name]
1431 p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1433 self.all_representations.add_fragment(state, representation, p)
1434 self.all_starting_models.add_pdb_fragment(p)
1436 def add_bead_element(self, state, name, start, end, num, hier,
1437 representation=
None):
1438 if self._is_excluded(name, start, end):
1440 if representation
is None:
1441 representation = self.default_representation
1442 asym = self.asym_units[name]
1443 pmi_offset = asym.entity.pmi_offset
1444 b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1446 self.all_representations.add_fragment(state, representation, b)
1448 def get_cross_link_group(self, pmi_restraint):
1449 r = _CrossLinkRestraint(pmi_restraint)
1450 self.system.restraints.append(r)
1451 self._add_restraint_dataset(r)
1454 def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1455 if c1
not in self._all_components
or c2
not in self._all_components:
1461 e1 = self.entities[c1]
1462 e2 = self.entities[c2]
1463 xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1464 residue2=e2.pmi_residue(r2))
1465 rsr.experimental_cross_links.append([xl])
1468 def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1471 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1472 d = ihm.restraint.UpperBoundDistanceRestraint(length)
1474 if _get_by_residue(p1)
and _get_by_residue(p2):
1475 cls = _ResidueCrossLink
1477 cls = _FeatureCrossLink
1478 xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1481 xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1482 rsr.cross_links.append(xl)
1484 def add_replica_exchange(self, state, rex):
1489 step = _ReplicaExchangeProtocolStep(state, rex)
1490 step.software = self.all_software.pmi
1491 self.all_protocols.add_step(step, state)
1493 def _add_simple_dynamics(self, num_models_end, method):
1495 state = self._last_state
1496 self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1499 def _add_protocol(self):
1501 state = self._last_state
1502 self.all_protocols.add_protocol(state)
1504 def _add_dataset(self, dataset):
1505 return self.all_datasets.add(self._last_state, dataset)
1507 def _add_restraint_dataset(self, restraint):
1508 return self.all_datasets.add_restraint(self._last_state, restraint)
1510 def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1512 state = self._last_state
1513 pp = ihm.analysis.ClusterStep(
'RMSD', num_models_begin, num_models_end)
1514 self.all_protocols.add_postproc(pp, state)
1517 def _add_no_postprocessing(self, num_models):
1519 state = self._last_state
1520 pp = ihm.analysis.EmptyStep()
1521 pp.num_models_begin = pp.num_models_end = num_models
1522 self.all_protocols.add_postproc(pp, state)
1525 def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1526 num_models_deposited, localization_densities,
1528 """Add an ensemble generated by ad hoc methods (not using PMI).
1529 This is currently only used by the Nup84 system."""
1531 state = self._last_state
1532 group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1533 state.add_model_group(group)
1535 self.system.locations.append(ensemble_file)
1536 e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1538 self.system.ensembles.append(e)
1539 for c
in state.all_modeled_components:
1540 den = localization_densities.get(c,
None)
1542 e.load_localization_density(state, c, self.asym_units[c], den)
1546 """Point a previously-created ensemble to an 'all-models' file.
1547 This could be a trajectory such as DCD, an RMF, or a multimodel
1549 self.system.locations.append(location)
1551 ind = i + self._state_ensemble_offset
1552 self.system.ensembles[ind].file = location
1554 def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1560 protocol = self.all_protocols.get_last_protocol(state)
1561 num_models = protocol.steps[-1].num_models_end
1562 pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1563 pp.software = self.all_software.pmi
1564 self.all_protocols.add_postproc(pp, state)
1565 for i
in range(rex._number_of_clusters):
1566 group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1567 'cluster %d' % (i + 1)))
1568 state.add_model_group(group)
1570 e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1571 self.system.ensembles.append(e)
1573 for fname, stuple
in sorted(density_custom_ranges.items()):
1574 e.load_localization_density(state, fname, stuple,
1576 for stats
in e.load_all_models(self, state):
1577 m = self.add_model(group)
1580 m.name =
'Best scoring model'
1583 for c
in state.all_modeled_components:
1586 def _get_subassembly(self, comps, name, description):
1587 """Get an Assembly consisting of the given components.
1588 `compdict` is a dictionary of the components to add, where keys
1589 are the component names and values are the sequence ranges (or
1590 None to use all residues in the component)."""
1592 for comp, seqrng
in comps.items():
1593 a = self.asym_units[comp]
1594 asyms.append(a
if seqrng
is None else a(*seqrng))
1596 a = ihm.Assembly(asyms, name=name, description=description)
1599 def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1601 """Add a basic FoXS fit. This is largely intended for use from the
1603 assembly = self._get_subassembly(
1605 name=
"SAXS subassembly",
1606 description=
"All components that fit SAXS data")
1607 r = ihm.restraint.SASRestraint(
1608 dataset, assembly, segment=
False,
1609 fitting_method=
'FoXS', fitting_atom_type=
'Heavy atoms',
1610 multi_state=
False, radius_of_gyration=rg, details=details)
1611 r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1612 self.system.restraints.append(r)
1613 self._add_restraint_dataset(r)
1615 def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1616 image_resolution, projection_number,
1617 micrographs_number):
1618 r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1619 image_resolution, projection_number,
1621 self.system.restraints.append(r)
1622 self._add_restraint_dataset(r)
1624 def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1626 r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1627 self.system.restraints.append(r)
1628 self._add_restraint_dataset(r)
1630 def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1631 sigma, pmi_restraint):
1632 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1633 sigma, pmi_restraint, self._xy_plane)
1635 def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1636 sigma, pmi_restraint):
1637 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1638 sigma, pmi_restraint, self._xz_plane)
1640 def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1641 sigma, pmi_restraint):
1642 self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1643 sigma, pmi_restraint, self._z_axis)
1645 def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1646 sigma, pmi_restraint, geom):
1647 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1648 r = _GeometricRestraint(
1649 self, state, pmi_restraint, geom, asym.get_feature(ps),
1650 ihm.restraint.LowerUpperBoundDistanceRestraint(lower_bound,
1653 self.system.restraints.append(r)
1654 self._add_restraint_dataset(r)
1656 def _get_membrane(self, tor_R, tor_r, tor_th):
1657 """Get an object representing a half-torus membrane"""
1658 if not hasattr(self,
'_seen_membranes'):
1659 self._seen_membranes = {}
1662 membrane_id = tuple(int(x * 100.)
for x
in (tor_R, tor_r, tor_th))
1663 if membrane_id
not in self._seen_membranes:
1664 m = ihm.geometry.HalfTorus(
1665 center=self._center_origin,
1666 transformation=self._identity_transform,
1667 major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1668 inner=
True, name=
'Membrane')
1669 self._seen_membranes[membrane_id] = m
1670 return self._seen_membranes[membrane_id]
1672 def add_membrane_surface_location_restraint(
1673 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1674 self._add_membrane_restraint(
1675 state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint,
1676 ihm.restraint.UpperBoundDistanceRestraint(0.))
1678 def add_membrane_exclusion_restraint(
1679 self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1680 self._add_membrane_restraint(
1681 state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint,
1682 ihm.restraint.LowerBoundDistanceRestraint(0.))
1684 def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1685 sigma, pmi_restraint, rsr):
1686 asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1687 r = _GeometricRestraint(
1688 self, state, pmi_restraint,
1689 self._get_membrane(tor_R, tor_r, tor_th), asym.get_feature(ps),
1691 self.system.restraints.append(r)
1692 self._add_restraint_dataset(r)
1694 def add_model(self, group, assembly=None, representation=None):
1695 state = self._last_state
1696 if representation
is None:
1697 representation = self.default_representation
1698 protocol = self.all_protocols.get_last_protocol(state)
1699 m = _Model(state.prot, self, protocol,
1700 assembly
if assembly
else state.modeled_assembly,
1707 """Get custom python-ihm dumpers for writing PMI to from mmCIF.
1708 This returns a list of custom dumpers that can be passed as all or
1709 part of the `dumpers` argument to ihm.dumper.write(). They add
1710 PMI-specific information to mmCIF or BinaryCIF files written out
1712 return [_ReplicaExchangeProtocolDumper]
1716 """Get custom python-ihm handlers for reading PMI data from mmCIF.
1717 This returns a list of custom handlers that can be passed as all or
1718 part of the `handlers` argument to ihm.reader.read(). They read
1719 PMI-specific information from mmCIF or BinaryCIF files read in
1721 return [_ReplicaExchangeProtocolHandler]
1725 """Extract metadata from an EM density GMM file."""
1728 """Extract metadata from `filename`.
1729 @return a dict with key `dataset` pointing to the GMM file and
1730 `number_of_gaussians` to the number of GMMs (or None)"""
1731 loc = ihm.location.InputFileLocation(
1733 details=
"Electron microscopy density map, "
1734 "represented as a Gaussian Mixture Model (GMM)")
1737 loc._allow_duplicates =
True
1738 d = ihm.dataset.EMDensityDataset(loc)
1739 ret = {
'dataset': d,
'number_of_gaussians':
None}
1741 with open(filename)
as fh:
1743 if line.startswith(
'# data_fn: '):
1744 p = ihm.metadata.MRCParser()
1745 fn = line[11:].rstrip(
'\r\n')
1746 dataset = p.parse_file(os.path.join(
1747 os.path.dirname(filename), fn))[
'dataset']
1748 ret[
'dataset'].parents.append(dataset)
1749 elif line.startswith(
'# ncenters: '):
1750 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 ...