IMP logo
IMP Reference Guide  2.20.2
The Integrative Modeling Platform
mmcif.py
1 """@namespace IMP.pmi.mmcif
2  @brief Support for the mmCIF file format.
3 
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.
10 """
11 
12 from __future__ import print_function
13 import copy
14 import RMF
15 import IMP.core
16 import IMP.algebra
17 import IMP.atom
18 import IMP.em
19 import IMP.isd
20 import IMP.pmi.tools
21 import IMP.pmi.alphabets
22 from IMP.pmi.tools import OrderedDict
23 import IMP.pmi.output
24 import IMP.mmcif.data
25 import re
26 import ast
27 import sys
28 import os
29 try:
30  from pathlib import Path
31 except ImportError: # Use bundled pathlib on Python 2 without pathlib
32  from IMP._compat_pathlib import Path # noqa: F401
33 import weakref
34 import ihm.location
35 import ihm.dataset
36 import ihm.dumper
37 import ihm.reader
38 import ihm.metadata
39 import ihm.startmodel
40 import ihm.model
41 import ihm.protocol
42 import ihm.analysis
43 import ihm.representation
44 import ihm.geometry
45 import ihm.cross_linkers
46 
47 
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'):
52  obj_by_id.append(obj)
53  obj.id = len(obj_by_id)
54  seen_objs[obj] = obj.id
55  else:
56  obj.id = seen_objs[obj]
57 
58 
59 def _get_by_residue(p):
60  """Determine whether the given particle represents a specific residue
61  or a more coarse-grained object."""
63 
64 
65 class _ComponentMapper(object):
66  """Map a Particle to a component name"""
67  def __init__(self, prot):
68  self.o = IMP.pmi.output.Output()
69  self.prot = 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)
74 
75  def __getitem__(self, p):
76  protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
77  return protname
78 
79 
80 class _AsymMapper(object):
81  """Map a Particle to an asym_unit"""
82  def __init__(self, simo, prot):
83  self.simo = simo
84  self._cm = _ComponentMapper(prot)
85  self._seen_ranges = {}
86 
87  def __getitem__(self, p):
88  protname = self._cm[p]
89  return self.simo.asym_units[protname]
90 
91  def get_feature(self, ps):
92  """Get an ihm.restraint.Feature that covers the given particles"""
93  # todo: handle things other than residues
94  rngs = []
95  for p in ps:
96  asym = self[p]
97  # todo: handle overlapping ranges
99  rind = IMP.atom.Residue(p).get_index()
100  rng = asym(rind, rind)
102  # PMI Fragments always contain contiguous residues
103  rinds = IMP.atom.Fragment(p).get_residue_indexes()
104  rng = asym(rinds[0], rinds[-1])
105  else:
106  raise ValueError("Unsupported particle type %s" % str(p))
107  # Join contiguous ranges
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],
111  rng.seq_id_range[1])
112  else:
113  rngs.append(rng)
114  # If an identical feature already exists, return that
115  # todo: python-ihm should handle this automatically for us
116  hrngs = tuple(rngs)
117  if hrngs in self._seen_ranges:
118  return self._seen_ranges[hrngs]
119  else:
120  feat = ihm.restraint.ResidueFeature(rngs)
121  self._seen_ranges[hrngs] = feat
122  return feat
123 
124 
125 class _AllSoftware(object):
126  def __init__(self, system):
127  self.system = 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')
141  # Only recent versions of python-ihm support adding citations for
142  # software
143  if hasattr(self.imp, 'citation'):
144  if sys.version_info[0] > 2:
145  # Don't include UTF8 characters in source; it confuses Python 2
146  javi = 'Vel\u00e1zquez-Muriel J'
147  else:
148  javi = 'Velazquez-Muriel J'
149  self.imp.citation = ihm.Citation(
150  pmid='22272186',
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',
155  year=2012,
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(
160  pmid='31396911',
161  title='Modeling Biological Complexes Using Integrative '
162  'Modeling Platform.',
163  journal='Methods Mol Biol', volume=2022, page_range=(353, 377),
164  year=2019,
165  authors=['Saltzberg D', 'Greenberg CH', 'Viswanath S',
166  'Chemmama I', 'Webb B', 'Pellarin R', 'Echeverria I',
167  'Sali A'],
168  doi='10.1007/978-1-4939-9608-7_15')
169  self.system.software.extend([self.pmi, self.imp])
170 
171  def set_modeller_used(self, version, date):
172  if self.modeller_used:
173  return
174  self.modeller_used = True
175  s = ihm.Software(
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(
183  pmid='8254673',
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')
189 
190  def set_phyre2_used(self):
191  if self.phyre2_used:
192  return
193  self.phyre2_used = True
194  s = ihm.Software(
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(
200  pmid='25950237',
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',
205  'Sternberg MJ'],
206  year=2015, doi='10.1038/nprot.2015.053')
207  self.system.software.append(s)
208 
209 
210 def _get_fragment_is_rigid(fragment):
211  """Determine whether a fragment is modeled rigidly"""
212  leaves = IMP.atom.get_leaves(fragment.hier)
213  # Assume all leaves are set up as rigid/flexible in the same way
214  return IMP.core.RigidMember.get_is_setup(leaves[0])
215 
216 
217 class _PDBFragment(ihm.representation.ResidueSegment):
218  """Record details about part of a PDB file used as input
219  for a component."""
220  def __init__(self, state, component, start, end, pdb_offset,
221  pdbname, chain, hier, asym_unit):
222  # start, end are PMI residue indexes (not IHM)
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
230  & IMP.atom.ChainPDBSelector([chain])
231  self.starting_hier = IMP.atom.read_pdb(pdbname, state.model, sel)
232 
233  rigid = property(lambda self: _get_fragment_is_rigid(self),
234  lambda self, val: None)
235 
236  def combine(self, other):
237  pass
238 
239 
240 class _BeadsFragment(ihm.representation.FeatureSegment):
241  """Record details about beads used to represent part of a component."""
242  chain = None
243 
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',
247  count=count)
248  self.state, self.component, self.hier = state, component, hier
249 
250  rigid = property(lambda self: _get_fragment_is_rigid(self),
251  lambda self, val: None)
252 
253  def combine(self, other):
254  # todo: don't combine if one fragment is rigid and the other flexible
255  if (type(other) == type(self) and # noqa: E721
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
261  return True
262 
263 
264 class _AllModelRepresentations(object):
265  def __init__(self, simo):
266  self.simo = simo
267  # dict of fragments, ordered by representation, then component name,
268  # then state
269  self.fragments = OrderedDict()
270  self._all_representations = {}
271 
272  def copy_component(self, state, name, original, asym_unit):
273  """Copy all representation for `original` in `state` to `name`"""
274  def copy_frag(f):
275  newf = copy.copy(f)
276  newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
277  return newf
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]]
284  # Assume representation for a component is the same in all
285  # states, so only write out the first one
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])
290 
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)
305  # Assume representation for a component is the same in all states,
306  # so only write out the first one
307  first_state = list(self.fragments[id_rep][comp].keys())[0]
308  if state is first_state:
309  representation.append(fragment)
310 
311 
312 class _AllDatasets(object):
313  """Track all datasets generated by PMI and add them to the ihm.System"""
314  def __init__(self, system):
315  self.system = system
316  self._datasets = []
317  self._datasets_by_state = {}
318  self._restraints_by_state = {}
319 
320  def get_all_group(self, state):
321  """Get a DatasetGroup encompassing all datasets so far in this state"""
322  # Note that if a restraint dataset is replaced after this is called,
323  # the group will still contain the old dataset - mark dataset as read
324  # only?
325  g = ihm.dataset.DatasetGroup(
326  self._datasets_by_state.get(state, [])
327  + [r.dataset for r in self._restraints_by_state.get(state, [])
328  if r.dataset])
329  return g
330 
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)
337  # Make sure that the dataset ends up in the mmCIF file even if nothing
338  # refers to it
339  self.system.orphan_datasets.append(dataset)
340  return dataset
341 
342  def add_restraint(self, state, restraint):
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)
347 
348 
349 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
350  """Restrain to a set of cross-links"""
351 
352  assembly = None
353  _label_map = {'wtDSS': 'DSS', 'scDSS': 'DSS', 'scEDC': 'EDC'}
354  _descriptor_map = {'DSS': ihm.cross_linkers.dss,
355  'EDC': ihm.cross_linkers.edc}
356 
357  def __init__(self, pmi_restraint):
358  self.pmi_restraint = pmi_restraint
359  # Use ChemDescriptor from PMI restraint if available, otherwise guess
360  # it using the label
361  linker = getattr(self.pmi_restraint, 'linker', None)
362  label = self.pmi_restraint.label
363  self.label = label
364  super(_CrossLinkRestraint, self).__init__(
365  dataset=self.pmi_restraint.dataset,
366  linker=linker or self._get_chem_descriptor(label))
367 
368  @classmethod
369  def _get_chem_descriptor(cls, label):
370  # Map commonly-used subtypes to more standard labels
371  label = cls._label_map.get(label, label)
372  if label not in cls._descriptor_map:
373  # If label is not a standard linker type, make a new chemical
374  # descriptor containing just the name. We don't know the chemistry
375  # so cannot specify a SMILES or INCHI string for it at this point
376  d = ihm.ChemDescriptor(label)
377  cls._descriptor_map[label] = d
378  return cls._descriptor_map[label]
379 
380  def _set_psi_sigma(self, model):
381  old_values = []
382  # Do nothing if given a different state than the restraint applies to
383  if model.m != self.pmi_restraint.model:
384  return
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()))
391  p.set_scale(sigma)
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()))
398  p.set_scale(psi)
399  # If the same particle is set multiple times, make sure we get the
400  # original value last
401  return list(reversed(old_values))
402 
403  def add_fits_from_model_statfile(self, model):
404  # Set psi/sigma for all particles
405  old_values = self._set_psi_sigma(model)
406  # If no stats were collected, we can't show the fit
407  if not old_values:
408  return
409  for xl in self.cross_links:
410  # Get current psi/sigma particle value for each XL
411  xl.fits[model] = ihm.restraint.CrossLinkFit(
412  psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
413  # Restore original psi/sigma
414  for p, scale in old_values:
415  p.set_scale(scale)
416 
417  # Have our dataset point to that in the original PMI restraint
418  def __set_dataset(self, val):
419  self.pmi_restraint.dataset = val
420  dataset = property(lambda self: self.pmi_restraint.dataset,
421  __set_dataset)
422 
423 
424 def get_asym_mapper_for_state(simo, state, asym_map):
425  asym = asym_map.get(state, None)
426  if asym is None:
427  asym = _AsymMapper(simo, state.prot)
428  asym_map[state] = asym
429  return asym
430 
431 
432 class _PMICrossLink(object):
433  # Query PMI particles to get psi and sigma values at time of use, not
434  # time of construction (since they may be set after creating the restraint)
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)
441 
442 
443 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
444  pass
445 
446 
447 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
448  pass
449 
450 
451 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
452  def __init__(self, state, pmi_restraint, image_number, resolution,
453  pixel_size, image_resolution, projection_number,
454  micrographs_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)
463 
464  # Have our dataset point to that in the original PMI restraint
465  def __get_dataset(self):
466  return self.pmi_restraint.datasets[self.image_number]
467 
468  def __set_dataset(self, val):
469  self.pmi_restraint.datasets[self.image_number] = val
470 
471  dataset = property(__get_dataset, __set_dataset)
472 
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,
480  rot_matrix=rm,
481  tr_vector=transform.get_translation())
482 
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])
490  for i in range(3)]
491  # If the model coordinates are transformed, need to back transform
492  # them first
493  inv = model.transform.get_inverse()
495  IMP.algebra.Vector3D(*t)) * inv
496 
497  def _get_cross_correlation(self, model):
498  """Get the cross correlation coefficient between the model projection
499  and the image"""
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)])
504 
505 
506 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
507 
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))
515 
516  # Have our dataset point to that in the original PMI restraint
517  def __set_dataset(self, val):
518  self.pmi_restraint.dataset = val
519  dataset = property(lambda self: self.pmi_restraint.dataset,
520  __set_dataset)
521 
522  def _get_assembly(self, densities, simo, state):
523  """Get the Assembly that this restraint acts on"""
524  cm = _ComponentMapper(state.prot)
525  components = {}
526  for d in densities:
527  components[cm[d]] = None # None == all residues in this component
528  a = simo._get_subassembly(
529  components, name="EM subassembly",
530  description="All components that fit the EM map")
531  return a
532 
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)
537 
538  def _get_cross_correlation(self, model):
539  """Get the cross correlation coefficient between the model
540  and the map"""
541  if model.stats is not None:
542  return float(model.stats['GaussianEMRestraint_%s_CCC'
543  % self.pmi_restraint.label])
544 
545 
546 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
547 
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,
555  restrain_all=True)
556 
557  # Have our dataset point to that in the original PMI restraint
558  def __set_dataset(self, val):
559  self.pmi_restraint.dataset = val
560  dataset = property(lambda self: self.pmi_restraint.dataset,
561  __set_dataset)
562 
563 
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'
568  else:
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,
577  dataset_group=None, # filled in by add_step()
578  method=method, name='Sampling',
579  num_models_begin=None, # filled in by add_step()
580  num_models_end=rex.vars["number_of_frames"],
581  multi_scale=True, multi_state=False, ordered=False, ensemble=True)
582 
583 
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():
594  for s in p.steps:
595  if isinstance(s, _ReplicaExchangeProtocolStep):
596  self._dump_step(p, s, lp)
597 
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)
605 
606 
607 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
608  category = '_imp_replica_exchange_protocol'
609 
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)
615  # todo: match step_id properly (don't assume contiguous)
616  s = p.steps[int(step_id)-1]
617  # Upgrade from plain ihm Step to IMP subclass
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)
625 
626 
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,
631  dataset_group=None, # filled in by add_step()
632  method=method, name='Sampling',
633  num_models_begin=None, # filled in by add_step()
634  num_models_end=num_models_end,
635  multi_scale=True, multi_state=False, ordered=False,
636  ensemble=True)
637 
638 
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
643  self.spheres = []
644  self.atoms = []
645 
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))
651  else:
652  self.atoms.append((xyz, atom_type, residue_type, residue_index,
653  all_indexes, radius))
654  orig_comp = property(lambda self: self.comp)
655 
656 
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
662 
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,
667  all_indexes, radius)
668  spheres = property(__get_spheres)
669 
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)
676 
677  entity = property(lambda self: self.orig_chain.entity)
678  orig_comp = property(lambda self: self.orig_chain.comp)
679 
680 
681 class _Excluder(object):
682  def __init__(self, component, simo):
683  self._seqranges = simo._exclude_coords.get(component, [])
684 
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]:
689  return True
690 
691 
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)
697  # Transformation from IMP coordinates into mmCIF coordinate system.
698  # Normally we pass through coordinates unchanged, but in some cases
699  # we may want to translate them (e.g. Nup84, where the deposited PDB
700  # files are all centered; we want the mmCIF files to match)
702  self.em2d_stats = None
703  self.stats = None
704  # True iff restraints act on this model
705  self._is_restrained = True
706  o = self.output = IMP.pmi.output.Output(atomistic=True)
707  name = 'cif-output'
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)
713  self.geometric_center = IMP.algebra.Vector3D(*self.geometric_center)
714  self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
715  self.rmsf = {}
716 
717  def all_chains(self, simo):
718  """Yield all chains, including transformed ones"""
719  chain_for_comp = {}
720  for c in self.chains:
721  yield c
722  chain_for_comp[c.comp] = c
723  for tc in simo._transformed_components:
724  orig_chain = chain_for_comp.get(tc.original, None)
725  if orig_chain:
726  asym = simo.asym_units[tc.name]
727  c = _TransformedChain(orig_chain, asym, tc.transform)
728  c.comp = tc.name
729  yield c
730 
731  def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
732  entity_for_chain = {}
733  comp_for_chain = {}
734  correct_asym = {}
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]
738  else:
739  # Remove copy number
740  pn = protname.split('.')[0]
741  entity_for_chain[chain_id] = simo.entities[pn]
742  comp_for_chain[chain_id] = protname
743  # When doing multi-state modeling, the chain ID returned here
744  # (assigned sequentially) might not be correct (states may have
745  # gaps in the chain IDs). Map it to the correct asym unit.
746  correct_asym[chain_id] = simo.asym_units[protname]
747 
748  # Gather by chain ID (should be sorted by chain ID already)
749  self.chains = []
750  chain = None
751  excluder = None
752 
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,
764  all_indexes, radius)
765 
766  def parse_rmsf_file(self, fname, component):
767  self.rmsf[component] = rmsf = {}
768  with open(str(fname)) as fh:
769  for line in fh:
770  resnum, blocknum, val = line.split()
771  rmsf[int(resnum)] = (int(blocknum), float(val))
772 
773  def get_rmsf(self, component, indexes):
774  """Get the RMSF value for the given residue indexes."""
775  if not self.rmsf:
776  return None
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"
781  % str(indexes))
782  return rmsf[indexes[0]][1]
783 
784  def get_atoms(self):
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(),
795  type_symbol=None, # todo: get element
796  x=pt[0], y=pt[1], z=pt[2])
797 
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))
813 
814 
815 class _AllProtocols(object):
816  def __init__(self, simo):
817  self.simo = weakref.proxy(simo)
818  # Lists of Protocols by state
819  self.protocols = OrderedDict()
820 
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)
828 
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
836  else:
837  step.num_models_begin = protocol.steps[-1].num_models_end
838  protocol.steps.append(step)
839  step.id = len(protocol.steps)
840  # Assume that protocol uses all currently-defined datasets
841  step.dataset_group = self.simo.all_datasets.get_all_group(state)
842 
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)
849 
850  def get_last_protocol(self, state):
851  """Return the most recently-added _Protocol"""
852  return self.protocols[state][-1]
853 
854 
855 class _AllStartingModels(object):
856  def __init__(self, simo):
857  self.simo = simo
858  # dict of starting models (entire PDB files), collected from fragments,
859  # ordered by component name and state
860  self.models = OrderedDict()
861  self.output = IMP.pmi.output.Output()
862 
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)
875  models.append(model)
876  else:
877  # Break circular ref between fragment and model
878  models[-1].fragments.append(weakref.proxy(fragment))
879  # Update residue range to cover all fragments
880  # We don't use fragment.offset here because fragment numbering
881  # is already in seq_id, not PDB numbering
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]
889 
890  def _add_model(self, f):
891  parser = ihm.metadata.PDBParser()
892  r = parser.parse_file(f.pdbname)
893 
894  self.simo._add_dataset(r['dataset'])
895  # We only want the templates that model the starting model chain
896  templates = r['templates'].get(f.chain, [])
897  for t in templates:
898  if t.alignment_file:
899  self.simo.system.locations.append(t.alignment_file)
900  if t.dataset:
901  self.simo._add_dataset(t.dataset)
902  source = r['entity_source'].get(f.chain)
903  if source:
904  f.asym_unit.entity.source = source
905  pmi_offset = f.asym_unit.entity.pmi_offset
906  m = _StartingModel(
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)]
914  return m
915 
916 
917 class _StartingModel(ihm.startmodel.StartingModel):
918  def get_seq_dif(self):
919  return self._seq_dif
920 
921  def get_atoms(self):
922  pmi_offset = self.asym_unit.entity.pmi_offset
923  mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
924  self.asym_unit)
925  for f in self.fragments:
926  sel = IMP.atom.Selection(
927  f.starting_hier,
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):
932  yield a
933  self._seq_dif = mh._seq_dif
934 
935 
936 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
937  """Post processing using AnalysisReplicaExchange0 macro"""
938 
939  def __init__(self, rex, num_models_begin):
940  self.rex = rex
941  num_models_end = 0
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)
948 
949  def get_stat_file(self, cluster_num):
950  return self.rex._outputdir / ("cluster.%d" % cluster_num) / 'stat.out'
951 
952  def get_all_stat_files(self):
953  for i in range(self.rex._number_of_clusters):
954  yield self.get_stat_file(i)
955 
956 
957 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
958  """Ensemble generated using AnalysisReplicaExchange0 macro"""
959 
960  num_models_deposited = None # Override base class property
961 
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
972 
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))
977 
978  def load_rmsf(self, model, component):
979  fname = self.get_rmsf_file(component)
980  if fname.exists():
981  model.parse_rmsf_file(fname, component)
982 
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))
987 
988  def load_localization_density(self, state, fname, select_tuple,
989  asym_units):
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),
995  details=details)
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])
999  else:
1000  asym = asym_units[s]
1001  den = ihm.model.LocalizationDensity(file=local_file,
1002  asym_unit=asym)
1003  self.densities.append(den)
1004 
1005  def load_all_models(self, simo, state):
1006  stat_fname = self.post_process.get_stat_file(self.cluster_num)
1007  model_num = 0
1008  with open(str(stat_fname)) as fh:
1009  stats = ast.literal_eval(fh.readline())
1010  # Correct path
1011  rmf_file = stat_fname.parent / ("%d.rmf3" % model_num)
1012  # todo: test with real PMI2 systems
1013  if rmf_file.exists():
1014  rh = RMF.open_rmf_file_read_only(str(rmf_file))
1015  system = state._pmi_object.system
1016  IMP.rmf.link_hierarchies(rh, [system.hier])
1017  IMP.rmf.load_frame(fh, RMF.FrameID(0))
1018  # todo: fill in other data from stat file, e.g. cross-link phi/psi
1019  yield stats
1020  model_num += 1
1021  if model_num >= self.num_models_deposited:
1022  return
1023 
1024  # todo: also support dRMS precision
1025  def _get_precision(self):
1026  precfile = (self.post_process.rex._outputdir /
1027  ("precision.%d.%d.out" % (self.cluster_num,
1028  self.cluster_num)))
1029  if not precfile.exists():
1030  return ihm.unknown
1031  # Fail if the precision.x.x.out file doesn't match the cluster
1032  r = re.compile(
1033  r'All .*/cluster.%d/ average centroid distance ([\d\.]+)'
1034  % self.cluster_num)
1035  with open(str(precfile)) as fh:
1036  for line in fh:
1037  m = r.match(line)
1038  if m:
1039  return float(m.group(1))
1040 
1041  precision = property(lambda self: self._get_precision(),
1042  lambda self, val: None)
1043 
1044 
1045 class _SimpleEnsemble(ihm.model.Ensemble):
1046  """Simple manually-created ensemble"""
1047 
1048  num_models_deposited = None # Override base class property
1049 
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
1057 
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)
1061 
1062 
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())
1068 
1069 
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
1073  share sequence."""
1074  def __init__(self, system):
1075  super(_EntityMapper, self).__init__()
1076  self._sequence_dict = {}
1077  self._entities = []
1078  self.system = system
1079 
1080  def _get_alphabet(self, alphabet):
1081  """Map a PMI alphabet to an IHM alphabet"""
1082  # todo: we should probably ignore alphabets entirely and use
1083  # residue types directly (e.g. we shouldn't compare one-letter
1084  # code sequence to determine if an Entity is unique)
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]
1091  else:
1092  raise TypeError("Don't know how to handle %s" % alphabet)
1093 
1094  def add(self, component_name, sequence, offset, alphabet):
1095  def entity_seq(sequence):
1096  # Map X to UNK
1097  if 'X' in sequence:
1098  return ['UNK' if s == 'X' else s for s in sequence]
1099  else:
1100  return sequence
1101  if sequence not in self._sequence_dict:
1102  # Use the name of the first component, stripped of any copy number,
1103  # as the description of the entity
1104  d = component_name.split("@")[0].split(".")[0]
1105  entity = Entity(entity_seq(sequence), description=d,
1106  pmi_offset=offset,
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]
1111 
1112 
1113 class _TransformedComponent(object):
1114  def __init__(self, name, original, transform):
1115  self.name, self.original, self.transform = name, original, transform
1116 
1117 
1118 class _SimpleRef(object):
1119  """Class with similar interface to weakref.ref, but keeps a strong ref"""
1120  def __init__(self, ref):
1121  self.ref = ref
1122 
1123  def __call__(self):
1124  return self.ref
1125 
1126 
1127 class _State(ihm.model.State):
1128  """Representation of a single state in the system."""
1129 
1130  def __init__(self, pmi_object, po):
1131  # Point to the PMI object for this state. Use a weak reference
1132  # since the state object typically points to us too, so we need
1133  # to break the reference cycle. In PMI1 this will be a
1134  # Representation object; in PMI2 it is the PMI2 State object itself.
1135  self._pmi_object = weakref.proxy(pmi_object)
1136  if hasattr(pmi_object, 'state'):
1137  # Need a strong ref to pmi_object.state to prevent it from being
1138  # cleaned up when doing PMI1 multistate modeling
1139  self._pmi_state = _SimpleRef(pmi_object.state)
1140  else:
1141  self._pmi_state = weakref.ref(pmi_object)
1142  # Preserve PMI state name
1143  old_name = self.name
1144  super(_State, self).__init__(experiment_type='Fraction of bulk')
1145  self.name = old_name
1146 
1147  # The assembly of all components modeled by IMP in this state.
1148  # This may be smaller than the complete assembly.
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)
1154 
1155  self.all_modeled_components = []
1156 
1157  def __hash__(self):
1158  return hash(self._pmi_state())
1159 
1160  def __eq__(self, other):
1161  return self._pmi_state() == other._pmi_state()
1162 
1163  def add_model_group(self, group):
1164  self.append(group)
1165  return group
1166 
1167  def get_prefixed_name(self, name):
1168  """Prefix the given name with the state name, if available."""
1169  if self.short_name:
1170  return self.short_name + ' ' + name
1171  else:
1172  # Can't use capitalize() since that un-capitalizes everything
1173  # but the first letter
1174  return name[0].upper() + name[1:] if name else ''
1175 
1176  def get_postfixed_name(self, name):
1177  """Postfix the given name with the state name, if available."""
1178  if self.short_name:
1179  return "%s in state %s" % (name, self.short_name)
1180  else:
1181  return name
1182 
1183  short_name = property(lambda self: self._pmi_state().short_name)
1184  long_name = property(lambda self: self._pmi_state().long_name)
1185 
1186  def __get_name(self):
1187  return self._pmi_state().long_name
1188 
1189  def __set_name(self, val):
1190  self._pmi_state().long_name = val
1191 
1192  name = property(__get_name, __set_name)
1193 
1194 
1195 class Entity(ihm.Entity):
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):
1203  # Offset between PMI numbering and IHM; <pmi_#> = <ihm_#> + pmi_offset
1204  # (pmi_offset is also the number of N-terminal gaps in the FASTA file)
1205  self.pmi_offset = pmi_offset
1206  super(Entity, self).__init__(sequence, *args, **keys)
1207 
1208  def pmi_residue(self, res_id):
1209  """Return a single IHM residue indexed using PMI numbering"""
1210  return self.residue(res_id - self.pmi_offset)
1211 
1212  def pmi_range(self, res_id_begin, res_id_end):
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)
1216 
1217 
1218 class AsymUnit(ihm.AsymUnit):
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"""
1225 
1226  def __init__(self, entity, *args, **keys):
1227  super(AsymUnit, self).__init__(
1228  entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1229 
1230  def pmi_residue(self, res_id):
1231  """Return a single IHM residue indexed using PMI numbering"""
1232  return self.residue(res_id - self.entity.pmi_offset)
1233 
1234  def pmi_range(self, res_id_begin, res_id_end):
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)
1238 
1239 
1241  """Class to encode a modeling protocol as mmCIF.
1242 
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.
1254 
1255  See also get_handlers(), get_dumpers().
1256  """ # noqa: E501
1257  def __init__(self):
1258  # Ultimately, collect data in an ihm.System object
1259  self.system = ihm.System()
1260  self._state_group = ihm.model.StateGroup()
1261  self.system.state_groups.append(self._state_group)
1262 
1263  self._state_ensemble_offset = 0
1264  self._main_script = os.path.abspath(sys.argv[0])
1265 
1266  # Point to the main modeling script
1267  loc = ihm.location.WorkflowFileLocation(
1268  path=self._main_script,
1269  details="The main integrative modeling script")
1270  self.system.locations.append(loc)
1271 
1272  self._states = {}
1273  self.__asym_states = {}
1274  self._working_directory = os.getcwd()
1275  self.default_representation = self.create_representation(
1276  "Default representation")
1277  self.entities = _EntityMapper(self.system)
1278  # Mapping from component names to ihm.AsymUnit
1279  self.asym_units = {}
1280  self._all_components = {}
1281  self.all_modeled_components = []
1282  self._transformed_components = []
1283  self.sequence_dict = {}
1284 
1285  # Common geometric objects used in PMI systems
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()
1291 
1292  # Coordinates to exclude
1293  self._exclude_coords = {}
1294 
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)
1299 
1300  self.all_software = _AllSoftware(self.system)
1301 
1302  def create_representation(self, name):
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)
1307  return r
1308 
1309  def exclude_coordinates(self, component, seqrange):
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)
1321 
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]:
1326  return True
1327 
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
1334  return s
1335 
1336  def _get_chain_for_component(self, name, output):
1337  """Get the chain ID for a component, if any."""
1338  # todo: handle multiple copies
1339  if name in self.asym_units:
1340  return self.asym_units[name]._id
1341  else:
1342  # A non-modeled component doesn't have a chain ID
1343  return None
1344 
1345  def _get_assembly_comps(self, assembly):
1346  """Get the names of the components in the given assembly"""
1347  # component name = asym_unit.details
1348  comps = {}
1349  for ca in assembly:
1350  comps[ca.details] = None
1351  return comps
1352 
1353  def create_transformed_component(self, state, name, original, transform):
1354  """Make a new component that's a transformed copy of another.
1355  All representation for the existing component is copied to the
1356  new one."""
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])
1368 
1369  def create_component(self, state, name, modeled, asym_name=None):
1370  if asym_name is None:
1371  asym_name = name
1372  new_comp = name not in self._all_components
1373  self._all_components[name] = None
1374  if modeled:
1375  state.all_modeled_components.append(name)
1376  if asym_name not in self.asym_units:
1377  # assign asym once we get sequence
1378  self.asym_units[asym_name] = None
1379  if new_comp:
1380  self.all_modeled_components.append(name)
1381 
1382  def add_component_sequence(self, state, name, seq, asym_name=None,
1383  alphabet=None):
1384  if asym_name is None:
1385  asym_name = name
1386 
1387  def get_offset(seq):
1388  # Count length of gaps at start of sequence
1389  for i in range(len(seq)):
1390  if seq[i] != '-':
1391  return seq[i:], i
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)
1396  else:
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:
1401  # Set up a new asymmetric unit for this component
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])
1407 
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)
1415 
1416  def finalize(self):
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()
1423 
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):
1427  return
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,
1432  hier, asym)
1433  self.all_representations.add_fragment(state, representation, p)
1434  self.all_starting_models.add_pdb_fragment(p)
1435 
1436  def add_bead_element(self, state, name, start, end, num, hier,
1437  representation=None):
1438  if self._is_excluded(name, start, end):
1439  return
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,
1445  num, hier, asym)
1446  self.all_representations.add_fragment(state, representation, b)
1447 
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) # so that all-dataset group works
1452  return r
1453 
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:
1456  # Crosslink refers to a component we didn't model
1457  # As a quick hack, just ignore it.
1458  # todo: need to add an entity for this anyway (so will need the
1459  # sequence, and add to struct_assembly)
1460  return None
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])
1466  return xl
1467 
1468  def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1469  psi_p, rsr):
1470  # Map p1/p2 to asym units
1471  asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1472  d = ihm.restraint.UpperBoundDistanceRestraint(length)
1473 
1474  if _get_by_residue(p1) and _get_by_residue(p2):
1475  cls = _ResidueCrossLink
1476  else:
1477  cls = _FeatureCrossLink
1478  xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1479  restrain_all=False)
1480  # Needed to look up fits later
1481  xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1482  rsr.cross_links.append(xl)
1483 
1484  def add_replica_exchange(self, state, rex):
1485  # todo: allow for metadata to say how many replicas were used in the
1486  # actual experiment, and how many independent runs were carried out
1487  # (use these as multipliers to get the correct total number of
1488  # output models)
1489  step = _ReplicaExchangeProtocolStep(state, rex)
1490  step.software = self.all_software.pmi
1491  self.all_protocols.add_step(step, state)
1492 
1493  def _add_simple_dynamics(self, num_models_end, method):
1494  # Always assumed that we're dealing with the last state
1495  state = self._last_state
1496  self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1497  method), state)
1498 
1499  def _add_protocol(self):
1500  # Always assumed that we're dealing with the last state
1501  state = self._last_state
1502  self.all_protocols.add_protocol(state)
1503 
1504  def _add_dataset(self, dataset):
1505  return self.all_datasets.add(self._last_state, dataset)
1506 
1507  def _add_restraint_dataset(self, restraint):
1508  return self.all_datasets.add_restraint(self._last_state, restraint)
1509 
1510  def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1511  # Always assumed that we're dealing with the last state
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)
1515  return pp
1516 
1517  def _add_no_postprocessing(self, num_models):
1518  # Always assumed that we're dealing with the last state
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)
1523  return pp
1524 
1525  def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1526  num_models_deposited, localization_densities,
1527  ensemble_file):
1528  """Add an ensemble generated by ad hoc methods (not using PMI).
1529  This is currently only used by the Nup84 system."""
1530  # Always assumed that we're dealing with the last state
1531  state = self._last_state
1532  group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1533  state.add_model_group(group)
1534  if ensemble_file:
1535  self.system.locations.append(ensemble_file)
1536  e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1537  ensemble_file)
1538  self.system.ensembles.append(e)
1539  for c in state.all_modeled_components:
1540  den = localization_densities.get(c, None)
1541  if den:
1542  e.load_localization_density(state, c, self.asym_units[c], den)
1543  return e
1544 
1545  def set_ensemble_file(self, i, location):
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
1548  PDB file."""
1549  self.system.locations.append(location)
1550  # Ensure that we point to an ensemble related to the current state
1551  ind = i + self._state_ensemble_offset
1552  self.system.ensembles[ind].file = location
1553 
1554  def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1555  # todo: add prefilter as an additional postprocess step (complication:
1556  # we don't know how many models it filtered)
1557  # todo: postpone rmsf/density localization extraction until after
1558  # relevant methods are called (currently it is done from the
1559  # constructor)
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)
1569  # todo: make # of models to deposit configurable somewhere
1570  e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1571  self.system.ensembles.append(e)
1572  # Add localization density info if available
1573  for fname, stuple in sorted(density_custom_ranges.items()):
1574  e.load_localization_density(state, fname, stuple,
1575  self.asym_units)
1576  for stats in e.load_all_models(self, state):
1577  m = self.add_model(group)
1578  # Since we currently only deposit 1 model, it is the
1579  # best scoring one
1580  m.name = 'Best scoring model'
1581  m.stats = stats
1582  # Add RMSF info if available
1583  for c in state.all_modeled_components:
1584  e.load_rmsf(m, c)
1585 
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)."""
1591  asyms = []
1592  for comp, seqrng in comps.items():
1593  a = self.asym_units[comp]
1594  asyms.append(a if seqrng is None else a(*seqrng))
1595 
1596  a = ihm.Assembly(asyms, name=name, description=description)
1597  return a
1598 
1599  def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1600  details):
1601  """Add a basic FoXS fit. This is largely intended for use from the
1602  NPC application."""
1603  assembly = self._get_subassembly(
1604  {comp: seqrange},
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) # so that all-dataset group works
1614 
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,
1620  micrographs_number)
1621  self.system.restraints.append(r)
1622  self._add_restraint_dataset(r) # so that all-dataset group works
1623 
1624  def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1625  # todo: need to set allow_duplicates on this dataset?
1626  r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1627  self.system.restraints.append(r)
1628  self._add_restraint_dataset(r) # so that all-dataset group works
1629 
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)
1634 
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)
1639 
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)
1644 
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,
1651  upper_bound),
1652  sigma)
1653  self.system.restraints.append(r)
1654  self._add_restraint_dataset(r) # so that all-dataset group works
1655 
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 = {}
1660  # If a membrane already exists with all dimensions within 0.01 of
1661  # this one, reuse it
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]
1671 
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.))
1677 
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.))
1683 
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),
1690  rsr, sigma)
1691  self.system.restraints.append(r)
1692  self._add_restraint_dataset(r) # so that all-dataset group works
1693 
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,
1701  representation)
1702  group.append(m)
1703  return m
1704 
1705 
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
1711  by python-ihm."""
1712  return [_ReplicaExchangeProtocolDumper]
1713 
1714 
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
1720  by python-ihm."""
1721  return [_ReplicaExchangeProtocolHandler]
1722 
1723 
1724 class GMMParser(ihm.metadata.Parser):
1725  """Extract metadata from an EM density GMM file."""
1726 
1727  def parse_file(self, filename):
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(
1732  filename,
1733  details="Electron microscopy density map, "
1734  "represented as a Gaussian Mixture Model (GMM)")
1735  # A 3DEM restraint's dataset ID uniquely defines the mmCIF
1736  # restraint, so we need to allow duplicates
1737  loc._allow_duplicates = True
1738  d = ihm.dataset.EMDensityDataset(loc)
1739  ret = {'dataset': d, 'number_of_gaussians': None}
1740 
1741  with open(filename) as fh:
1742  for line in 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:])
1751  return ret
Select non water and non hydrogen atoms.
Definition: pdb.h:331
def get_handlers
Get custom python-ihm handlers for reading PMI data from mmCIF.
Definition: mmcif.py:1715
def create_representation
Create a new Representation and return it.
Definition: mmcif.py:1302
def create_transformed_component
Make a new component that's a transformed copy of another.
Definition: mmcif.py:1353
Simple 3D transformation class.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:158
def exclude_coordinates
Don't record coordinates for the given domain.
Definition: mmcif.py:1312
A decorator to associate a particle with a part of a protein/DNA/RNA.
Definition: Fragment.h:20
Class to encode a modeling protocol as mmCIF.
Definition: mmcif.py:1253
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: atom/Atom.h:245
def parse_file
Extract metadata from filename.
Definition: mmcif.py:1728
def pmi_range
Return a range of IHM residues indexed using PMI numbering.
Definition: mmcif.py:1212
Miscellaneous utilities.
Definition: tools.py:1
Extract metadata from an EM density GMM file.
Definition: mmcif.py:1724
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:749
def pmi_residue
Return a single IHM residue indexed using PMI numbering.
Definition: mmcif.py:1230
void read_pdb(TextInput input, int model, Hierarchy h)
A single asymmetric unit in the system.
Definition: mmcif.py:1222
A single entity in the system.
Definition: mmcif.py:1195
def set_ensemble_file
Point a previously-created ensemble to an 'all-models' file.
Definition: mmcif.py:1545
Classes to represent data structures used in mmCIF.
Definition: data.py:1
def pmi_range
Return a range of IHM residues indexed using PMI numbering.
Definition: mmcif.py:1234
void add_restraint(RMF::FileHandle fh, Restraint *hs)
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def finalize
Do any final processing on the class hierarchy.
Definition: mmcif.py:1416
Base class for capturing a modeling protocol.
Definition: output.py:45
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.
Definition: mmcif.py:1706
Class for easy writing of PDBs, RMFs, and stat files.
Definition: output.py:203
3D rotation class.
Definition: Rotation3D.h:52
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Classes for writing output files and processing them.
Definition: output.py:1
A decorator for a residue.
Definition: Residue.h:137
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.
Definition: mmcif.py:1208
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
VectorD< 3 > Vector3D
Definition: VectorD.h:425
Mapping between FASTA one-letter codes and residue types.
Definition: alphabets.py:1
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.
Definition: Selection.h:70
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:256
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...