IMP logo
IMP Reference Guide  2.20.0
The Integrative Modeling Platform
mmcif/restraint.py
1 """@namespace IMP.mmcif.restraint
2  @brief Map IMP restraints to mmCIF categories
3 """
4 
6 import ihm.restraint
7 import operator
8 
9 
10 def _get_by_residue(m, p):
11  """Determine whether the given particle represents a specific residue
12  or a more coarse-grained object."""
13  return (IMP.atom.Residue.get_is_setup(m, p)
15 
16 
17 def _get_scale(m, p):
18  """Get the numerical value of the Scale particle"""
19  if not IMP.isd.Scale.get_is_setup(m, p):
20  raise ValueError("not scale particle")
21  return IMP.isd.Scale(m, p).get_scale()
22 
23 
24 class _AsymMapper(object):
25  """Map ParticleIndexes to ihm.AsymUnit"""
26  def __init__(self, m, components, ignore_non_structural=False):
27  self.m = m
28  self.components = components
29  self._seen_ranges = {}
30  self._ignore_non_structural = ignore_non_structural
31 
32  def __getitem__(self, pi):
33  m = self.m
34  # Walk up the hierarchy until we find the Chain,
35  # then map that to AsymUnit
37  h = IMP.atom.Hierarchy(m, pi)
38  while h:
40  return self.components[IMP.atom.Chain(h)].asym_unit
41  h = h.get_parent()
42  if not self._ignore_non_structural:
43  raise KeyError("Could not find top-level Chain for "
44  + m.get_particle_name(pi))
45 
46  def get_feature(self, ps):
47  """Get an ihm.restraint.Feature that covers the given particles"""
48  # todo: handle things other than residues
49  m = self.m
50  rngs = []
51  for p in ps:
52  asym = self[p]
53  # todo: handle overlapping ranges
55  rind = IMP.atom.Residue(m, p).get_index()
56  rng = asym(rind, rind)
58  # PMI Fragments always contain contiguous residues
59  rinds = IMP.atom.Fragment(m, p).get_residue_indexes()
60  rng = asym(rinds[0], rinds[-1])
61  else:
62  raise ValueError("Unsupported particle type %s" % str(p))
63  # Join contiguous ranges
64  if len(rngs) > 0 and rngs[-1].asym == asym \
65  and rngs[-1].seq_id_range[1] == rng.seq_id_range[0] - 1:
66  rngs[-1].seq_id_range = (rngs[-1].seq_id_range[0],
67  rng.seq_id_range[1])
68  else:
69  rngs.append(rng)
70  # If an identical feature already exists, return that
71  # todo: python-ihm should handle this automatically for us
72  hrngs = tuple(rngs)
73  if hrngs in self._seen_ranges:
74  return self._seen_ranges[hrngs]
75  else:
76  feat = ihm.restraint.ResidueFeature(rngs)
77  self._seen_ranges[hrngs] = feat
78  return feat
79 
80 
81 def _parse_restraint_info(info):
82  """Convert RestraintInfo object to Python dict"""
83  d = {}
84  if info is None:
85  return d
86  info.set_was_used(True)
87  for typ in ('int', 'float', 'string', 'filename', 'floats', 'filenames',
88  'particle_indexes'):
89  for i in range(getattr(info, 'get_number_of_' + typ)()):
90  key = getattr(info, 'get_%s_key' % typ)(i)
91  value = getattr(info, 'get_%s_value' % typ)(i)
92  d[key] = value
93  return d
94 
95 
96 def _get_restraint_assembly(imp_restraint, components):
97  """Get the assembly corresponding to all input particles for
98  the restraint"""
99  asym_map = _AsymMapper(imp_restraint.get_model(), components,
100  ignore_non_structural=True)
101  asyms = frozenset(
102  asym_map[p]
103  for p in IMP.get_input_particles(imp_restraint.get_inputs()))
104  asyms = sorted((a for a in asyms if a is not None),
105  key=operator.attrgetter('id'))
106  return asyms
107 
108 
109 class _EM3DRestraint(ihm.restraint.EM3DRestraint):
110  def __init__(self, imp_restraint, info, components, system):
111  # If a subunit contains any density, add the entire subunit to
112  # this restraint's assembly
113  asyms = _get_restraint_assembly(imp_restraint, components)
114 
115  assembly = ihm.Assembly(
116  asyms, name="EM subassembly",
117  description="All components that fit the EM map")
118 
119  self._filename = info['filename']
120  self._asyms = tuple(asyms)
121  p = IMP.mmcif.metadata._GMMParser()
122  r = p.parse_file(info['filename'])
123  super(_EM3DRestraint, self).__init__(
124  dataset=r['dataset'], assembly=assembly,
125  number_of_gaussians=r['number_of_gaussians'],
126  fitting_method='Gaussian mixture model')
127 
128  def _get_signature(self):
129  return ("EM3DRestraint", self._filename, self._asyms,
130  self.number_of_gaussians, self.fitting_method)
131 
132  def add_model_fit(self, imp_restraint, model):
133  info = _parse_restraint_info(imp_restraint.get_dynamic_info())
134  self.fits[model] = ihm.restraint.EM3DRestraintFit(
135  cross_correlation_coefficient=info['cross correlation'])
136 
137 
138 def _make_em2d_restraint(imp_restraint, info, components, system):
139  for i in range(len(info['image files'])):
140  yield _EM2DRestraint(imp_restraint, info, components, i)
141 
142 
143 class _EM2DRestraint(ihm.restraint.EM2DRestraint):
144  def __init__(self, imp_restraint, info, components, image_number):
145  asyms = _get_restraint_assembly(imp_restraint, components)
146 
147  assembly = ihm.Assembly(
148  asyms, name="2D EM subassembly",
149  description="All components that fit the EM images")
150 
151  self._image_number = image_number
152  self._filename = info['image files'][image_number]
153  self._asyms = tuple(asyms)
154 
155  loc = ihm.location.InputFileLocation(
156  self._filename,
157  details="Electron microscopy class average")
158  dataset = ihm.dataset.EM2DClassDataset(loc)
159 
160  super(_EM2DRestraint, self).__init__(
161  dataset=dataset, assembly=assembly,
162  segment=False,
163  number_raw_micrographs=info['micrographs number'] or None,
164  pixel_size_width=info['pixel size'],
165  pixel_size_height=info['pixel size'],
166  image_resolution=info['resolution'],
167  number_of_projections=info['projection number'])
168 
169  def _get_signature(self):
170  return ("EM2DRestraint", self._filename, self._asyms,
171  self.number_raw_micrographs, self.pixel_size_width,
172  self.pixel_size_height, self.image_resolution,
173  self.number_of_projections)
174 
175  def add_model_fit(self, imp_restraint, model):
176  info = _parse_restraint_info(imp_restraint.get_dynamic_info())
177  ccc = info['cross correlation'][self._image_number]
178  transform = self._get_transformation(model, info, self._image_number)
179  rot = transform.get_rotation()
180  rm = [[e for e in rot.get_rotation_matrix_row(i)] for i in range(3)]
181  self.fits[model] = ihm.restraint.EM2DRestraintFit(
182  cross_correlation_coefficient=ccc, rot_matrix=rm,
183  tr_vector=transform.get_translation())
184 
185  def _get_transformation(self, model, info, nimage):
186  """Get the transformation that places the model on image nimage"""
187  r = info['rotation'][nimage * 4: nimage * 4 + 4]
188  t = info['translation'][nimage * 3: nimage * 3 + 3]
191 
192 
193 class _SAXSRestraint(ihm.restraint.SASRestraint):
194  def __init__(self, imp_restraint, info, components, system):
195  asyms = _get_restraint_assembly(imp_restraint, components)
196 
197  assembly = ihm.Assembly(
198  asyms, name="SAXS subassembly",
199  description="All components that fit the SAXS profile")
200 
201  self._filename = info['filename']
202  self._asyms = tuple(asyms)
203  loc = ihm.location.InputFileLocation(
204  info['filename'], details='SAXS profile')
205  dataset = ihm.dataset.SASDataset(loc)
206  super(_SAXSRestraint, self).__init__(
207  dataset=dataset, assembly=assembly,
208  segment=False, fitting_method='IMP SAXS restraint',
209  fitting_atom_type=info['form factor type'],
210  multi_state=False)
211 
212  def _get_signature(self):
213  return ("SAXSRestraint", self._filename, self._asyms,
214  self.segment, self.fitting_method, self.fitting_atom_type,
215  self.multi_state)
216 
217  def add_model_fit(self, imp_restraint, model):
218  # We don't know the chi value; we only report a score
219  self.fits[model] = ihm.restraint.SASRestraintFit(chi_value=None)
220 
221 
222 class _CrossLinkRestraint(ihm.restraint.CrossLinkRestraint):
223  def __init__(self, imp_restraint, info, components, system):
224  self._info = info
225  loc = ihm.location.InputFileLocation(
226  info['filename'], details='Crosslinks')
227  dataset = ihm.dataset.CXMSDataset(loc)
228  linker = ihm.ChemDescriptor(
229  auth_name=info['linker author name'],
230  chemical_name=info.get('linker chemical name'),
231  smiles=info.get('linker smiles'),
232  smiles_canonical=info.get('linker smiles canonical'),
233  inchi=info.get('linker inchi'),
234  inchi_key=info.get('linker inchi key'))
235  super(_CrossLinkRestraint, self).__init__(
236  dataset=dataset, linker=linker)
237  # Map from IMP/RMF chain names to ihm.Entity
238  cmap = {e.description: e for e in system.entities}
239  dist = ihm.restraint.UpperBoundDistanceRestraint(info['linker length'])
240  asym = _AsymMapper(imp_restraint.get_model(), components)
241  self._add_all_links(IMP.RestraintSet.get_from(imp_restraint), cmap,
242  asym, dist)
243 
244  def _add_all_links(self, rset, cmap, asym, dist):
245  """Add info for each cross-link in the given RestraintSet"""
246  for link in rset.restraints:
247  # Recurse into any child RestraintSets
248  try:
249  child_rs = IMP.RestraintSet.get_from(link)
250  except ValueError:
251  child_rs = None
252  if child_rs:
253  self._add_all_links(child_rs, cmap, asym, dist)
254  else:
255  info = _parse_restraint_info(link.get_static_info())
256  # todo: handle ambiguous cross-links, fix residue numbering
257  r1 = cmap[info['protein1']].residue(info['residue1'])
258  r2 = cmap[info['protein2']].residue(info['residue2'])
259  ex_xl = ihm.restraint.ExperimentalCrossLink(residue1=r1,
260  residue2=r2)
261  self.experimental_cross_links.append([ex_xl])
262  # todo: handle multiple contributions
263  m = link.get_model()
264  endp1, endp2 = info['endpoints']
265  asym1 = asym[endp1]
266  asym2 = asym[endp2]
267  if _get_by_residue(m, endp1) and _get_by_residue(m, endp2):
268  cls = ihm.restraint.ResidueCrossLink
269  else:
270  cls = ihm.restraint.FeatureCrossLink
271  xl = cls(ex_xl, asym1=asym1, asym2=asym2, distance=dist,
272  restrain_all=False,
273  psi=_get_scale(m, info['psis'][0]),
274  sigma1=_get_scale(m, info['sigmas'][0]),
275  sigma2=_get_scale(m, info['sigmas'][1]))
276  self.cross_links.append(xl)
277 
278  def _get_signature(self):
279  # Assume that if we have the same restraint info (linker, csv file)
280  # and the same number of links, it is the same restraint.
281  # dict is not hashable, but a tuple of its items is.
282  return ("CXMSRestraint", tuple(self._info.items()),
283  len(self.cross_links), len(self.experimental_cross_links))
284 
285  def add_model_fit(self, imp_restraint, model):
286  pass # todo
287 
288 
289 class _GeometricRestraint(ihm.restraint.GeometricRestraint):
290  """Base for all geometric restraints"""
291 
292  def __init__(self, imp_restraint, info, components, system):
293  self._info = info
294  asym = _AsymMapper(imp_restraint.get_model(), components)
295  super(_GeometricRestraint, self).__init__(
296  dataset=None,
297  geometric_object=self._geom_object,
298  feature=asym.get_feature(
299  IMP.get_input_particles(imp_restraint.get_inputs())),
300  distance=self._get_distance(info),
301  harmonic_force_constant=1. / info['sigma'],
302  restrain_all=True)
303 
304  def _get_distance(self, info):
305  pass
306 
307  def add_model_fit(self, imp_restraint, model):
308  pass
309 
310  def _get_signature(self):
311  return ("GeometricRestraint", self.feature, self.geometric_object,
312  tuple(self._info.items()))
313 
314 
315 class _ZAxialRestraint(_GeometricRestraint):
316  """Handle an IMP.npc.ZAxialRestraint"""
317  _geom_object = ihm.geometry.XYPlane()
318 
319  def _get_distance(self, info):
320  return ihm.restraint.LowerUpperBoundDistanceRestraint(
321  info['lower bound'], info['upper bound'])
322 
323 
324 class _AllRestraints(object):
325  """Map IMP restraints to mmCIF objects"""
326  _typemap = {
327  "IMP.isd.GaussianEMRestraint": _EM3DRestraint,
328  "IMP.pmi.CrossLinkingMassSpectrometryRestraint": _CrossLinkRestraint,
329  "IMP.em2d.PCAFitRestraint": _make_em2d_restraint,
330  "IMP.npc.ZAxialPositionRestraint": _ZAxialRestraint,
331  "IMP.saxs.Restraint": _SAXSRestraint}
332 
333  def __init__(self, system, components):
334  self._system = system
335  self._components = components
336  self._seen_restraints = {}
337 
338  def add(self, restraint):
339  """Add and return a new restraint"""
340  sig = restraint._get_signature()
341  if sig not in self._seen_restraints:
342  self._system.restraints.append(restraint)
343  self._seen_restraints[sig] = restraint
344  return self._seen_restraints[sig]
345 
346  def handle(self, restraint, ihm_models):
347  """Handle an individual IMP restraint.
348  Yield wrapped version(s) of the restraint if it is handled in
349  mmCIF (one IMP restraint may map to multiple mmCIF restraints).
350  These may be new or existing restraint objects."""
351  info = _parse_restraint_info(restraint.get_static_info())
352  if 'type' in info and info['type'] in self._typemap:
353  cif_rs = self._typemap[info['type']](restraint,
354  info, self._components,
355  self._system)
356  try:
357  _ = iter(cif_rs)
358  except TypeError:
359  cif_rs = [cif_rs]
360  for cif_r in cif_rs:
361  cif_r = self.add(cif_r)
362  for model in ihm_models:
363  cif_r.add_model_fit(restraint, model)
364  yield cif_r
Simple 3D transformation class.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:158
A decorator to associate a particle with a part of a protein/DNA/RNA.
Definition: Fragment.h:20
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: atom/Atom.h:245
Add scale parameter to particle.
Definition: Scale.h:24
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Scale.h:29
3D rotation class.
Definition: Rotation3D.h:52
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
A decorator for a residue.
Definition: Residue.h:137
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Chain.h:82
static bool get_is_setup(Model *m, ParticleIndex p)
Check if the particle has the needed attributes for a cast to succeed.
VectorD< 3 > Vector3D
Definition: VectorD.h:425
Store info for a chain of a protein.
Definition: Chain.h:61
Classes to extract metadata for various input files.