IMP logo
IMP Reference Guide  2.13.0
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 import textwrap
30 import weakref
31 import operator
32 import itertools
33 import ihm.location
34 import ihm.dataset
35 import ihm.dumper
36 import ihm.reader
37 import ihm.metadata
38 import ihm.startmodel
39 import ihm.model
40 import ihm.protocol
41 import ihm.analysis
42 import ihm.representation
43 import ihm.geometry
44 import ihm.cross_linkers
45 
46 def _assign_id(obj, seen_objs, obj_by_id):
47  """Assign a unique ID to obj, and track all ids in obj_by_id."""
48  if obj not in seen_objs:
49  if not hasattr(obj, 'id'):
50  obj_by_id.append(obj)
51  obj.id = len(obj_by_id)
52  seen_objs[obj] = obj.id
53  else:
54  obj.id = seen_objs[obj]
55 
56 
57 def _get_by_residue(p):
58  """Determine whether the given particle represents a specific residue
59  or a more coarse-grained object."""
61 
62 
63 class _ComponentMapper(object):
64  """Map a Particle to a component name"""
65  def __init__(self, prot):
66  self.o = IMP.pmi.output.Output()
67  self.prot = prot
68  self.name = 'cif-output'
69  self.o.dictionary_pdbs[self.name] = self.prot
70  self.o._init_dictchain(self.name, self.prot, multichar_chain=True)
71 
72  def __getitem__(self, p):
73  protname, is_a_bead = self.o.get_prot_name_from_particle(self.name, p)
74  return protname
75 
76 
77 class _AsymMapper(object):
78  """Map a Particle to an asym_unit"""
79  def __init__(self, simo, prot):
80  self.simo = simo
81  self._cm = _ComponentMapper(prot)
82  self._seen_ranges = {}
83 
84  def __getitem__(self, p):
85  protname = self._cm[p]
86  return self.simo.asym_units[protname]
87 
88  def get_feature(self, ps):
89  """Get an ihm.restraint.Feature that covers the given particles"""
90  # todo: handle things other than residues
91  rngs = []
92  for p in ps:
93  asym = self[p]
94  # todo: handle overlapping ranges
96  rind = IMP.atom.Residue(p).get_index()
97  rng = asym(rind, rind)
99  # PMI Fragments always contain contiguous residues
100  rinds = IMP.atom.Fragment(p).get_residue_indexes()
101  rng = asym(rinds[0], rinds[-1])
102  else:
103  raise ValueError("Unsupported particle type %s" % str(p))
104  # Join contiguous ranges
105  if len(rngs) > 0 and rngs[-1].asym == asym \
106  and rngs[-1].seq_id_range[1] == rng.seq_id_range[0] - 1:
107  rngs[-1].seq_id_range = (rngs[-1].seq_id_range[0],
108  rng.seq_id_range[1])
109  else:
110  rngs.append(rng)
111  # If an identical feature already exists, return that
112  # todo: python-ihm should handle this automatically for us
113  hrngs = tuple(rngs)
114  if hrngs in self._seen_ranges:
115  return self._seen_ranges[hrngs]
116  else:
117  feat = ihm.restraint.ResidueFeature(rngs)
118  self._seen_ranges[hrngs] = feat
119  return feat
120 
121 
122 class _AllSoftware(object):
123  def __init__(self, system):
124  self.system = system
125  self.modeller_used = self.phyre2_used = False
126  self.pmi = ihm.Software(
127  name="IMP PMI module",
128  version=IMP.pmi.__version__,
129  classification="integrative model building",
130  description="integrative model building",
131  location='https://integrativemodeling.org')
132  self.imp = ihm.Software(
133  name="Integrative Modeling Platform (IMP)",
134  version=IMP.__version__,
135  classification="integrative model building",
136  description="integrative model building",
137  location='https://integrativemodeling.org')
138  self.system.software.extend([self.pmi, self.imp])
139 
140  def set_modeller_used(self, version, date):
141  if self.modeller_used:
142  return
143  self.modeller_used = True
144  self.system.software.append(ihm.Software(
145  name='MODELLER', classification='comparative modeling',
146  description='Comparative modeling by satisfaction '
147  'of spatial restraints, build ' + date,
148  location='https://salilab.org/modeller/',
149  version=version))
150 
151  def set_phyre2_used(self):
152  if self.phyre2_used:
153  return
154  self.phyre2_used = True
155  self.system.software.append(ihm.Software(
156  name='Phyre2', classification='protein homology modeling',
157  description='Protein Homology/analogY Recognition '
158  'Engine V 2.0',
159  version='2.0',
160  location='http://www.sbg.bio.ic.ac.uk/~phyre2/'))
161 
162 def _get_fragment_is_rigid(fragment):
163  """Determine whether a fragment is modeled rigidly"""
164  leaves = IMP.atom.get_leaves(fragment.hier)
165  # Assume all leaves are set up as rigid/flexible in the same way
166  return IMP.core.RigidMember.get_is_setup(leaves[0])
167 
168 class _PDBFragment(ihm.representation.ResidueSegment):
169  """Record details about part of a PDB file used as input
170  for a component."""
171  def __init__(self, state, component, start, end, pdb_offset,
172  pdbname, chain, hier, asym_unit):
173  # start, end are PMI residue indexes (not IHM)
174  super(_PDBFragment, self).__init__(
175  asym_unit=asym_unit.pmi_range(start, end),
176  rigid=None, primitive='sphere')
177  self.component, self.start, self.end, self.offset, self.pdbname \
178  = component, start, end, pdb_offset, pdbname
179  self.state, self.chain, self.hier = state, chain, hier
182  self.starting_hier = IMP.atom.read_pdb(pdbname, state.model, sel)
183 
184  rigid = property(lambda self: _get_fragment_is_rigid(self),
185  lambda self, val: None)
186 
187  def combine(self, other):
188  pass
189 
190 
191 class _BeadsFragment(ihm.representation.FeatureSegment):
192  """Record details about beads used to represent part of a component."""
193  chain = None
194  def __init__(self, state, component, start, end, count, hier, asym_unit):
195  super(_BeadsFragment, self).__init__(asym_unit=asym_unit(start, end),
196  rigid=None, primitive='sphere', count=count)
197  self.state, self.component, self.hier \
198  = state, component, hier
199 
200  rigid = property(lambda self: _get_fragment_is_rigid(self),
201  lambda self, val: None)
202 
203  def combine(self, other):
204  # todo: don't combine if one fragment is rigid and the other flexible
205  if type(other) == type(self) and \
206  other.asym_unit.seq_id_range[0] == self.asym_unit.seq_id_range[1] + 1:
207  self.asym_unit.seq_id_range = (self.asym_unit.seq_id_range[0],
208  other.asym_unit.seq_id_range[1])
209  self.count += other.count
210  return True
211 
212 
213 class _AllModelRepresentations(object):
214  def __init__(self, simo):
215  self.simo = simo
216  # dict of fragments, ordered by representation, then component name,
217  # then state
218  self.fragments = OrderedDict()
219  self._all_representations = {}
220 
221  def copy_component(self, state, name, original, asym_unit):
222  """Copy all representation for `original` in `state` to `name`"""
223  def copy_frag(f):
224  newf = copy.copy(f)
225  newf.asym_unit = asym_unit(*f.asym_unit.seq_id_range)
226  return newf
227  for rep in self.fragments:
228  if original in self.fragments[rep]:
229  if name not in self.fragments[rep]:
230  self.fragments[rep][name] = OrderedDict()
231  self.fragments[rep][name][state] = [copy_frag(f)
232  for f in self.fragments[rep][original][state]]
233  # Assume representation for a component is the same in all
234  # states, so only write out the first one
235  first_state = list(self.fragments[rep][name].keys())[0]
236  if state is first_state:
237  representation = self._all_representations[rep]
238  representation.extend(self.fragments[rep][name][state])
239 
240  def add_fragment(self, state, representation, fragment):
241  """Add a model fragment."""
242  comp = fragment.component
243  id_rep = id(representation)
244  self._all_representations[id_rep] = representation
245  if id_rep not in self.fragments:
246  self.fragments[id_rep] = OrderedDict()
247  if comp not in self.fragments[id_rep]:
248  self.fragments[id_rep][comp] = OrderedDict()
249  if state not in self.fragments[id_rep][comp]:
250  self.fragments[id_rep][comp][state] = []
251  fragments = self.fragments[id_rep][comp][state]
252  if len(fragments) == 0 or not fragments[-1].combine(fragment):
253  fragments.append(fragment)
254  # Assume representation for a component is the same in all states,
255  # so only write out the first one
256  first_state = list(self.fragments[id_rep][comp].keys())[0]
257  if state is first_state:
258  representation.append(fragment)
259 
260 
261 class _AllDatasets(object):
262  """Track all datasets generated by PMI and add them to the ihm.System"""
263  def __init__(self, system):
264  self.system = system
265  self._datasets = []
266  self._datasets_by_state = {}
267  self._restraints_by_state = {}
268 
269  def get_all_group(self, state):
270  """Get a DatasetGroup encompassing all datasets so far in this state"""
271  # Note that if a restraint dataset is replaced after this is called,
272  # the group will still contain the old dataset - mark dataset as read
273  # only?
274  g = ihm.dataset.DatasetGroup(self._datasets_by_state.get(state, [])
275  + [r.dataset for r in self._restraints_by_state.get(state, [])
276  if r.dataset])
277  return g
278 
279  def add(self, state, dataset):
280  """Add a new dataset."""
281  self._datasets.append(dataset)
282  if state not in self._datasets_by_state:
283  self._datasets_by_state[state] = []
284  self._datasets_by_state[state].append(dataset)
285  # Make sure that the dataset ends up in the mmCIF file even if nothing
286  # refers to it
287  self.system.orphan_datasets.append(dataset)
288  return dataset
289 
290  def add_restraint(self, state, restraint):
291  """Add the dataset for a restraint"""
292  if state not in self._restraints_by_state:
293  self._restraints_by_state[state] = []
294  self._restraints_by_state[state].append(restraint)
295 
296 
297 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
298  """Restrain to a set of cross-links"""
299 
300  assembly = None
301  _label_map = {'wtDSS':'DSS', 'scDSS':'DSS', 'scEDC':'EDC'}
302  _descriptor_map = {'DSS': ihm.cross_linkers.dss,
303  'EDC': ihm.cross_linkers.edc}
304 
305  def __init__(self, pmi_restraint):
306  self.pmi_restraint = pmi_restraint
307  # Use ChemDescriptor from PMI restraint if available, otherwise guess
308  # it using the label
309  linker = getattr(self.pmi_restraint, 'linker', None)
310  label = self.pmi_restraint.label
311  self.label = label
312  super(_CrossLinkRestraint, self).__init__(
313  dataset=self.pmi_restraint.dataset,
314  linker=linker or self._get_chem_descriptor(label))
315 
316  @classmethod
317  def _get_chem_descriptor(cls, label):
318  # Map commonly-used subtypes to more standard labels
319  label = cls._label_map.get(label, label)
320  if label not in cls._descriptor_map:
321  # If label is not a standard linker type, make a new chemical
322  # descriptor containing just the name. We don't know the chemistry
323  # so cannot specify a SMILES or INCHI string for it at this point
324  d = ihm.ChemDescriptor(label)
325  cls._descriptor_map[label] = d
326  return cls._descriptor_map[label]
327 
328  def _set_psi_sigma(self, model):
329  old_values = []
330  # Do nothing if given a different state than the restraint applies to
331  if model.m != self.pmi_restraint.model:
332  return
333  for resolution in self.pmi_restraint.sigma_dictionary:
334  statname = 'ISDCrossLinkMS_Sigma_%s_%s' % (resolution, self.label)
335  if model.stats and statname in model.stats:
336  sigma = float(model.stats[statname])
337  p = self.pmi_restraint.sigma_dictionary[resolution][0]
338  old_values.append((p, p.get_scale()))
339  p.set_scale(sigma)
340  for psiindex in self.pmi_restraint.psi_dictionary:
341  statname = 'ISDCrossLinkMS_Psi_%s_%s' % (psiindex, self.label)
342  if model.stats and statname in model.stats:
343  psi = float(model.stats[statname])
344  p = self.pmi_restraint.psi_dictionary[psiindex][0]
345  old_values.append((p, p.get_scale()))
346  p.set_scale(psi)
347  # If the same particle is set multiple times, make sure we get the
348  # original value last
349  return list(reversed(old_values))
350 
351  def add_fits_from_model_statfile(self, model):
352  # Set psi/sigma for all particles
353  old_values = self._set_psi_sigma(model)
354  # If no stats were collected, we can't show the fit
355  if not old_values:
356  return
357  for xl in self.cross_links:
358  # Get current psi/sigma particle value for each XL
359  xl.fits[model] = ihm.restraint.CrossLinkFit(
360  psi=xl.psi, sigma1=xl.sigma1, sigma2=xl.sigma2)
361  # Restore original psi/sigma
362  for p, scale in old_values:
363  p.set_scale(scale)
364 
365  # Have our dataset point to that in the original PMI restraint
366  def __set_dataset(self, val):
367  self.pmi_restraint.dataset = val
368  dataset = property(lambda self: self.pmi_restraint.dataset,
369  __set_dataset)
370 
371 
372 def get_asym_mapper_for_state(simo, state, asym_map):
373  asym = asym_map.get(state, None)
374  if asym is None:
375  asym = _AsymMapper(simo, state.prot)
376  asym_map[state] = asym
377  return asym
378 
379 class _PMICrossLink(object):
380  # Query PMI particles to get psi and sigma values at time of use, not
381  # time of construction (since they may be set after creating the restraint)
382  psi = property(lambda self: self.psi_p.get_scale(),
383  lambda self, val: None)
384  sigma1 = property(lambda self: self.sigma1_p.get_scale(),
385  lambda self, val: None)
386  sigma2 = property(lambda self: self.sigma2_p.get_scale(),
387  lambda self, val: None)
388 
389 
390 class _ResidueCrossLink(ihm.restraint.ResidueCrossLink, _PMICrossLink):
391  pass
392 
393 
394 class _FeatureCrossLink(ihm.restraint.FeatureCrossLink, _PMICrossLink):
395  pass
396 
397 
398 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
399  def __init__(self, state, pmi_restraint, image_number, resolution,
400  pixel_size, image_resolution, projection_number,
401  micrographs_number):
402  self.pmi_restraint, self.image_number = pmi_restraint, image_number
403  super(_EM2DRestraint, self).__init__(
404  dataset=pmi_restraint.datasets[image_number],
405  assembly=state.modeled_assembly,
406  segment=False, number_raw_micrographs=micrographs_number,
407  pixel_size_width=pixel_size, pixel_size_height=pixel_size,
408  image_resolution=image_resolution,
409  number_of_projections=projection_number)
410 
411  # Have our dataset point to that in the original PMI restraint
412  def __get_dataset(self):
413  return self.pmi_restraint.datasets[self.image_number]
414  def __set_dataset(self, val):
415  self.pmi_restraint.datasets[self.image_number] = val
416  dataset = property(__get_dataset, __set_dataset)
417 
418  def add_fits_from_model_statfile(self, model):
419  ccc = self._get_cross_correlation(model)
420  transform = self._get_transformation(model)
421  rot = transform.get_rotation()
422  rm = [[e for e in rot.get_rotation_matrix_row(i)] for i in range(3)]
423  self.fits[model] = ihm.restraint.EM2DRestraintFit(
424  cross_correlation_coefficient=ccc,
425  rot_matrix=rm,
426  tr_vector=transform.get_translation())
427 
428  def _get_transformation(self, model):
429  """Get the transformation that places the model on the image"""
430  stats = model.em2d_stats or model.stats
431  prefix = 'ElectronMicroscopy2D_%s_Image%d' % (self.pmi_restraint.label,
432  self.image_number + 1)
433  r = [float(stats[prefix + '_Rotation%d' % i]) for i in range(4)]
434  t = [float(stats[prefix + '_Translation%d' % i])
435  for i in range(3)]
436  # If the model coordinates are transformed, need to back transform
437  # them first
438  inv = model.transform.get_inverse()
440  IMP.algebra.Vector3D(*t)) * inv
441 
442  def _get_cross_correlation(self, model):
443  """Get the cross correlation coefficient between the model projection
444  and the image"""
445  stats = model.em2d_stats or model.stats
446  return float(stats['ElectronMicroscopy2D_%s_Image%d_CCC'
447  % (self.pmi_restraint.label,
448  self.image_number + 1)])
449 
450 
451 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
452 
453  def __init__(self, simo, state, pmi_restraint, target_ps, densities):
454  self.pmi_restraint = pmi_restraint
455  super(_EM3DRestraint, self).__init__(
456  dataset=pmi_restraint.dataset,
457  assembly=self._get_assembly(densities, simo, state),
458  fitting_method='Gaussian mixture models',
459  number_of_gaussians=len(target_ps))
460 
461  # Have our dataset point to that in the original PMI restraint
462  def __set_dataset(self, val):
463  self.pmi_restraint.dataset = val
464  dataset = property(lambda self: self.pmi_restraint.dataset,
465  __set_dataset)
466 
467  def _get_assembly(self, densities, simo, state):
468  """Get the Assembly that this restraint acts on"""
469  cm = _ComponentMapper(state.prot)
470  components = {}
471  for d in densities:
472  components[cm[d]] = None # None == all residues in this component
473  a = simo._get_subassembly(components,
474  name="EM subassembly",
475  description="All components that fit the EM map")
476  return a
477 
478  def add_fits_from_model_statfile(self, model):
479  ccc = self._get_cross_correlation(model)
480  self.fits[model] = ihm.restraint.EM3DRestraintFit(
481  cross_correlation_coefficient=ccc)
482 
483  def _get_cross_correlation(self, model):
484  """Get the cross correlation coefficient between the model
485  and the map"""
486  if model.stats is not None:
487  return float(model.stats['GaussianEMRestraint_%s_CCC'
488  % self.pmi_restraint.label])
489 
490 
491 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
492 
493  def __init__(self, simo, state, pmi_restraint, geometric_object,
494  feature, distance, sigma):
495  self.pmi_restraint = pmi_restraint
496  super(_GeometricRestraint, self).__init__(
497  dataset=pmi_restraint.dataset,
498  geometric_object=geometric_object, feature=feature,
499  distance=distance, harmonic_force_constant = 1. / sigma,
500  restrain_all=True)
501 
502  # Have our dataset point to that in the original PMI restraint
503  def __set_dataset(self, val):
504  self.pmi_restraint.dataset = val
505  dataset = property(lambda self: self.pmi_restraint.dataset,
506  __set_dataset)
507 
508 
509 class _ReplicaExchangeProtocolStep(ihm.protocol.Step):
510  def __init__(self, state, rex):
511  if rex.monte_carlo_sample_objects is not None:
512  method = 'Replica exchange monte carlo'
513  else:
514  method = 'Replica exchange molecular dynamics'
515  self.monte_carlo_temperature = rex.vars['monte_carlo_temperature']
516  self.replica_exchange_minimum_temperature = \
517  rex.vars['replica_exchange_minimum_temperature']
518  self.replica_exchange_maximum_temperature = \
519  rex.vars['replica_exchange_maximum_temperature']
520  super(_ReplicaExchangeProtocolStep, self).__init__(
521  assembly=state.modeled_assembly,
522  dataset_group=None, # filled in by add_step()
523  method=method, name='Sampling',
524  num_models_begin=None, # filled in by add_step()
525  num_models_end=rex.vars["number_of_frames"],
526  multi_scale=True, multi_state=False, ordered=False)
527 
528 
529 class _ReplicaExchangeProtocolDumper(ihm.dumper.Dumper):
530  """Write IMP-specific information about replica exchange to mmCIF.
531  Note that IDs will have already been assigned by python-ihm's
532  standard modeling protocol dumper."""
533  def dump(self, system, writer):
534  with writer.loop("_imp_replica_exchange_protocol",
535  ["protocol_id", "step_id", "monte_carlo_temperature",
536  "replica_exchange_minimum_temperature",
537  "replica_exchange_maximum_temperature"]) as l:
538  for p in system._all_protocols():
539  for s in p.steps:
540  if isinstance(s, _ReplicaExchangeProtocolStep):
541  self._dump_step(p, s, l)
542 
543  def _dump_step(self, p, s, l):
544  l.write(protocol_id=p._id, step_id=s._id,
545  monte_carlo_temperature=s.monte_carlo_temperature,
546  replica_exchange_minimum_temperature= \
547  s.replica_exchange_minimum_temperature,
548  replica_exchange_maximum_temperature= \
549  s.replica_exchange_maximum_temperature)
550 
551 
552 class _ReplicaExchangeProtocolHandler(ihm.reader.Handler):
553  category = '_imp_replica_exchange_protocol'
554 
555  """Read IMP-specific information about replica exchange from mmCIF."""
556  def __call__(self, protocol_id, step_id, monte_carlo_temperature,
557  replica_exchange_minimum_temperature,
558  replica_exchange_maximum_temperature):
559  p = self.sysr.protocols.get_by_id(protocol_id)
560  # todo: match step_id properly (don't assume contiguous)
561  s = p.steps[int(step_id)-1]
562  # Upgrade from plain ihm Step to IMP subclass
563  s.__class__ = _ReplicaExchangeProtocolStep
564  s.monte_carlo_temperature = \
565  self.get_float(monte_carlo_temperature)
566  s.replica_exchange_minimum_temperature = \
567  self.get_float(replica_exchange_minimum_temperature)
568  s.replica_exchange_maximum_temperature = \
569  self.get_float(replica_exchange_maximum_temperature)
570 
571 
572 class _SimpleProtocolStep(ihm.protocol.Step):
573  def __init__(self, state, num_models_end, method):
574  super(_SimpleProtocolStep, self).__init__(
575  assembly=state.modeled_assembly,
576  dataset_group=None, # filled in by add_step()
577  method=method, name='Sampling',
578  num_models_begin=None, # filled in by add_step()
579  num_models_end=num_models_end,
580  multi_scale=True, multi_state=False, ordered=False)
581 
582 
583 class _Chain(object):
584  """Represent a single chain in a Model"""
585  def __init__(self, pmi_chain_id, asym_unit):
586  self.pmi_chain_id, self.asym_unit = pmi_chain_id, asym_unit
587  self.spheres = []
588  self.atoms = []
589 
590  def add(self, xyz, atom_type, residue_type, residue_index,
591  all_indexes, radius):
592  if atom_type is None:
593  self.spheres.append((xyz, residue_type, residue_index,
594  all_indexes, radius))
595  else:
596  self.atoms.append((xyz, atom_type, residue_type, residue_index,
597  all_indexes, radius))
598  orig_comp = property(lambda self: self.comp)
599 
600 class _TransformedChain(object):
601  """Represent a chain that is a transformed version of another"""
602  def __init__(self, orig_chain, asym_unit, transform):
603  self.orig_chain, self.asym_unit = orig_chain, asym_unit
604  self.transform = transform
605 
606  def __get_spheres(self):
607  for (xyz, residue_type, residue_index, all_indexes,
608  radius) in self.orig_chain.spheres:
609  yield (self.transform * xyz, residue_type, residue_index,
610  all_indexes, radius)
611  spheres = property(__get_spheres)
612 
613  def __get_atoms(self):
614  for (xyz, atom_type, residue_type, residue_index, all_indexes,
615  radius) in self.orig_chain.atoms:
616  yield (self.transform * xyz, atom_type, residue_type, residue_index,
617  all_indexes, radius)
618  atoms = property(__get_atoms)
619 
620  entity = property(lambda self: self.orig_chain.entity)
621  orig_comp = property(lambda self: self.orig_chain.comp)
622 
623 class _Excluder(object):
624  def __init__(self, component, simo):
625  self._seqranges = simo._exclude_coords.get(component, [])
626 
627  def is_excluded(self, indexes):
628  """Return True iff the given sequence range is excluded."""
629  for seqrange in self._seqranges:
630  if indexes[0] >= seqrange[0] and indexes[-1] <= seqrange[1]:
631  return True
632 
633 
634 class _Model(ihm.model.Model):
635  def __init__(self, prot, simo, protocol, assembly, representation):
636  super(_Model, self).__init__(assembly=assembly, protocol=protocol,
637  representation=representation)
638  self.simo = weakref.proxy(simo)
639  # Transformation from IMP coordinates into mmCIF coordinate system.
640  # Normally we pass through coordinates unchanged, but in some cases
641  # we may want to translate them (e.g. Nup84, where the deposited PDB
642  # files are all centered; we want the mmCIF files to match)
644  self.em2d_stats = None
645  self.stats = None
646  # True iff restraints act on this model
647  self._is_restrained = True
648  o = self.output = IMP.pmi.output.Output(atomistic=True)
649  name = 'cif-output'
650  self.m = prot.get_model()
651  o.dictionary_pdbs[name] = prot
652  o._init_dictchain(name, prot, multichar_chain=True)
653  (particle_infos_for_pdb,
654  self.geometric_center) = o.get_particle_infos_for_pdb_writing(name)
655  self.geometric_center = IMP.algebra.Vector3D(*self.geometric_center)
656  self._make_spheres_atoms(particle_infos_for_pdb, o, name, simo)
657  self.rmsf = {}
658 
659  def all_chains(self, simo):
660  """Yield all chains, including transformed ones"""
661  chain_for_comp = {}
662  for c in self.chains:
663  yield c
664  chain_for_comp[c.comp] = c
665  for tc in simo._transformed_components:
666  orig_chain = chain_for_comp.get(tc.original, None)
667  if orig_chain:
668  asym = simo.asym_units[tc.name]
669  c = _TransformedChain(orig_chain, asym, tc.transform)
670  c.comp = tc.name
671  yield c
672 
673  def _make_spheres_atoms(self, particle_infos_for_pdb, o, name, simo):
674  entity_for_chain = {}
675  comp_for_chain = {}
676  correct_asym = {}
677  for protname, chain_id in o.dictchain[name].items():
678  if protname in simo.entities:
679  entity_for_chain[chain_id] = simo.entities[protname]
680  else:
681  # Remove copy number
682  pn = protname.split('.')[0]
683  entity_for_chain[chain_id] = simo.entities[pn]
684  comp_for_chain[chain_id] = protname
685  # When doing multi-state modeling, the chain ID returned here
686  # (assigned sequentially) might not be correct (states may have
687  # gaps in the chain IDs). Map it to the correct asym unit.
688  correct_asym[chain_id] = simo.asym_units[protname]
689 
690  # Gather by chain ID (should be sorted by chain ID already)
691  self.chains = []
692  chain = None
693  excluder = None
694 
695  for (xyz, atom_type, residue_type, chain_id, residue_index,
696  all_indexes, radius) in particle_infos_for_pdb:
697  if chain is None or chain.pmi_chain_id != chain_id:
698  chain = _Chain(chain_id, correct_asym[chain_id])
699  chain.entity = entity_for_chain[chain_id]
700  chain.comp = comp_for_chain[chain_id]
701  self.chains.append(chain)
702  excluder = _Excluder(chain.comp, simo)
703  if not excluder.is_excluded(all_indexes if all_indexes
704  else [residue_index]):
705  chain.add(xyz, atom_type, residue_type, residue_index,
706  all_indexes, radius)
707 
708  def parse_rmsf_file(self, fname, component):
709  self.rmsf[component] = rmsf = {}
710  with open(fname) as fh:
711  for line in fh:
712  resnum, blocknum, val = line.split()
713  rmsf[int(resnum)] = (int(blocknum), float(val))
714 
715  def get_rmsf(self, component, indexes):
716  """Get the RMSF value for the given residue indexes."""
717  if not self.rmsf:
718  return None
719  rmsf = self.rmsf[component]
720  blocknums = dict.fromkeys(rmsf[ind][0] for ind in indexes)
721  if len(blocknums) != 1:
722  raise ValueError("Residue indexes %s aren't all in the same block"
723  % str(indexes))
724  return rmsf[indexes[0]][1]
725 
726  def get_atoms(self):
727  for chain in self.all_chains(self.simo):
728  pmi_offset = chain.asym_unit.entity.pmi_offset
729  for atom in chain.atoms:
730  (xyz, atom_type, residue_type, residue_index,
731  all_indexes, radius) = atom
732  pt = self.transform * xyz
733  yield ihm.model.Atom(asym_unit=chain.asym_unit,
734  seq_id=residue_index - pmi_offset,
735  atom_id=atom_type.get_string(),
736  type_symbol=None, # todo: get element
737  x=pt[0], y=pt[1], z=pt[2])
738 
739  def get_spheres(self):
740  for chain in self.all_chains(self.simo):
741  pmi_offset = chain.asym_unit.entity.pmi_offset
742  for sphere in chain.spheres:
743  (xyz, residue_type, residue_index,
744  all_indexes, radius) = sphere
745  if all_indexes is None:
746  all_indexes = (residue_index,)
747  pt = self.transform * xyz
748  yield ihm.model.Sphere(asym_unit=chain.asym_unit,
749  seq_id_range=(all_indexes[0] - pmi_offset,
750  all_indexes[-1] - pmi_offset),
751  x=pt[0], y=pt[1], z=pt[2], radius=radius,
752  rmsf=self.get_rmsf(chain.orig_comp, all_indexes))
753 
754 
755 class _AllProtocols(object):
756  def __init__(self, simo):
757  self.simo = weakref.proxy(simo)
758  # Lists of Protocols by state
759  self.protocols = OrderedDict()
760 
761  def add_protocol(self, state):
762  """Add a new Protocol"""
763  if state not in self.protocols:
764  self.protocols[state] = []
765  p = ihm.protocol.Protocol()
766  self.simo.system.orphan_protocols.append(p)
767  self.protocols[state].append(p)
768 
769  def add_step(self, step, state):
770  """Add a ProtocolStep to the last Protocol of the given State"""
771  if state not in self.protocols:
772  self.add_protocol(state)
773  protocol = self.get_last_protocol(state)
774  if len(protocol.steps) == 0:
775  step.num_models_begin = 0
776  else:
777  step.num_models_begin = protocol.steps[-1].num_models_end
778  protocol.steps.append(step)
779  step.id = len(protocol.steps)
780  # Assume that protocol uses all currently-defined datasets
781  step.dataset_group = self.simo.all_datasets.get_all_group(state)
782 
783  def add_postproc(self, step, state):
784  """Add a postprocessing step to the last protocol"""
785  protocol = self.get_last_protocol(state)
786  if len(protocol.analyses) == 0:
787  protocol.analyses.append(ihm.analysis.Analysis())
788  protocol.analyses[-1].steps.append(step)
789 
790  def get_last_protocol(self, state):
791  """Return the most recently-added _Protocol"""
792  return self.protocols[state][-1]
793 
794 
795 class _AllStartingModels(object):
796  def __init__(self, simo):
797  self.simo = simo
798  # dict of starting models (entire PDB files), collected from fragments,
799  # ordered by component name and state
800  self.models = OrderedDict()
801  self.output = IMP.pmi.output.Output()
802 
803  def add_pdb_fragment(self, fragment):
804  """Add a starting model PDB fragment."""
805  comp = fragment.component
806  state = fragment.state
807  if comp not in self.models:
808  self.models[comp] = OrderedDict()
809  if state not in self.models[comp]:
810  self.models[comp][state] = []
811  models = self.models[comp][state]
812  if len(models) == 0 \
813  or models[-1].fragments[0].pdbname != fragment.pdbname:
814  model = self._add_model(fragment)
815  models.append(model)
816  else:
817  # Break circular ref between fragment and model
818  models[-1].fragments.append(weakref.proxy(fragment))
819  # Update residue range to cover all fragments
820  pmi_offset = models[-1].asym_unit.entity.pmi_offset
821  sid_begin = min(fragment.start + fragment.offset - pmi_offset,
822  models[-1].asym_unit.seq_id_range[0])
823  sid_end = max(fragment.end + fragment.offset - pmi_offset,
824  models[-1].asym_unit.seq_id_range[1])
825  models[-1].asym_unit = fragment.asym_unit.asym(sid_begin, sid_end)
826  fragment.starting_model = models[-1]
827 
828  def _add_model(self, f):
829  parser = ihm.metadata.PDBParser()
830  r = parser.parse_file(f.pdbname)
831 
832  self.simo._add_dataset(r['dataset'])
833  # We only want the templates that model the starting model chain
834  templates = r['templates'].get(f.chain, [])
835  for t in templates:
836  if t.alignment_file:
837  self.simo.system.locations.append(t.alignment_file)
838  if t.dataset:
839  self.simo._add_dataset(t.dataset)
840  source = r['entity_source'].get(f.chain)
841  if source:
842  f.asym_unit.entity.source = source
843  pmi_offset = f.asym_unit.entity.pmi_offset
844  m = _StartingModel(
845  asym_unit=f.asym_unit.asym.pmi_range(f.start, f.end),
846  dataset=r['dataset'], asym_id=f.chain,
847  templates=templates, offset=f.offset + pmi_offset,
848  metadata=r['metadata'],
849  software=r['software'][0] if r['software'] else None,
850  script_file=r['script'])
851  m.fragments = [weakref.proxy(f)]
852  return m
853 
854 
855 class _StartingModel(ihm.startmodel.StartingModel):
856  def get_seq_dif(self):
857  return self._seq_dif
858 
859  def get_atoms(self):
860  pmi_offset = self.asym_unit.entity.pmi_offset
861  mh = IMP.mmcif.data._StartingModelAtomHandler(self.templates,
862  self.asym_unit)
863  for f in self.fragments:
864  sel = IMP.atom.Selection(f.starting_hier,
865  residue_indexes=list(range(f.start - f.offset,
866  f.end - f.offset + 1)))
867  for a in mh.get_ihm_atoms(sel.get_selected_particles(),
868  f.offset - pmi_offset):
869  yield a
870  self._seq_dif = mh._seq_dif
871 
872 
873 class _ReplicaExchangeAnalysisPostProcess(ihm.analysis.ClusterStep):
874  """Post processing using AnalysisReplicaExchange0 macro"""
875 
876  def __init__(self, rex, num_models_begin):
877  self.rex = rex
878  num_models_end = 0
879  for fname in self.get_all_stat_files():
880  with open(fname) as fh:
881  num_models_end += len(fh.readlines())
882  super(_ReplicaExchangeAnalysisPostProcess, self).__init__(
883  feature='RMSD', num_models_begin=num_models_begin,
884  num_models_end=num_models_end)
885 
886  def get_stat_file(self, cluster_num):
887  return os.path.join(self.rex._outputdir, "cluster.%d" % cluster_num,
888  'stat.out')
889 
890  def get_all_stat_files(self):
891  for i in range(self.rex._number_of_clusters):
892  yield self.get_stat_file(i)
893 
894 
895 class _ReplicaExchangeAnalysisEnsemble(ihm.model.Ensemble):
896  """Ensemble generated using AnalysisReplicaExchange0 macro"""
897 
898  num_models_deposited = None # Override base class property
899 
900  def __init__(self, pp, cluster_num, model_group, num_deposit):
901  with open(pp.get_stat_file(cluster_num)) as fh:
902  num_models = len(fh.readlines())
903  super(_ReplicaExchangeAnalysisEnsemble, self).__init__(
904  num_models=num_models,
905  model_group=model_group, post_process=pp,
906  clustering_feature=pp.feature,
907  name=model_group.name)
908  self.cluster_num = cluster_num
909  self.num_models_deposited = num_deposit
910 
911  def get_rmsf_file(self, component):
912  return os.path.join(self.post_process.rex._outputdir,
913  'cluster.%d' % self.cluster_num,
914  'rmsf.%s.dat' % component)
915 
916  def load_rmsf(self, model, component):
917  fname = self.get_rmsf_file(component)
918  if os.path.exists(fname):
919  model.parse_rmsf_file(fname, component)
920 
921  def get_localization_density_file(self, fname):
922  return os.path.join(self.post_process.rex._outputdir,
923  'cluster.%d' % self.cluster_num,
924  '%s.mrc' % fname)
925 
926  def load_localization_density(self, state, fname, select_tuple, asym_units):
927  fullpath = self.get_localization_density_file(fname)
928  if os.path.exists(fullpath):
929  details = "Localization density for %s %s" \
930  % (fname, self.model_group.name)
931  local_file = ihm.location.OutputFileLocation(fullpath,
932  details=details)
933  for s in select_tuple:
934  if isinstance(s, tuple) and len(s) == 3:
935  asym = asym_units[s[2]].pmi_range(s[0], s[1])
936  else:
937  asym = asym_units[s]
938  den = ihm.model.LocalizationDensity(file=local_file,
939  asym_unit=asym)
940  self.densities.append(den)
941 
942  def load_all_models(self, simo, state):
943  stat_fname = self.post_process.get_stat_file(self.cluster_num)
944  model_num = 0
945  with open(stat_fname) as fh:
946  stats = ast.literal_eval(fh.readline())
947  # Correct path
948  rmf_file = os.path.join(os.path.dirname(stat_fname),
949  "%d.rmf3" % model_num)
950  # todo: test with real PMI2 systems
951  if os.path.exists(rmf_file):
952  rh = RMF.open_rmf_file_read_only(rmf_file)
953  system = state._pmi_object.system
954  IMP.rmf.link_hierarchies(rh, [system.hier])
955  IMP.rmf.load_frame(fh, RMF.FrameID(0))
956  # todo: fill in other data from stat file, e.g. cross-link phi/psi
957  yield stats
958  model_num += 1
959  if model_num >= self.num_models_deposited:
960  return
961 
962  # todo: also support dRMS precision
963  def _get_precision(self):
964  precfile = os.path.join(self.post_process.rex._outputdir,
965  "precision.%d.%d.out" % (self.cluster_num,
966  self.cluster_num))
967  if not os.path.exists(precfile):
968  return ihm.unknown
969  # Fail if the precision.x.x.out file doesn't match the cluster
970  r = re.compile('All .*/cluster.%d/ average centroid distance ([\d\.]+)'
971  % self.cluster_num)
972  with open(precfile) as fh:
973  for line in fh:
974  m = r.match(line)
975  if m:
976  return float(m.group(1))
977 
978  precision = property(lambda self: self._get_precision(),
979  lambda self, val: None)
980 
981 class _SimpleEnsemble(ihm.model.Ensemble):
982  """Simple manually-created ensemble"""
983 
984  num_models_deposited = None # Override base class property
985 
986  def __init__(self, pp, model_group, num_models, drmsd,
987  num_models_deposited, ensemble_file):
988  super(_SimpleEnsemble, self).__init__(
989  model_group=model_group, post_process=pp, num_models=num_models,
990  file=ensemble_file, precision=drmsd, name=model_group.name,
991  clustering_feature='dRMSD')
992  self.num_models_deposited = num_models_deposited
993 
994  def load_localization_density(self, state, component, asym, local_file):
995  den = ihm.model.LocalizationDensity(file=local_file, asym_unit=asym)
996  self.densities.append(den)
997 
998 
999 class _CustomDNAAlphabet(object):
1000  """Custom DNA alphabet that maps A,C,G,T (rather than DA,DC,DG,DT
1001  as in python-ihm)"""
1002  _comps = dict([cc.code_canonical, cc]
1003  for cc in ihm.DNAAlphabet._comps.values())
1004 
1005 
1006 class _EntityMapper(dict):
1007  """Handle mapping from IMP components (without copy number) to CIF entities.
1008  Multiple components may map to the same entity if they share sequence."""
1009  def __init__(self, system):
1010  super(_EntityMapper, self).__init__()
1011  self._sequence_dict = {}
1012  self._entities = []
1013  self.system = system
1014 
1015  def _get_alphabet(self, alphabet):
1016  """Map a PMI alphabet to an IHM alphabet"""
1017  # todo: we should probably ignore alphabets entirely and use
1018  # residue types directly (e.g. we shouldn't compare one-letter
1019  # code sequence to determine if an Entity is unique)
1020  alphabet_map = {None: ihm.LPeptideAlphabet,
1021  IMP.pmi.alphabets.amino_acid: ihm.LPeptideAlphabet,
1022  IMP.pmi.alphabets.rna: ihm.RNAAlphabet,
1023  IMP.pmi.alphabets.dna: _CustomDNAAlphabet}
1024  if alphabet in alphabet_map:
1025  return alphabet_map[alphabet]
1026  else:
1027  raise TypeError("Don't know how to handle %s" % alphabet)
1028 
1029  def add(self, component_name, sequence, offset, alphabet):
1030  def entity_seq(sequence):
1031  # Map X to UNK
1032  if 'X' in sequence:
1033  return ['UNK' if s == 'X' else s for s in sequence]
1034  else:
1035  return sequence
1036  if sequence not in self._sequence_dict:
1037  # Use the name of the first component, stripped of any copy number,
1038  # as the description of the entity
1039  d = component_name.split("@")[0].split(".")[0]
1040  entity = Entity(entity_seq(sequence), description=d,
1041  pmi_offset=offset,
1042  alphabet=self._get_alphabet(alphabet))
1043  self.system.entities.append(entity)
1044  self._sequence_dict[sequence] = entity
1045  self[component_name] = self._sequence_dict[sequence]
1046 
1047 
1048 class _TransformedComponent(object):
1049  def __init__(self, name, original, transform):
1050  self.name, self.original, self.transform = name, original, transform
1051 
1052 
1053 class _SimpleRef(object):
1054  """Class with similar interface to weakref.ref, but keeps a strong ref"""
1055  def __init__(self, ref):
1056  self.ref = ref
1057  def __call__(self):
1058  return self.ref
1059 
1060 
1061 class _State(ihm.model.State):
1062  """Representation of a single state in the system."""
1063 
1064  def __init__(self, pmi_object, po):
1065  # Point to the PMI object for this state. Use a weak reference
1066  # since the state object typically points to us too, so we need
1067  # to break the reference cycle. In PMI1 this will be a
1068  # Representation object; in PMI2 it is the PMI2 State object itself.
1069  self._pmi_object = weakref.proxy(pmi_object)
1070  if hasattr(pmi_object, 'state'):
1071  # Need a strong ref to pmi_object.state to prevent it from being
1072  # cleaned up when doing PMI1 multistate modeling
1073  self._pmi_state = _SimpleRef(pmi_object.state)
1074  else:
1075  self._pmi_state = weakref.ref(pmi_object)
1076  # Preserve PMI state name
1077  old_name = self.name
1078  super(_State, self).__init__(experiment_type='Fraction of bulk')
1079  self.name = old_name
1080 
1081  # The assembly of all components modeled by IMP in this state.
1082  # This may be smaller than the complete assembly.
1083  self.modeled_assembly = ihm.Assembly(
1084  name="Modeled assembly",
1085  description=self.get_postfixed_name(
1086  "All components modeled by IMP"))
1087  po.system.orphan_assemblies.append(self.modeled_assembly)
1088 
1089  self.all_modeled_components = []
1090 
1091  def __hash__(self):
1092  return hash(self._pmi_state())
1093  def __eq__(self, other):
1094  return self._pmi_state() == other._pmi_state()
1095 
1096  def add_model_group(self, group):
1097  self.append(group)
1098  return group
1099 
1100  def get_prefixed_name(self, name):
1101  """Prefix the given name with the state name, if available."""
1102  if self.short_name:
1103  return self.short_name + ' ' + name
1104  else:
1105  # Can't use capitalize() since that un-capitalizes everything
1106  # but the first letter
1107  return name[0].upper() + name[1:] if name else ''
1108 
1109  def get_postfixed_name(self, name):
1110  """Postfix the given name with the state name, if available."""
1111  if self.short_name:
1112  return "%s in state %s" % (name, self.short_name)
1113  else:
1114  return name
1115 
1116  short_name = property(lambda self: self._pmi_state().short_name)
1117  long_name = property(lambda self: self._pmi_state().long_name)
1118  def __get_name(self):
1119  return self._pmi_state().long_name
1120  def __set_name(self, val):
1121  self._pmi_state().long_name = val
1122  name = property(__get_name, __set_name)
1123 
1124 
1125 class Entity(ihm.Entity):
1126  """A single entity in the system.
1127  This functions identically to the base ihm.Entity class, but it
1128  allows identifying residues by either the IHM numbering scheme
1129  (seq_id, which is always contiguous starting at 1) or by the PMI
1130  scheme (which is similar, but may not start at 1 if the sequence in
1131  the FASTA file has one or more N-terminal gaps"""
1132  def __init__(self, sequence, pmi_offset, *args, **keys):
1133  # Offset between PMI numbering and IHM; <pmi_#> = <ihm_#> + pmi_offset
1134  # (pmi_offset is also the number of N-terminal gaps in the FASTA file)
1135  self.pmi_offset = pmi_offset
1136  super(Entity, self).__init__(sequence, *args, **keys)
1137 
1138  def pmi_residue(self, res_id):
1139  """Return a single IHM residue indexed using PMI numbering"""
1140  return self.residue(res_id - self.pmi_offset)
1141 
1142  def pmi_range(self, res_id_begin, res_id_end):
1143  """Return a range of IHM residues indexed using PMI numbering"""
1144  off = self.pmi_offset
1145  return self(res_id_begin - off, res_id_end - off)
1146 
1147 
1148 class AsymUnit(ihm.AsymUnit):
1149  """A single asymmetric unit in the system.
1150  This functions identically to the base ihm.AsymUnit class, but it
1151  allows identifying residues by either the IHM numbering scheme
1152  (seq_id, which is always contiguous starting at 1) or by the PMI
1153  scheme (which is similar, but may not start at 1 if the sequence in
1154  the FASTA file has one or more N-terminal gaps"""
1155 
1156  def __init__(self, entity, *args, **keys):
1157  super(AsymUnit, self).__init__(
1158  entity, auth_seq_id_map=entity.pmi_offset, *args, **keys)
1159 
1160  def pmi_residue(self, res_id):
1161  """Return a single IHM residue indexed using PMI numbering"""
1162  return self.residue(res_id - self.entity.pmi_offset)
1163 
1164  def pmi_range(self, res_id_begin, res_id_end):
1165  """Return a range of IHM residues indexed using PMI numbering"""
1166  off = self.entity.pmi_offset
1167  return self(res_id_begin - off, res_id_end - off)
1168 
1169 
1171  """Class to encode a modeling protocol as mmCIF.
1172 
1173  IMP has basic support for writing out files in mmCIF format, for
1174  deposition in [PDB-Dev](https://pdb-dev.wwpdb.org/).
1175  After creating an instance of this class, attach it to an
1176  IMP.pmi.topology.System object. After this, any
1177  generated models and metadata are automatically collected in the
1178  `system` attribute, which is an
1179  [ihm.System](https://python-ihm.readthedocs.io/en/latest/main.html#ihm.System) object.
1180  Once the protocol is complete, call finalize() to make sure `system`
1181  contains everything, then use the
1182  [python-ihm API](https://python-ihm.readthedocs.io/en/latest/dumper.html#ihm.dumper.write)
1183  to write out files in mmCIF or BinaryCIF format.
1184 
1185  See also get_handlers(), get_dumpers().
1186  """
1187  def __init__(self, fh=None):
1188  # Ultimately, collect data in an ihm.System object
1189  self.system = ihm.System()
1190  self._state_group = ihm.model.StateGroup()
1191  self.system.state_groups.append(self._state_group)
1192 
1193  self.fh = fh
1194  if fh:
1196  "Passing a file handle to ProtocolOutput is deprecated. "
1197  "Create it with no arguments, then use the python-ihm API to "
1198  "output files (after calling finalize())")
1199  self._state_ensemble_offset = 0
1200  self._main_script = os.path.abspath(sys.argv[0])
1201 
1202  # Point to the main modeling script
1203  loc = ihm.location.WorkflowFileLocation(path=self._main_script,
1204  details="The main integrative modeling script")
1205  self.system.locations.append(loc)
1206 
1207  self._states = {}
1208  self.__asym_states = {}
1209  self._working_directory = os.getcwd()
1210  self.default_representation = self.create_representation(
1211  "Default representation")
1212  self.entities = _EntityMapper(self.system)
1213  # Mapping from component names to ihm.AsymUnit
1214  self.asym_units = {}
1215  self._all_components = {}
1216  self.all_modeled_components = []
1217  self._transformed_components = []
1218  self.sequence_dict = {}
1219 
1220  # Common geometric objects used in PMI systems
1221  self._xy_plane = ihm.geometry.XYPlane()
1222  self._xz_plane = ihm.geometry.XZPlane()
1223  self._z_axis = ihm.geometry.ZAxis()
1224  self._center_origin = ihm.geometry.Center(0,0,0)
1225  self._identity_transform = ihm.geometry.Transformation.identity()
1226 
1227  # Coordinates to exclude
1228  self._exclude_coords = {}
1229 
1230  self.all_representations = _AllModelRepresentations(self)
1231  self.all_protocols = _AllProtocols(self)
1232  self.all_datasets = _AllDatasets(self.system)
1233  self.all_starting_models = _AllStartingModels(self)
1234 
1235  self.all_software = _AllSoftware(self.system)
1236 
1237  def create_representation(self, name):
1238  """Create a new Representation and return it. This can be
1239  passed to add_model(), add_bead_element() or add_pdb_element()."""
1240  r = ihm.representation.Representation(name=name)
1241  self.system.orphan_representations.append(r)
1242  return r
1243 
1244  def exclude_coordinates(self, component, seqrange):
1245  """Don't record coordinates for the given domain.
1246  Coordinates for the given domain (specified by a component name
1247  and a 2-element tuple giving the start and end residue numbers)
1248  will be excluded from the mmCIF file. This can be used to exclude
1249  parts of the structure that weren't well resolved in modeling.
1250  Any bead or residue that lies wholly within this range will be
1251  excluded. Multiple ranges for a given component can be excluded
1252  by calling this method multiple times."""
1253  if component not in self._exclude_coords:
1254  self._exclude_coords[component] = []
1255  self._exclude_coords[component].append(seqrange)
1256 
1257  def _is_excluded(self, component, start, end):
1258  """Return True iff this chunk of sequence should be excluded"""
1259  for seqrange in self._exclude_coords.get(component, ()):
1260  if start >= seqrange[0] and end <= seqrange[1]:
1261  return True
1262 
1263  def _add_state(self, state):
1264  """Create a new state and return a pointer to it."""
1265  self._state_ensemble_offset = len(self.system.ensembles)
1266  s = _State(state, self)
1267  self._state_group.append(s)
1268  self._last_state = s
1269  return s
1270 
1271  def _get_chain_for_component(self, name, output):
1272  """Get the chain ID for a component, if any."""
1273  # todo: handle multiple copies
1274  if name in self.asym_units:
1275  return self.asym_units[name]._id
1276  else:
1277  # A non-modeled component doesn't have a chain ID
1278  return None
1279 
1280  def _get_assembly_comps(self, assembly):
1281  """Get the names of the components in the given assembly"""
1282  # component name = asym_unit.details
1283  comps = {}
1284  for ca in assembly:
1285  comps[ca.details] = None
1286  return comps
1287 
1288  def create_transformed_component(self, state, name, original, transform):
1289  """Make a new component that's a transformed copy of another.
1290  All representation for the existing component is copied to the
1291  new one."""
1292  assembly_comps = self._get_assembly_comps(state.modeled_assembly)
1293  if name in assembly_comps:
1294  raise ValueError("Component %s already exists" % name)
1295  elif original not in assembly_comps:
1296  raise ValueError("Original component %s does not exist" % original)
1297  self.create_component(state, name, True)
1298  self.add_component_sequence(state, name, self.sequence_dict[original])
1299  self._transformed_components.append(_TransformedComponent(
1300  name, original, transform))
1301  self.all_representations.copy_component(state, name, original,
1302  self.asym_units[name])
1303 
1304  def create_component(self, state, name, modeled, asym_name=None):
1305  if asym_name is None:
1306  asym_name = name
1307  new_comp = name not in self._all_components
1308  self._all_components[name] = None
1309  if modeled:
1310  state.all_modeled_components.append(name)
1311  if asym_name not in self.asym_units:
1312  # assign asym once we get sequence
1313  self.asym_units[asym_name] = None
1314  if new_comp:
1315  self.all_modeled_components.append(name)
1316 
1317  def add_component_sequence(self, state, name, seq, asym_name=None,
1318  alphabet=None):
1319  if asym_name is None:
1320  asym_name = name
1321  def get_offset(seq):
1322  # Count length of gaps at start of sequence
1323  for i in range(len(seq)):
1324  if seq[i] != '-':
1325  return seq[i:], i
1326  seq, offset = get_offset(seq)
1327  if name in self.sequence_dict:
1328  if self.sequence_dict[name] != seq:
1329  raise ValueError("Sequence mismatch for component %s" % name)
1330  else:
1331  self.sequence_dict[name] = seq
1332  self.entities.add(name, seq, offset, alphabet)
1333  if asym_name in self.asym_units:
1334  if self.asym_units[asym_name] is None:
1335  # Set up a new asymmetric unit for this component
1336  entity = self.entities[name]
1337  asym = AsymUnit(entity, details=asym_name)
1338  self.system.asym_units.append(asym)
1339  self.asym_units[asym_name] = asym
1340  state.modeled_assembly.append(self.asym_units[asym_name])
1341 
1342  def _add_restraint_model_fits(self):
1343  """Add fits to restraints for all known models"""
1344  for group, m in self.system._all_models():
1345  if m._is_restrained:
1346  for r in self.system.restraints:
1347  if hasattr(r, 'add_fits_from_model_statfile'):
1348  r.add_fits_from_model_statfile(m)
1349 
1351  "Use finalize() then the python-ihm API instead")
1352  def flush(self, format='mmCIF'):
1353  """Write out all information to the file.
1354  Information can be written in any format supported by
1355  the ihm library (typically this is 'mmCIF' or 'BCIF')."""
1356  self.finalize()
1357  ihm.dumper.write(self.fh, [self.system], format, dumpers=get_dumpers())
1358 
1359  def finalize(self):
1360  """Do any final processing on the class hierarchy.
1361  After calling this method, the `system` member (an instance
1362  of `ihm.System`) completely reproduces the PMI modeling, and
1363  can be written out to an mmCIF file with `ihm.dumper.write`,
1364  and/or modified using the ihm API."""
1365  self._add_restraint_model_fits()
1366 
1367  def add_pdb_element(self, state, name, start, end, offset, pdbname,
1368  chain, hier, representation=None):
1369  if self._is_excluded(name, start, end):
1370  return
1371  if representation is None:
1372  representation = self.default_representation
1373  asym = self.asym_units[name]
1374  p = _PDBFragment(state, name, start, end, offset, pdbname, chain,
1375  hier, asym)
1376  self.all_representations.add_fragment(state, representation, p)
1377  self.all_starting_models.add_pdb_fragment(p)
1378 
1379  def add_bead_element(self, state, name, start, end, num, hier,
1380  representation=None):
1381  if self._is_excluded(name, start, end):
1382  return
1383  if representation is None:
1384  representation = self.default_representation
1385  asym = self.asym_units[name]
1386  pmi_offset = asym.entity.pmi_offset
1387  b = _BeadsFragment(state, name, start - pmi_offset, end - pmi_offset,
1388  num, hier, asym)
1389  self.all_representations.add_fragment(state, representation, b)
1390 
1391  def get_cross_link_group(self, pmi_restraint):
1392  r = _CrossLinkRestraint(pmi_restraint)
1393  self.system.restraints.append(r)
1394  self._add_restraint_dataset(r) # so that all-dataset group works
1395  return r
1396 
1397  def add_experimental_cross_link(self, r1, c1, r2, c2, rsr):
1398  if c1 not in self._all_components or c2 not in self._all_components:
1399  # Crosslink refers to a component we didn't model
1400  # As a quick hack, just ignore it.
1401  # todo: need to add an entity for this anyway (so will need the
1402  # sequence, and add to struct_assembly)
1403  return None
1404  e1 = self.entities[c1]
1405  e2 = self.entities[c2]
1406  xl = ihm.restraint.ExperimentalCrossLink(residue1=e1.pmi_residue(r1),
1407  residue2=e2.pmi_residue(r2))
1408  rsr.experimental_cross_links.append([xl])
1409  return xl
1410 
1411  def add_cross_link(self, state, ex_xl, p1, p2, length, sigma1_p, sigma2_p,
1412  psi_p, rsr):
1413  # Map p1/p2 to asym units
1414  asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1415  d = ihm.restraint.UpperBoundDistanceRestraint(length)
1416 
1417  if _get_by_residue(p1) and _get_by_residue(p2):
1418  cls = _ResidueCrossLink
1419  else:
1420  cls = _FeatureCrossLink
1421  xl = cls(ex_xl, asym1=asym[p1], asym2=asym[p2], distance=d,
1422  restrain_all=True)
1423  # Needed to look up fits later
1424  xl.psi_p, xl.sigma1_p, xl.sigma2_p = psi_p, sigma1_p, sigma2_p
1425  rsr.cross_links.append(xl)
1426 
1427  def add_replica_exchange(self, state, rex):
1428  # todo: allow for metadata to say how many replicas were used in the
1429  # actual experiment, and how many independent runs were carried out
1430  # (use these as multipliers to get the correct total number of
1431  # output models)
1432  step = _ReplicaExchangeProtocolStep(state, rex)
1433  step.software = self.all_software.pmi
1434  self.all_protocols.add_step(step, state)
1435 
1436  def _add_simple_dynamics(self, num_models_end, method):
1437  # Always assumed that we're dealing with the last state
1438  state = self._last_state
1439  self.all_protocols.add_step(_SimpleProtocolStep(state, num_models_end,
1440  method), state)
1441 
1442  def _add_protocol(self):
1443  # Always assumed that we're dealing with the last state
1444  state = self._last_state
1445  self.all_protocols.add_protocol(state)
1446 
1447  def _add_dataset(self, dataset):
1448  return self.all_datasets.add(self._last_state, dataset)
1449 
1450  def _add_restraint_dataset(self, restraint):
1451  return self.all_datasets.add_restraint(self._last_state, restraint)
1452 
1453  def _add_simple_postprocessing(self, num_models_begin, num_models_end):
1454  # Always assumed that we're dealing with the last state
1455  state = self._last_state
1456  pp = ihm.analysis.ClusterStep('RMSD', num_models_begin, num_models_end)
1457  self.all_protocols.add_postproc(pp, state)
1458  return pp
1459 
1460  def _add_no_postprocessing(self, num_models):
1461  # Always assumed that we're dealing with the last state
1462  state = self._last_state
1463  pp = ihm.analysis.EmptyStep()
1464  pp.num_models_begin = pp.num_models_end = num_models
1465  self.all_protocols.add_postproc(pp, state)
1466  return pp
1467 
1468  def _add_simple_ensemble(self, pp, name, num_models, drmsd,
1469  num_models_deposited, localization_densities,
1470  ensemble_file):
1471  """Add an ensemble generated by ad hoc methods (not using PMI).
1472  This is currently only used by the Nup84 system."""
1473  # Always assumed that we're dealing with the last state
1474  state = self._last_state
1475  group = ihm.model.ModelGroup(name=state.get_postfixed_name(name))
1476  state.add_model_group(group)
1477  if ensemble_file:
1478  self.system.locations.append(ensemble_file)
1479  e = _SimpleEnsemble(pp, group, num_models, drmsd, num_models_deposited,
1480  ensemble_file)
1481  self.system.ensembles.append(e)
1482  for c in state.all_modeled_components:
1483  den = localization_densities.get(c, None)
1484  if den:
1485  e.load_localization_density(state, c, self.asym_units[c], den)
1486  return e
1487 
1488  def set_ensemble_file(self, i, location):
1489  """Point a previously-created ensemble to an 'all-models' file.
1490  This could be a trajectory such as DCD, an RMF, or a multimodel
1491  PDB file."""
1492  self.system.locations.append(location)
1493  # Ensure that we point to an ensemble related to the current state
1494  ind = i + self._state_ensemble_offset
1495  self.system.ensembles[ind].file = location
1496 
1497  def add_replica_exchange_analysis(self, state, rex, density_custom_ranges):
1498  # todo: add prefilter as an additional postprocess step (complication:
1499  # we don't know how many models it filtered)
1500  # todo: postpone rmsf/density localization extraction until after
1501  # relevant methods are called (currently it is done from the
1502  # constructor)
1503  protocol = self.all_protocols.get_last_protocol(state)
1504  num_models = protocol.steps[-1].num_models_end
1505  pp = _ReplicaExchangeAnalysisPostProcess(rex, num_models)
1506  pp.software = self.all_software.pmi
1507  self.all_protocols.add_postproc(pp, state)
1508  for i in range(rex._number_of_clusters):
1509  group = ihm.model.ModelGroup(name=state.get_prefixed_name(
1510  'cluster %d' % (i + 1)))
1511  state.add_model_group(group)
1512  # todo: make # of models to deposit configurable somewhere
1513  e = _ReplicaExchangeAnalysisEnsemble(pp, i, group, 1)
1514  self.system.ensembles.append(e)
1515  # Add localization density info if available
1516  for fname, stuple in sorted(density_custom_ranges.items()):
1517  e.load_localization_density(state, fname, stuple,
1518  self.asym_units)
1519  for stats in e.load_all_models(self, state):
1520  m = self.add_model(group)
1521  # Since we currently only deposit 1 model, it is the
1522  # best scoring one
1523  m.name = 'Best scoring model'
1524  m.stats = stats
1525  # Add RMSF info if available
1526  for c in state.all_modeled_components:
1527  e.load_rmsf(m, c)
1528 
1529  def _get_subassembly(self, comps, name, description):
1530  """Get an Assembly consisting of the given components.
1531  `compdict` is a dictionary of the components to add, where keys
1532  are the component names and values are the sequence ranges (or
1533  None to use all residues in the component)."""
1534  asyms = []
1535  for comp, seqrng in comps.items():
1536  a = self.asym_units[comp]
1537  asyms.append(a if seqrng is None else a(*seqrng))
1538 
1539  a = ihm.Assembly(asyms, name=name, description=description)
1540  return a
1541 
1542  def _add_foxs_restraint(self, model, comp, seqrange, dataset, rg, chi,
1543  details):
1544  """Add a basic FoXS fit. This is largely intended for use from the
1545  NPC application."""
1546  assembly = self._get_subassembly({comp:seqrange},
1547  name="SAXS subassembly",
1548  description="All components that fit SAXS data")
1549  r = ihm.restraint.SASRestraint(dataset, assembly, segment=False,
1550  fitting_method='FoXS', fitting_atom_type='Heavy atoms',
1551  multi_state=False, radius_of_gyration=rg, details=details)
1552  r.fits[model] = ihm.restraint.SASRestraintFit(chi_value=chi)
1553  self.system.restraints.append(r)
1554  self._add_restraint_dataset(r) # so that all-dataset group works
1555 
1556  def add_em2d_restraint(self, state, r, i, resolution, pixel_size,
1557  image_resolution, projection_number,
1558  micrographs_number):
1559  r = _EM2DRestraint(state, r, i, resolution, pixel_size,
1560  image_resolution, projection_number,
1561  micrographs_number)
1562  self.system.restraints.append(r)
1563  self._add_restraint_dataset(r) # so that all-dataset group works
1564 
1565  def add_em3d_restraint(self, state, target_ps, densities, pmi_restraint):
1566  # todo: need to set allow_duplicates on this dataset?
1567  r = _EM3DRestraint(self, state, pmi_restraint, target_ps, densities)
1568  self.system.restraints.append(r)
1569  self._add_restraint_dataset(r) # so that all-dataset group works
1570 
1571  def add_zaxial_restraint(self, state, ps, lower_bound, upper_bound,
1572  sigma, pmi_restraint):
1573  self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1574  sigma, pmi_restraint, self._xy_plane)
1575 
1576  def add_yaxial_restraint(self, state, ps, lower_bound, upper_bound,
1577  sigma, pmi_restraint):
1578  self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1579  sigma, pmi_restraint, self._xz_plane)
1580 
1581  def add_xyradial_restraint(self, state, ps, lower_bound, upper_bound,
1582  sigma, pmi_restraint):
1583  self._add_geometric_restraint(state, ps, lower_bound, upper_bound,
1584  sigma, pmi_restraint, self._z_axis)
1585 
1586  def _add_geometric_restraint(self, state, ps, lower_bound, upper_bound,
1587  sigma, pmi_restraint, geom):
1588  asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1589  r = _GeometricRestraint(self, state, pmi_restraint, geom,
1590  asym.get_feature(ps),
1591  ihm.restraint.LowerUpperBoundDistanceRestraint(
1592  lower_bound, upper_bound),
1593  sigma)
1594  self.system.restraints.append(r)
1595  self._add_restraint_dataset(r) # so that all-dataset group works
1596 
1597  def _get_membrane(self, tor_R, tor_r, tor_th):
1598  """Get an object representing a half-torus membrane"""
1599  if not hasattr(self, '_seen_membranes'):
1600  self._seen_membranes = {}
1601  # If a membrane already exists with all dimensions within 0.01 of
1602  # this one, reuse it
1603  membrane_id = tuple(int(x * 100.) for x in (tor_R, tor_r, tor_th))
1604  if membrane_id not in self._seen_membranes:
1605  m = ihm.geometry.HalfTorus(center=self._center_origin,
1606  transformation=self._identity_transform,
1607  major_radius=tor_R, minor_radius=tor_r, thickness=tor_th,
1608  inner=True, name='Membrane')
1609  self._seen_membranes[membrane_id] = m
1610  return self._seen_membranes[membrane_id]
1611 
1612  def add_membrane_surface_location_restraint(
1613  self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1614  self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1615  pmi_restraint, ihm.restraint.UpperBoundDistanceRestraint(0.))
1616 
1617  def add_membrane_exclusion_restraint(
1618  self, state, ps, tor_R, tor_r, tor_th, sigma, pmi_restraint):
1619  self._add_membrane_restraint(state, ps, tor_R, tor_r, tor_th, sigma,
1620  pmi_restraint, ihm.restraint.LowerBoundDistanceRestraint(0.))
1621 
1622  def _add_membrane_restraint(self, state, ps, tor_R, tor_r, tor_th,
1623  sigma, pmi_restraint, rsr):
1624  asym = get_asym_mapper_for_state(self, state, self.__asym_states)
1625  r = _GeometricRestraint(self, state, pmi_restraint,
1626  self._get_membrane(tor_R, tor_r, tor_th),
1627  asym.get_feature(ps), rsr, sigma)
1628  self.system.restraints.append(r)
1629  self._add_restraint_dataset(r) # so that all-dataset group works
1630 
1631  def add_model(self, group, assembly=None, representation=None):
1632  state = self._last_state
1633  if representation is None:
1634  representation = self.default_representation
1635  protocol = self.all_protocols.get_last_protocol(state)
1636  m = _Model(state.prot, self, protocol,
1637  assembly if assembly else state.modeled_assembly,
1638  representation)
1639  group.append(m)
1640  return m
1641 
1642 
1644  """Get custom python-ihm dumpers for writing PMI to from mmCIF.
1645  This returns a list of custom dumpers that can be passed as all or
1646  part of the `dumpers` argument to ihm.dumper.write(). They add
1647  PMI-specific information to mmCIF or BinaryCIF files written out
1648  by python-ihm."""
1649  return [_ReplicaExchangeProtocolDumper]
1650 
1651 
1653  """Get custom python-ihm handlers for reading PMI data from mmCIF.
1654  This returns a list of custom handlers that can be passed as all or
1655  part of the `handlers` argument to ihm.reader.read(). They read
1656  PMI-specific information from mmCIF or BinaryCIF files read in
1657  by python-ihm."""
1658  return [_ReplicaExchangeProtocolHandler]
1659 
1660 
1662  "Use ihm.reader.read() instead with handlers=get_handlers()")
1663 def read(fh, model_class=ihm.model.Model, format='mmCIF', handlers=[]):
1664  """Read data from the mmCIF file handle `fh`.
1665  This is a simple wrapper around `ihm.reader.read()`, which also
1666  reads PMI-specific information from the mmCIF or BinaryCIF file."""
1667  return ihm.reader.read(fh, model_class=model_class, format=format,
1668  handlers=get_handlers() + handlers)
1669 
1670 
1671 class GMMParser(ihm.metadata.Parser):
1672  """Extract metadata from an EM density GMM file."""
1673 
1674  def parse_file(self, filename):
1675  """Extract metadata from `filename`.
1676  @return a dict with key `dataset` pointing to the GMM file and
1677  `number_of_gaussians` to the number of GMMs (or None)"""
1678  l = ihm.location.InputFileLocation(filename,
1679  details="Electron microscopy density map, "
1680  "represented as a Gaussian Mixture Model (GMM)")
1681  # A 3DEM restraint's dataset ID uniquely defines the mmCIF restraint, so
1682  # we need to allow duplicates
1683  l._allow_duplicates = True
1684  d = ihm.dataset.EMDensityDataset(l)
1685  ret = {'dataset':d, 'number_of_gaussians':None}
1686 
1687  with open(filename) as fh:
1688  for line in fh:
1689  if line.startswith('# data_fn: '):
1690  p = ihm.metadata.MRCParser()
1691  fn = line[11:].rstrip('\r\n')
1692  dataset = p.parse_file(os.path.join(
1693  os.path.dirname(filename), fn))['dataset']
1694  ret['dataset'].parents.append(dataset)
1695  elif line.startswith('# ncenters: '):
1696  ret['number_of_gaussians'] = int(line[12:])
1697  return ret
Select non water and non hydrogen atoms.
Definition: pdb.h:245
def get_handlers
Get custom python-ihm handlers for reading PMI data from mmCIF.
Definition: mmcif.py:1652
def create_representation
Create a new Representation and return it.
Definition: mmcif.py:1237
def create_transformed_component
Make a new component that's a transformed copy of another.
Definition: mmcif.py:1288
Simple 3D transformation class.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:156
def exclude_coordinates
Don't record coordinates for the given domain.
Definition: mmcif.py:1247
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:1183
def flush
Write out all information to the file.
Definition: mmcif.py:1350
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: atom/Atom.h:241
def parse_file
Extract metadata from filename.
Definition: mmcif.py:1675
def pmi_range
Return a range of IHM residues indexed using PMI numbering.
Definition: mmcif.py:1142
Miscellaneous utilities.
Definition: tools.py:1
void handle_use_deprecated(std::string message)
def deprecated_function
Python decorator to mark a function as deprecated.
Definition: __init__.py:9880
Extract metadata from an EM density GMM file.
Definition: mmcif.py:1671
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:731
def pmi_residue
Return a single IHM residue indexed using PMI numbering.
Definition: mmcif.py:1160
void read_pdb(TextInput input, int model, Hierarchy h)
A single asymmetric unit in the system.
Definition: mmcif.py:1152
A single entity in the system.
Definition: mmcif.py:1125
def set_ensemble_file
Point a previously-created ensemble to an 'all-models' file.
Definition: mmcif.py:1488
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:1164
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 deprecated_method
Python decorator to mark a method as deprecated.
Definition: __init__.py:9868
def finalize
Do any final processing on the class hierarchy.
Definition: mmcif.py:1359
Base class for capturing a modeling protocol.
Definition: output.py:44
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:1643
Class for easy writing of PDBs, RMFs, and stat files.
Definition: output.py:198
3D rotation class.
Definition: Rotation3D.h:51
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:135
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:1138
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
def read
Read data from the mmCIF file handle fh.
Definition: mmcif.py:1661
VectorD< 3 > Vector3D
Definition: VectorD.h:421
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:66
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:189
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...