IMP logo
IMP Reference Guide  2.19.0
The Integrative Modeling Platform
restraints/crosslinking.py
1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling cross-linking data.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 import IMP.algebra
9 import IMP.atom
10 import IMP.isd
11 import IMP.container
12 import IMP.pmi.tools
13 import IMP.pmi.output
15 import IMP.pmi.restraints
16 from math import log
17 from collections import defaultdict
18 import itertools
19 import operator
20 import warnings
21 
22 
23 class CrossLinkingMassSpectrometryRestraint(IMP.pmi.restraints.RestraintBase):
24  """Setup cross-link distance restraints from mass spectrometry data.
25  The noise in the data and the structural uncertainty of cross-linked
26  amino-acids is inferred using Bayes theory of probability
27  @note Wraps an IMP::isd::CrossLinkMSRestraint
28  """
29  _include_in_rmf = True
30 
31  def __init__(self, root_hier, database=None, length=10.0, resolution=None,
32  slope=0.02, label=None, filelabel="None",
33  attributes_for_label=None, linker=None, weight=1.):
34  """Constructor.
35  @param root_hier The canonical hierarchy containing all the states
36  @param database The IMP.pmi.io.crosslink.CrossLinkDataBase
37  object that contains the cross-link dataset
38  @param length maximal cross-linker length (including the residue
39  sidechains)
40  @param resolution what representation resolution should the cross-link
41  restraint be applied to.
42  @param slope The slope of a distance-linear scoring function that
43  funnels the score when the particles are
44  too far away. Suggested value 0.02.
45  @param label the extra text to label the restraint so that it is
46  searchable in the output
47  @param filelabel automatically generated file containing
48  missing/included/excluded
49  cross-links will be labeled using this text
50  @param attributes_for_label
51  @param weight Weight of restraint
52  @param linker description of the chemistry of the linker itself, as
53  an ihm.ChemDescriptor object
54  (see https://python-ihm.readthedocs.io/en/latest/main.html#ihm.ChemDescriptor).
55  Common cross-linkers can be found in the `ihm.cross_linkers`
56  module.
57  """ # noqa: E501
58 
59  model = root_hier.get_model()
60 
61  super(CrossLinkingMassSpectrometryRestraint, self).__init__(
62  model, weight=weight, label=label,
63  restraint_set_class=IMP.pmi.CrossLinkRestraintSet)
64 
65  if database is None:
66  raise Exception("You must pass a database")
67  if not isinstance(database, IMP.pmi.io.crosslink.CrossLinkDataBase):
68  raise TypeError(
69  "CrossLinkingMassSpectrometryRestraint: database should "
70  "be an IMP.pmi.io.crosslink.CrossLinkDataBase object")
71  self.database = database
72 
73  if resolution == 0 or resolution is None:
74  raise Exception("You must pass a resolution and it can't be zero")
75 
76  indb = open("included." + filelabel + ".xl.db", "w")
77  exdb = open("excluded." + filelabel + ".xl.db", "w") # noqa: F841
78  midb = open("missing." + filelabel + ".xl.db", "w")
79 
80  self.linker = linker
81  if linker is None:
82  raise ValueError(
83  "No linker chemistry specified. A linker must be given, as an "
84  "ihm.ChemDescriptor object (see the "
85  "CrossLinkingMassSpectrometryRestraint documentation).")
86  self.rs.set_name(self.rs.get_name() + "_Data")
87  self.rspsi = self._create_restraint_set("PriorPsi")
88  self.rssig = self._create_restraint_set("PriorSig")
89  self.rslin = self._create_restraint_set("Linear")
90  # Add custom metadata (will be saved in RMF output)
91  self.rs.set_metadata(self.database.name, length, slope)
92  if linker:
93  self.rs.set_linker_auth_name(linker.auth_name)
94  for attr in ('chemical_name', 'smiles', 'smiles_canonical',
95  'inchi', 'inchi_key'):
96  val = getattr(linker, attr)
97  if val:
98  getattr(self.rs, "set_linker_" + attr)(val)
99 
100  # dummy linear restraint used for Chimera display
101  self.linear = IMP.core.Linear(0, 0.0)
102  self.linear.set_slope(0.0)
103  dps2 = IMP.core.DistancePairScore(self.linear)
104 
105  self.psi_is_sampled = True
106  self.sigma_is_sampled = True
107  self.psi_dictionary = {}
108  self.sigma_dictionary = {}
109  self.xl_list = []
110  self.outputlevel = "low"
111 
112  restraints = []
113 
114  xl_groups = [p.get_cross_link_group(self)
115  for p, state in IMP.pmi.tools._all_protocol_outputs(
116  root_hier)]
117 
118  # first add all the molecule copies as clones to the database
119  copies_to_add = defaultdict(int)
120  print('gathering copies')
121  for xlid in self.database.xlid_iterator():
122  for xl in self.database[xlid]:
123  r1 = xl[self.database.residue1_key]
124  c1 = xl[self.database.protein1_key]
125  r2 = xl[self.database.residue2_key]
126  c2 = xl[self.database.protein2_key]
127  for c, r in ((c1, r1), (c2, r2)):
128  if c in copies_to_add:
129  continue
130  sel = IMP.atom.Selection(
131  root_hier, state_index=0, molecule=c, residue_index=r,
132  resolution=resolution).get_selected_particles()
133  if len(sel) > 0:
134  copies_to_add[c] = len(sel)-1
135  print(copies_to_add)
136  for molname in copies_to_add:
137  if copies_to_add[molname] == 0:
138  continue
140  self.database.protein1_key, operator.eq, molname)
141  self.database.set_value(self.database.protein1_key,
142  molname + '.0', fo1)
144  self.database.protein2_key, operator.eq, molname)
145  self.database.set_value(self.database.protein2_key,
146  molname + '.0', fo2)
147  for ncopy in range(copies_to_add[molname]):
148  self.database.clone_protein('%s.0' % molname,
149  '%s.%i' % (molname, ncopy + 1))
150  print('done pmi2 prelims')
151 
152  for xlid in self.database.xlid_iterator():
153  new_contribution = True
154  for xl in self.database[xlid]:
155 
156  r1 = xl[self.database.residue1_key]
157  c1 = xl[self.database.protein1_key]
158  r2 = xl[self.database.residue2_key]
159  c2 = xl[self.database.protein2_key]
160 
161  name1 = c1
162  name2 = c2
163  copy1 = 0
164  copy2 = 0
165  if '.' in c1:
166  name1, copy1 = c1.split('.')
167  if '.' in c2:
168  name2, copy2 = c2.split('.')
169 
170  # todo: check that offset is handled correctly
171  ex_xls = [
172  (p[0].add_experimental_cross_link(r1, name1, r2, name2,
173  group), group)
174  for p, group in zip(
175  IMP.pmi.tools._all_protocol_outputs(root_hier),
176  xl_groups)]
177 
178  iterlist = range(len(IMP.atom.get_by_type(
179  root_hier, IMP.atom.STATE_TYPE)))
180  for nstate, r in enumerate(iterlist):
181  # loop over every state
182  xl[self.database.state_key] = nstate
183  xl[self.database.data_set_name_key] = self.label
184 
185  ps1 = IMP.atom.Selection(
186  root_hier, state_index=nstate, molecule=name1,
187  copy_index=int(copy1), residue_index=r1,
188  resolution=resolution).get_selected_particles()
189  ps2 = IMP.atom.Selection(
190  root_hier, state_index=nstate, molecule=name2,
191  copy_index=int(copy2), residue_index=r2,
192  resolution=resolution).get_selected_particles()
193 
194  ps1 = [IMP.atom.Hierarchy(p) for p in ps1]
195  ps2 = [IMP.atom.Hierarchy(p) for p in ps2]
196 
197  if len(ps1) > 1:
198  raise ValueError(
199  "residue %d of chain %s selects multiple "
200  "particles %s" % (r1, c1, str(ps1)))
201  elif len(ps1) == 0:
202  warnings.warn(
203  "CrossLinkingMassSpectrometryRestraint: "
204  "residue %d of chain %s is not there" % (r1, c1),
206  midb.write(str(xl) + "\n")
207  continue
208 
209  if len(ps2) > 1:
210  raise ValueError(
211  "residue %d of chain %s selects multiple "
212  "particles %s" % (r2, c2, str(ps2)))
213  elif len(ps2) == 0:
214  warnings.warn(
215  "CrossLinkingMassSpectrometryRestraint: "
216  "residue %d of chain %s is not there" % (r2, c2),
218  midb.write(str(xl) + "\n")
219  continue
220 
221  p1 = ps1[0]
222  p2 = ps2[0]
223 
224  if p1 == p2 and r1 == r2:
225  warnings.warn(
226  "CrossLinkingMassSpectrometryRestraint: "
227  "same particle and same residue, skipping "
228  "cross-link", IMP.pmi.StructureWarning)
229  continue
230 
231  if new_contribution:
232  print("generating a new cross-link restraint")
233  new_contribution = False
235  self.model,
236  length,
237  slope)
238  dr.set_source_protein1(c1)
239  dr.set_source_protein2(c2)
240  dr.set_source_residue1(r1)
241  dr.set_source_residue2(r2)
242  restraints.append(dr)
243 
244  if self.database.sigma1_key not in xl.keys():
245  sigma1name = "SIGMA"
246  xl[self.database.sigma1_key] = sigma1name
247  else:
248  sigma1name = xl[self.database.sigma1_key]
249  sigma1 = self.create_sigma(sigma1name)
250 
251  if self.database.sigma2_key not in xl.keys():
252  sigma2name = "SIGMA"
253  xl[self.database.sigma2_key] = sigma2name
254  else:
255  sigma2name = xl[self.database.sigma2_key]
256  sigma2 = self.create_sigma(sigma2name)
257 
258  if self.database.psi_key not in xl.keys():
259  psiname = "PSI"
260  xl[self.database.psi_key] = psiname
261  else:
262  psiname = xl[self.database.psi_key]
263 
264  psi = self.create_psi(psiname)
265 
266  p1i = p1.get_particle_index()
267  xl["Particle1"] = p1
268  p2i = p2.get_particle_index()
269  xl["Particle2"] = p2
270  s1i = sigma1.get_particle().get_index()
271  xl["Particle_sigma1"] = sigma1
272  s2i = sigma2.get_particle().get_index()
273  xl["Particle_sigma2"] = sigma2
274  psii = psi.get_particle().get_index()
275  xl["Particle_psi"] = psi
276 
277  dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
278  xl["Restraint"] = dr
279 
280  print("--------------")
281  print("CrossLinkingMassSpectrometryRestraint: generating "
282  "cross-link restraint between")
283  print("CrossLinkingMassSpectrometryRestraint: residue %d "
284  "of chain %s and residue %d of chain %s"
285  % (r1, c1, r2, c2))
286  print("CrossLinkingMassSpectrometryRestraint: with sigma1 "
287  "%s sigma2 %s psi %s"
288  % (sigma1name, sigma2name, psiname))
289  print("CrossLinkingMassSpectrometryRestraint: between "
290  "particles %s and %s"
291  % (p1.get_name(), p2.get_name()))
292  print("==========================================\n")
293  for p, ex_xl in zip(IMP.pmi.tools._all_protocol_outputs(
294  root_hier),
295  ex_xls):
296  p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
297  sigma1, sigma2, psi, ex_xl[1])
298 
299  # check if the two residues belong to the same rigid body
302  IMP.core.RigidMember(p1).get_rigid_body() ==
303  IMP.core.RigidMember(p2).get_rigid_body()):
304  xl["IntraRigidBody"] = True
305  else:
306  xl["IntraRigidBody"] = False
307 
308  xl_label = self.database.get_short_cross_link_string(xl)
309  xl["ShortLabel"] = xl_label
310  dr.set_name(xl_label)
311 
312  if p1i != p2i:
313  pr = IMP.core.PairRestraint(self.model, dps2,
314  (p1i, p2i))
315  pr.set_name(xl_label)
316  self.rslin.add_restraint(pr)
317 
318  self.xl_list.append(xl)
319 
320  indb.write(str(xl) + "\n")
321 
322  if len(self.xl_list) == 0:
323  raise SystemError(
324  "CrossLinkingMassSpectrometryRestraint: no cross-link "
325  "was constructed")
326  self.xl_restraints = restraints
327  lw = IMP.isd.LogWrapper(restraints, 1.0)
328  self.rs.add_restraint(lw)
329  indb.close()
330  exdb.close()
331  midb.close()
332 
333  def __set_dataset(self, ds):
334  self.database.dataset = ds
335  dataset = property(lambda self: self.database.dataset, __set_dataset)
336 
337  def get_hierarchies(self):
338  """ get the hierarchy """
339  return self.prot
340 
342  """ get the restraint set """
343  return self.rs
344 
345  def get_restraints(self):
346  """ get the restraints in a list """
347  return self.xl_restraints
348 
350  """ get the dummy restraints to be displayed in the rmf file """
351  return self.rs, self.rssig, self.rspsi
352 
354  """ Get a list of tuples containing the particle pairs """
355  ppairs = []
356  for i in range(len(self.pairs)):
357  p0 = self.pairs[i][0]
358  p1 = self.pairs[i][1]
359  ppairs.append((p0, p1))
360  return ppairs
361 
362  def set_output_level(self, level="low"):
363  """ Set the output level of the output """
364  self.outputlevel = level
365 
366  def set_psi_is_sampled(self, is_sampled=True):
367  """ Switch on/off the sampling of psi particles """
368  self.psi_is_sampled = is_sampled
369 
370  def set_sigma_is_sampled(self, is_sampled=True):
371  """ Switch on/off the sampling of sigma particles """
372  self.sigma_is_sampled = is_sampled
373 
374  def create_sigma(self, name):
375  """ This is called internally. Creates a nuisance
376  on the structural uncertainty """
377  if name in self.sigma_dictionary:
378  return self.sigma_dictionary[name][0]
379 
380  sigmainit = 2.0
381  sigmaminnuis = 0.0000001
382  sigmamaxnuis = 1000.0
383  sigmamin = 0.01
384  sigmamax = 100.0
385  sigmatrans = 0.5
386  sigma = IMP.pmi.tools.SetupNuisance(
387  self.model, sigmainit, sigmaminnuis, sigmamaxnuis,
388  self.sigma_is_sampled, name=name).get_particle()
389  self.sigma_dictionary[name] = (
390  sigma,
391  sigmatrans,
392  self.sigma_is_sampled)
393  self.rssig.add_restraint(
395  self.model,
396  sigma,
397  1000000000.0,
398  sigmamax,
399  sigmamin))
400  return sigma
401 
402  def create_psi(self, name):
403  """ This is called internally. Creates a nuisance
404  on the data uncertainty """
405  if name in self.psi_dictionary:
406  return self.psi_dictionary[name][0]
407 
408  psiinit = 0.25
409  psiminnuis = 0.0000001
410  psimaxnuis = 0.4999999
411  psimin = 0.01
412  psimax = 0.49
413  psitrans = 0.1
414  psi = IMP.pmi.tools.SetupNuisance(
415  self.model, psiinit, psiminnuis, psimaxnuis,
416  self.psi_is_sampled, name=name).get_particle()
417  self.psi_dictionary[name] = (psi, psitrans, self.psi_is_sampled)
418 
419  self.rspsi.add_restraint(
420  IMP.isd.UniformPrior(self.model, psi, 1000000000.0,
421  psimax, psimin))
422 
423  self.rspsi.add_restraint(IMP.isd.JeffreysRestraint(self.model, psi))
424  return psi
425 
426  def get_output(self):
427  """Get the output of the restraint to be used by the IMP.pmi.output
428  object"""
429  output = super(CrossLinkingMassSpectrometryRestraint,
430  self).get_output()
431 
432  for xl in self.xl_list:
433 
434  xl_label = xl["ShortLabel"]
435  ln = xl["Restraint"]
436  p0 = xl["Particle1"]
437  p1 = xl["Particle2"]
438  output["CrossLinkingMassSpectrometryRestraint_Score_" +
439  xl_label] = str(-log(ln.unprotected_evaluate(None)))
440 
441  d0 = IMP.core.XYZ(p0)
442  d1 = IMP.core.XYZ(p1)
443  output["CrossLinkingMassSpectrometryRestraint_Distance_" +
444  xl_label] = str(IMP.core.get_distance(d0, d1))
445 
446  for psiname in self.psi_dictionary:
447  output["CrossLinkingMassSpectrometryRestraint_Psi_" +
448  str(psiname) + self._label_suffix] = str(
449  self.psi_dictionary[psiname][0].get_scale())
450 
451  for sigmaname in self.sigma_dictionary:
452  output["CrossLinkingMassSpectrometryRestraint_Sigma_" +
453  str(sigmaname) + self._label_suffix] = str(
454  self.sigma_dictionary[sigmaname][0].get_scale())
455 
456  return output
457 
458  def get_movers(self):
459  """ Get all need data to construct a mover in IMP.pmi.dof class"""
460  movers = []
461  if self.sigma_is_sampled:
462  for sigmaname in self.sigma_dictionary:
463  mover_name = \
464  "Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" \
465  + str(sigmaname) + "_" + self.label
466  particle = self.sigma_dictionary[sigmaname][0]
467  maxstep = (self.sigma_dictionary[sigmaname][1])
469  [particle], IMP.FloatKeys([IMP.FloatKey("nuisance")]),
470  maxstep)
471  mv.set_name(mover_name)
472  movers.append(mv)
473 
474  if self.psi_is_sampled:
475  for psiname in self.psi_dictionary:
476  mover_name = \
477  "Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" + \
478  str(psiname) + "_" + self.label
479  particle = self.psi_dictionary[psiname][0]
480  maxstep = (self.psi_dictionary[psiname][1])
482  [particle], IMP.FloatKeys([IMP.FloatKey("nuisance")]),
483  maxstep)
484  mv.set_name(mover_name)
485  movers.append(mv)
486 
487  return movers
488 
490  """ Get the particles to be sampled by the IMP.pmi.sampler object """
491  ps = {}
492  if self.sigma_is_sampled:
493  for sigmaname in self.sigma_dictionary:
494  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" +
495  str(sigmaname) + self._label_suffix] =\
496  ([self.sigma_dictionary[sigmaname][0]],
497  self.sigma_dictionary[sigmaname][1])
498 
499  if self.psi_is_sampled:
500  for psiname in self.psi_dictionary:
501  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
502  str(psiname) + self._label_suffix] = \
503  ([self.psi_dictionary[psiname][0]],
504  self.psi_dictionary[psiname][1])
505 
506  return ps
507 
508 
509 class AtomicCrossLinkMSRestraint(IMP.pmi.restraints.RestraintBase):
510  """Setup cross-link distance restraints at atomic level
511  The "atomic" aspect is that it models the particle uncertainty with
512  a Gaussian. The noise in the data and the structural uncertainty of
513  cross-linked amino-acids is inferred using Bayes' theory of probability
514  @note Wraps an IMP::isd::AtomicCrossLinkMSRestraint
515  """
516 
517  _include_in_rmf = True
518 
519  def __init__(self, root_hier, xldb, atom_type="NZ", length=10.0,
520  slope=0.01, nstates=None, label=None,
521  nuisances_are_optimized=True, sigma_init=5.0,
522  psi_init=0.01, one_psi=True, filelabel=None, weight=1.):
523  """Constructor.
524  Automatically creates one "sigma" per cross-linked residue and one
525  "psis" per pair.
526  Other nuisance options are available.
527  @note Will return an error if the data+extra_sel don't specify
528  two particles per XL pair.
529  @param root_hier The root hierarchy on which you'll do selection
530  @param xldb CrossLinkDataBase object
531  @param atom_type Can either be "NZ" or "CA"
532  @param length The XL linker length
533  @param slope Linear term to add to the restraint function to
534  help when far away
535  @param nstates The number of states to model. Defaults to the
536  number of states in root.
537  @param label The output label for the restraint
538  @param nuisances_are_optimized Whether to optimize nuisances
539  @param sigma_init The initial value for all the sigmas
540  @param psi_init The initial value for all the psis
541  @param one_psi Use a single psi for all restraints (if False,
542  creates one per XL)
543  @param filelabel automatically generated file containing
544  missing/included/excluded
545  cross-links will be labeled using this text
546  @param weight Weight of restraint
547 
548  """
549 
550  # basic params
551  self.root = root_hier
552  rname = "AtomicXLRestraint"
553  super(AtomicCrossLinkMSRestraint, self).__init__(
554  self.root.get_model(), name="AtomicXLRestraint", label=label,
555  weight=weight)
556  self.xldb = xldb
557  self.length = length
558  self.sigma_is_sampled = nuisances_are_optimized
559  self.psi_is_sampled = nuisances_are_optimized
560 
561  if atom_type not in ("NZ", "CA"):
562  raise Exception("AtomicXLRestraint: atom_type must be NZ or CA")
563  self.atom_type = atom_type
564  self.nstates = nstates
565  if nstates is None:
566  self.nstates = len(IMP.atom.get_by_type(self.root,
567  IMP.atom.STATE_TYPE))
568  elif nstates != len(IMP.atom.get_by_type(self.root,
569  IMP.atom.STATE_TYPE)):
570  print("Warning: nstates is not the same as the number of states "
571  "in root")
572 
573  self.rs_psi = self._create_restraint_set("psi")
574  self.rs_sig = self._create_restraint_set("sigma")
575  self.rs_lin = self._create_restraint_set("linear")
576 
577  self.psi_dictionary = {}
578  self.sigma_dictionary = {}
579 
580  self.particles = defaultdict(set)
581  self.one_psi = one_psi
582  self.bonded_pairs = []
583  if self.one_psi:
584  print('creating a single psi for all XLs')
585  else:
586  print('creating one psi for each XL')
587 
588  # output logging file
589  if filelabel is not None:
590  indb = open("included." + filelabel + ".xl.db", "w") # noqa: F841
591  exdb = open("excluded." + filelabel + ".xl.db", "w") # noqa: F841
592  midb = open("missing." + filelabel + ".xl.db", "w")
593 
594  # Setup nuisances
595  self._create_sigma('sigma', sigma_init)
596  if one_psi:
597  self._create_psi('psi', psi_init)
598  else:
599  for xlid in self.xldb.xlid_iterator():
600  self._create_psi(xlid, psi_init)
601 
602  # create all the XLs
603  xlrs = []
604  for xlid in self.xldb.xlid_iterator():
605  # create restraint for this data point
606  if one_psi:
607  psip = self.psi_dictionary['psi'][0].get_particle_index()
608  else:
609  psip = self.psi_dictionary[
610  unique_id][0].get_particle_index() # noqa: F821
611  r = IMP.isd.AtomicCrossLinkMSRestraint(self.model,
612  self.length,
613  psip,
614  slope,
615  True)
616  num_contributions = 0
617 
618  # add a contribution for each XL ambiguity option within each state
619  for nstate in self.nstates:
620  for xl in self.xldb[xlid]:
621  r1 = xl[self.xldb.residue1_key]
622  c1 = xl[self.xldb.protein1_key].strip()
623  r2 = xl[self.xldb.residue2_key]
624  c2 = xl[self.xldb.protein2_key].strip()
625 
626  # perform selection. these may contain multiples if
627  # Copies are used
628  ps1 = IMP.atom.Selection(
629  self.root, state_index=nstate,
630  atom_type=IMP.atom.AtomType(self.atom_type),
631  molecule=c1,
632  residue_index=r1).get_selected_particles()
633  ps2 = IMP.atom.Selection(
634  self.root, state_index=nstate,
635  atom_type=IMP.atom.AtomType(self.atom_type),
636  molecule=c2,
637  residue_index=r2).get_selected_particles()
638  if len(ps1) == 0:
639  warnings.warn("AtomicXLRestraint: residue %d of "
640  "chain %s is not there" % (r1, c1),
642  if filelabel is not None:
643  midb.write(str(xl) + "\n")
644  continue
645 
646  if len(ps2) == 0:
647  warnings.warn("AtomicXLRestraint: residue %d of "
648  "chain %s is not there" % (r2, c2),
650  if filelabel is not None:
651  midb.write(str(xl) + "\n")
652  continue
653 
654  # figure out sig1 and sig2 based on num XLs
655  '''
656  num1=num_xls_per_res[str(xl.r1)]
657  num2=num_xls_per_res[str(xl.r2)]
658  if num1<sig_threshold:
659  sig1=self.sig_low
660  else:
661  sig1=self.sig_high
662  if num2<sig_threshold:
663  sig2=self.sig_low
664  else:
665  sig2=self.sig_high
666  '''
667  sig1 = self.sigma_dictionary['sigma'][0]
668  sig2 = self.sigma_dictionary['sigma'][0]
669 
670  # add each copy contribution to restraint
671  for p1, p2 in itertools.product(ps1, ps2):
672  if p1 == p2:
673  continue
674  self.particles[nstate] |= set((p1, p2))
675  r.add_contribution(
676  [p1.get_index(), p2.get_index()],
677  [sig1.get_particle_index(),
678  sig2.get_particle_index()])
679  num_contributions += 1
680 
681  if num_contributions > 0:
682  print('XL:', xlid, 'num contributions:', num_contributions)
683  xlrs.append(r)
684  if len(xlrs) == 0:
685  raise Exception("You didn't create any XL restraints")
686  print('created', len(xlrs), 'XL restraints')
687  rname = self.rs.get_name()
688  self.rs = IMP.isd.LogWrapper(xlrs, self.weight)
689  self.rs.set_name(rname)
690  self.rs.set_weight(self.weight)
691  self.restraint_sets = [self.rs] + self.restraint_sets[1:]
692 
693  def get_hierarchy(self):
694  return self.prot
695 
696  def _create_sigma(self, name, sigmainit):
697  """ This is called internally. Creates a nuisance
698  on the structural uncertainty """
699  if name in self.sigma_dictionary:
700  return self.sigma_dictionary[name][0]
701 
702  sigmaminnuis = 0.0000001
703  sigmamaxnuis = 1000.0
704  sigmamin = 0.01
705  sigmamax = 100.0
706  sigmatrans = 0.5
707  sigma = IMP.pmi.tools.SetupNuisance(
708  self.model, sigmainit, sigmaminnuis, sigmamaxnuis,
709  self.sigma_is_sampled).get_particle()
710  self.sigma_dictionary[name] = (
711  sigma, sigmatrans, self.sigma_is_sampled)
712  self.rs_sig.add_restraint(IMP.isd.UniformPrior(
713  self.model, sigma, 1000000000.0, sigmamax, sigmamin))
714  return sigma
715 
716  def _create_psi(self, name, psiinit):
717  """ This is called internally. Creates a nuisance
718  on the data uncertainty """
719  if name in self.psi_dictionary:
720  return self.psi_dictionary[name][0]
721 
722  psiminnuis = 0.0000001
723  psimaxnuis = 0.4999999
724  psimin = 0.01
725  psimax = 0.49
726  psitrans = 0.1
727  psi = IMP.pmi.tools.SetupNuisance(self.model,
728  psiinit,
729  psiminnuis,
730  psimaxnuis,
731  self.psi_is_sampled).get_particle()
732  self.psi_dictionary[name] = (
733  psi,
734  psitrans,
735  self.psi_is_sampled)
736 
737  self.rs_psi.add_restraint(
739  self.model,
740  psi,
741  1000000000.0,
742  psimax,
743  psimin))
744 
745  self.rs_psi.add_restraint(IMP.isd.JeffreysRestraint(self.model, psi))
746  return psi
747 
749  """ create dummy harmonic restraints for each XL but don't add to model
750  Makes it easy to see each contribution to each XL in RMF
751  """
752  class MyGetRestraint(object):
753  def __init__(self, rs):
754  self.rs = rs
755 
756  def get_restraint_for_rmf(self):
757  return self.rs
758 
759  dummy_model = IMP.Model()
760  hps = IMP.core.HarmonicDistancePairScore(self.length, 1.0)
761  dummy_rs = []
762  for nxl in range(self.get_number_of_restraints()):
763  xl = IMP.isd.AtomicCrossLinkMSRestraint.get_from(
764  self.rs.get_restraint(nxl))
765  rs = IMP.RestraintSet(dummy_model, 'atomic_xl_'+str(nxl))
766  for ncontr in range(xl.get_number_of_contributions()):
767  ps = xl.get_contribution(ncontr)
769  hps, [self.model.get_particle(p) for p in ps],
770  'xl%i_contr%i' % (nxl, ncontr))
771  rs.add_restraint(dr)
772  dummy_rs.append(MyGetRestraint(rs))
773  return dummy_rs
774 
776  """ Get the particles to be sampled by the IMP.pmi.sampler object """
777  ps = {}
778  if self.sigma_is_sampled:
779  for sigmaname in self.sigma_dictionary:
780  ps["Nuisances_AtomicCrossLinkingMSRestraint_Sigma_" +
781  str(sigmaname) + self._label_suffix] = \
782  ([self.sigma_dictionary[sigmaname][0]],
783  self.sigma_dictionary[sigmaname][1])
784  if self.psi_is_sampled:
785  for psiname in self.psi_dictionary:
786  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
787  str(psiname) + self._label_suffix] = \
788  ([self.psi_dictionary[psiname][0]],
789  self.psi_dictionary[psiname][1])
790  return ps
791 
792  def get_bonded_pairs(self):
793  return self.bonded_pairs
794 
795  def get_number_of_restraints(self):
796  return self.rs.get_number_of_restraints()
797 
798  def __repr__(self):
799  return 'XL restraint with ' + \
800  str(len(self.rs.get_restraint(0).get_number_of_restraints())) \
801  + ' data points'
802 
803  def load_nuisances_from_stat_file(self, in_fn, nframe):
804  """Read a stat file and load all the sigmas.
805  This is potentially quite stupid.
806  It's also a hack since the sigmas should be stored in the RMF file.
807  Also, requires one sigma and one psi for ALL XLs.
808  """
809  import subprocess
810  sig_val = float(subprocess.check_output(
811  ["process_output.py", "-f", in_fn, "-s",
812  "AtomicXLRestraint_sigma"]).split('\n>')[1+nframe])
813  psi_val = float(subprocess.check_output(
814  ["process_output.py", "-f", in_fn, "-s",
815  "AtomicXLRestraint_psi"]).split('\n>')[1+nframe])
816  for nxl in range(self.rs.get_number_of_restraints()):
817  xl = IMP.isd.AtomicCrossLinkMSRestraint.get_from(
818  self.rs.get_restraint(nxl))
819  psip = xl.get_psi()
820  IMP.isd.Scale(self.model, psip).set_scale(psi_val)
821  for contr in range(xl.get_number_of_contributions()):
822  sig1, sig2 = xl.get_contribution_sigmas(contr)
823  IMP.isd.Scale(self.model, sig1).set_scale(sig_val)
824 
825  print('loaded nuisances from file')
826 
827  def plot_violations(self, out_prefix, max_prob_for_violation=0.1,
828  min_dist_for_violation=1e9, coarsen=False,
829  limit_to_chains=None, exclude_chains=''):
830  """Create CMM files, one for each state, of all cross-links.
831  will draw in GREEN if non-violated in all states (or if only one state)
832  will draw in PURPLE if non-violated only in a subset of states
833  (draws nothing elsewhere)
834  will draw in RED in ALL states if all violated
835  (if only one state, you'll only see green and red)
836 
837  @param out_prefix Output xlink files prefix
838  @param max_prob_for_violation It's a violation if the probability
839  is below this
840  @param min_dist_for_violation It's a violation if the min dist
841  is above this
842  @param coarsen Use CA positions
843  @param limit_to_chains Try to visualize just these chains
844  @param exclude_to_chains Try to NOT visualize these chains
845  """
846  print('going to calculate violations and plot CMM files')
847  all_stats = self.get_best_stats()
848  all_dists = [s["low_dist"] for s in all_stats]
849 
850  # prepare one output file per state
851  out_fns = []
852  out_nvs = []
853  state_info = []
854  cmds = defaultdict(set)
855  for nstate in range(self.nstates):
856  outf = open(out_prefix+str(nstate)+'.cmm', 'w')
857  outf.write('<marker_set name="xlinks_state%i"> \n' % nstate)
858  out_fns.append(outf)
859  out_nvs.append(0)
860  print('will limit to', limit_to_chains)
861  print('will exclude', exclude_chains)
862  state_info.append(self.get_best_stats(nstate,
863  limit_to_chains,
864  exclude_chains))
865 
866  for nxl in range(self.rs.get_number_of_restraints()):
867  # for this XL, check which states passed
868  npass = []
869  nviol = []
870  for nstate in range(self.nstates):
871  prob = state_info[nstate][nxl]["prob"]
872  low_dist = state_info[nstate][nxl]["low_dist"]
873  if (prob < max_prob_for_violation or
874  low_dist > min_dist_for_violation):
875  nviol.append(nstate)
876  else:
877  npass.append(nstate)
878 
879  # special coloring when all pass or all fail
880  all_pass = False
881  all_viol = False
882  if len(npass) == self.nstates:
883  all_pass = True
884  elif len(nviol) == self.nstates:
885  all_viol = True
886 
887  # finally, color based on above info
888  print(nxl, 'state dists:',
889  [state_info[nstate][nxl]["low_dist"]
890  for nstate in range(self.nstates)],
891  'viol states:', nviol, 'all viol?', all_viol)
892  for nstate in range(self.nstates):
893  if all_pass:
894  r = 0.365
895  g = 0.933
896  b = 0.365
897  continue
898  elif all_viol:
899  r = 0.980
900  g = 0.302
901  b = 0.247
902  continue
903  else:
904  if nstate in nviol:
905  continue
906  else:
907  r = 0.365
908  g = 0.933
909  b = 0.365
910  # now only showing if UNIQUELY PASSING in this state
911  pp = state_info[nstate][nxl]["low_pp"]
912  c1 = IMP.core.XYZ(self.model, pp[0]).get_coordinates()
913  c2 = IMP.core.XYZ(self.model, pp[1]).get_coordinates()
914 
916  IMP.atom.Atom(self.model, pp[0])).get_index()
917  ch1 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model, pp[0]))
919  IMP.atom.Atom(self.model, pp[0])).get_index()
920  ch2 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model, pp[0]))
921 
922  cmds[nstate].add((ch1, r1))
923  cmds[nstate].add((ch2, r2))
924 
925  outf = out_fns[nstate]
926  nv = out_nvs[nstate]
927  outf.write('<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" '
928  'radius="0.8" r="%.2f" g="%.2f" b="%.2f"/> \n'
929  % (nv, c1[0], c1[1], c1[2], r, g, b))
930  outf.write('<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" '
931  'radius="0.8" r="%.2f" g="%.2f" b="%.2f"/> \n'
932  % (nv+1, c2[0], c2[1], c2[2], r, g, b))
933  outf.write('<link id1= "%d" id2="%d" radius="0.8" '
934  'r="%.2f" g="%.2f" b="%.2f"/> \n'
935  % (nv, nv+1, r, g, b))
936  out_nvs[nstate] += 2
937 
938  for nstate in range(self.nstates):
939  out_fns[nstate].write('</marker_set>\n')
940  out_fns[nstate].close()
941  cmd = ''
942  for ch, r in cmds[nstate]:
943  cmd += '#%i:%i.%s ' % (nstate, r, ch)
944  print(cmd)
945 
946  return all_dists
947 
948  def _get_contribution_info(self, xl, ncontr, use_CA=False):
949  """Return the particles at that contribution. If requested will
950  return CA's instead"""
951  idx1 = xl.get_contribution(ncontr)[0]
952  idx2 = xl.get_contribution(ncontr)[1]
953  if use_CA:
954  sel = IMP.atom.Selection(
955  IMP.atom.get_residue(IMP.atom.Atom(self.model, idx1)),
956  atom_type=IMP.atom.AtomType("CA"))
957  idx1 = sel.get_selected_particle_indexes()[0]
958  sel = IMP.atom.Selection(
959  IMP.atom.get_residue(IMP.atom.Atom(self.model, idx2)),
960  atom_type=IMP.atom.AtomType("CA"))
961  idx2 = sel.get_selected_particle_indexes()[0]
963  IMP.core.XYZ(self.model, idx1).get_coordinates(),
964  IMP.core.XYZ(self.model, idx2).get_coordinates())
965  return idx1, idx2, dist
966 
967  def get_best_stats(self, limit_to_state=None, limit_to_chains=None,
968  exclude_chains='', use_CA=False):
969  ''' return the probability, best distance, two coords, and possibly
970  the psi for each xl
971  @param limit_to_state Only examine contributions from one state
972  @param limit_to_chains Returns the particles for certain "easy
973  to visualize" chains
974  @param exclude_chains Even if you limit, don't let one end be in
975  this list. Only works if you also limit chains
976  @param use_CA Limit to CA particles
977  '''
978  ret = []
979  for nxl in range(self.rs.get_number_of_restraints()):
980  this_info = {}
981  xl = IMP.isd.AtomicCrossLinkMSRestraint.get_from(
982  self.rs.get_restraint(nxl))
983  low_dist = 1e6
984  low_pp = None
985  state_contrs = []
986  low_pp_lim = None
987  low_dist_lim = 1e6
988  for contr in range(xl.get_number_of_contributions()):
989  pp = xl.get_contribution(contr)
990  if use_CA:
991  sel = IMP.atom.Selection(
992  IMP.atom.get_residue(IMP.atom.Atom(self.model, pp[0])),
993  atom_type=IMP.atom.AtomType("CA"))
994  idx1 = sel.get_selected_particle_indexes()[0]
995  sel = IMP.atom.Selection(
996  IMP.atom.get_residue(IMP.atom.Atom(self.model, pp[1])),
997  atom_type=IMP.atom.AtomType("CA"))
998  idx2 = sel.get_selected_particle_indexes()[0]
999  pp = [idx1, idx2]
1000  if limit_to_state is not None:
1001  nstate = IMP.atom.get_state_index(
1002  IMP.atom.Atom(self.model, pp[0]))
1003  if nstate != limit_to_state:
1004  continue
1005  state_contrs.append(contr)
1006  dist = IMP.core.get_distance(IMP.core.XYZ(self.model, pp[0]),
1007  IMP.core.XYZ(self.model, pp[1]))
1008  if limit_to_chains is not None:
1009  c1 = IMP.atom.get_chain_id(
1010  IMP.atom.Atom(self.model, pp[0]))
1011  c2 = IMP.atom.get_chain_id(
1012  IMP.atom.Atom(self.model, pp[1]))
1013  if (c1 in limit_to_chains or c2 in limit_to_chains) and (
1014  c1 not in exclude_chains
1015  and c2 not in exclude_chains):
1016  if dist < low_dist_lim:
1017  low_dist_lim = dist
1018  low_pp_lim = pp
1019  if dist < low_dist:
1020  low_dist = dist
1021  low_pp = pp
1022  if limit_to_state is not None:
1023  this_info["prob"] = xl.evaluate_for_contributions(
1024  state_contrs, None)
1025  else:
1026  this_info["prob"] = xl.unprotected_evaluate(None)
1027  if limit_to_chains is not None:
1028  this_info["low_pp"] = low_pp_lim
1029  else:
1030  this_info["low_pp"] = low_pp
1031 
1032  this_info["low_dist"] = low_dist
1033  if not self.one_psi:
1034  pval = IMP.isd.Scale(self.model, xl.get_psi()).get_scale()
1035  this_info["psi"] = pval
1036  ret.append(this_info)
1037  return ret
1038 
1039  def print_stats(self):
1040  stats = self.get_best_stats()
1041  for nxl, s in enumerate(stats):
1042  print(s["low_dist"])
1043 
1044  def get_output(self):
1045  output = super(AtomicCrossLinkMSRestraint, self).get_output()
1046 
1047  # HACK to make it easier to see the few sigmas
1048  # output["AtomicXLRestraint_sigma"] = self.sigma.get_scale()
1049  # if self.one_psi:
1050  # output["AtomicXLRestraint_psi"] = self.psi.get_scale()
1051  #
1052 
1053  # count distances above length
1054  bad_count = 0
1055  stats = self.get_best_stats()
1056  for nxl, s in enumerate(stats):
1057  if s['low_dist'] > 20.0:
1058  bad_count += 1
1059  output["AtomicXLRestraint_%i_%s" % (nxl, "Prob")] = str(s['prob'])
1060  output["AtomicXLRestraint_%i_%s"
1061  % (nxl, "BestDist")] = str(s['low_dist'])
1062  if not self.one_psi:
1063  output["AtomicXLRestraint_%i_%s" % (nxl, "psi")] \
1064  = str(s['psi'])
1065  output["AtomicXLRestraint_NumViol"] = str(bad_count)
1066  return output
1067 
1068 
1069 class CysteineCrossLinkRestraint(object):
1070  def __init__(self, root_hier, filename, cbeta=False,
1071  betatuple=(0.03, 0.1),
1072  disttuple=(0.0, 25.0, 1000),
1073  omegatuple=(1.0, 1000.0, 50),
1074  sigmatuple=(0.3, 0.3, 1),
1075  betaissampled=True,
1076  sigmaissampled=False,
1077  weightissampled=True,
1078  epsilonissampled=True):
1079  # the file must have residue1 chain1 residue2 chain2 fractionvalue
1080  # epsilonname
1081  # epsilonname is a name for the epsilon particle that must be used
1082  # for that particular residue pair, eg, "Epsilon-Intra-Solvent", or
1083  # "Epsilon-Solvent-Membrane", etc.
1084 
1085  self.m = root_hier.get_model()
1086  self.rs = IMP.RestraintSet(self.m, 'Cysteine_Crosslink')
1087  self.cbeta = cbeta
1088  self.epsilonmaxtrans = 0.01
1089  self.sigmamaxtrans = 0.1
1090  self.betamaxtrans = 0.01
1091  self.weightmaxtrans = 0.1
1092  self.label = "None"
1093  self.outputlevel = "low"
1094  self.betaissampled = betaissampled
1095  self.sigmaissampled = sigmaissampled
1096  self.weightissampled = weightissampled
1097  self.epsilonissampled = epsilonissampled
1098 
1099  betalower = betatuple[0]
1100  betaupper = betatuple[1]
1101  beta = 0.04
1102  sigma = 10.0
1103  betangrid = 30
1104 
1105  # beta
1106  self.beta = IMP.pmi.tools.SetupNuisance(
1107  self.m,
1108  beta,
1109  betalower,
1110  betaupper,
1111  betaissampled).get_particle(
1112  )
1113  # sigma
1114  self.sigma = IMP.pmi.tools.SetupNuisance(
1115  self.m,
1116  sigma,
1117  sigmatuple[0],
1118  sigmatuple[1],
1119  sigmaissampled).get_particle()
1120  # population particle
1121  self.weight = IMP.pmi.tools.SetupWeight(
1122  self.m,
1123  isoptimized=False).get_particle(
1124  )
1125 
1126  # read the file
1127  fl = IMP.pmi.tools.open_file_or_inline_text(filename)
1128  # determine the upperlimit for epsilon
1129  # and also how many epsilon are needed
1130  self.epsilons = {}
1131  data = []
1132  for line in fl:
1133  t = line.split()
1134  if t[0][0] == "#":
1135  continue
1136  if t[5] in self.epsilons:
1137  if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
1138  self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
1139  else:
1140  self.epsilons[t[5]] = IMP.pmi.tools.SetupNuisance(
1141  self.m, 0.01, 0.01, 1.0 - float(t[4]),
1142  epsilonissampled).get_particle()
1143  up = self.epsilons[t[5]].get_upper()
1144  low = self.epsilons[t[5]].get_lower()
1145  if up < low:
1146  self.epsilons[t[5]].set_upper(low)
1147 
1148  data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
1149 
1150  # create CrossLinkData
1151  if not self.cbeta:
1152  crossdata = IMP.pmi.tools.get_cross_link_data(
1153  "cysteine", "cysteine_CA_FES.txt.standard",
1154  disttuple, omegatuple, sigmatuple, disttuple[1],
1155  disttuple[1], 1)
1156  else:
1157  crossdata = IMP.pmi.tools.get_cross_link_data(
1158  "cysteine", "cysteine_CB_FES.txt.standard",
1159  disttuple, omegatuple, sigmatuple, disttuple[1],
1160  disttuple[1], 1)
1161 
1162  # create grids needed by CysteineCrossLinkData
1163  fmod_grid = IMP.pmi.tools.get_grid(0.0, 1.0, 300, True)
1164  omega2_grid = IMP.pmi.tools.get_log_grid(0.001, 10000.0, 100)
1165  beta_grid = IMP.pmi.tools.get_log_grid(betalower, betaupper, betangrid)
1166 
1167  for d in data:
1168  print("--------------")
1169  print("CysteineCrossLink: attempting to create a restraint "
1170  + str(d))
1171  resid1 = d[0]
1172  chain1 = d[1]
1173  resid2 = d[2]
1174  chain2 = d[3]
1175  fexp = d[4]
1176  epslabel = d[5]
1177 
1178  # CysteineCrossLinkData
1179 
1181  fexp,
1182  fmod_grid,
1183  omega2_grid,
1184  beta_grid)
1185 
1187  self.m,
1188  self.beta,
1189  self.sigma,
1190  self.epsilons[epslabel],
1191  self.weight,
1192  crossdata,
1193  ccldata)
1194 
1195  failed = False
1196  if not self.cbeta:
1197  p1 = None
1198  p2 = None
1199 
1200  p1 = IMP.atom.Selection(root_hier, resolution=1,
1201  molecule=chain1, residue_index=resid1,
1202  copy_index=0)
1203  p1 = p1.get_selected_particles()
1204  if len(p1) > 0:
1205  p1 = p1[0]
1206  else:
1207  failed = True
1208 
1209  p2 = IMP.atom.Selection(root_hier, resolution=1,
1210  molecule=chain2, residue_index=resid2,
1211  copy_index=0)
1212  p2 = p2.get_selected_particles()
1213  if len(p2) > 0:
1214  p2 = p2[0]
1215  else:
1216  failed = True
1217 
1218  else:
1219  # use cbetas
1220  p1 = []
1221  p2 = []
1222  for t in range(-1, 2):
1223  p = IMP.atom.Selection(root_hier, resolution=1,
1224  molecule=chain1, copy_index=0,
1225  residue_index=resid1 + t)
1226  p = p.get_selected_particles()
1227  if len(p) == 1:
1228  p1 += p
1229  else:
1230  failed = True
1231  warnings.warn(
1232  "CysteineCrossLink: missing representation for "
1233  "residue %d of chain %s" % (resid1 + t, chain1),
1235 
1236  p = IMP.atom.Selection(root_hier, resolution=1,
1237  molecule=chain2, copy_index=0,
1238  residue_index=resid2 + t)
1239  p = p.get_selected_particles()
1240  if len(p) == 1:
1241  p2 += p
1242  else:
1243  failed = True
1244  warnings.warn(
1245  "CysteineCrossLink: missing representation for "
1246  "residue %d of chain %s" % (resid2 + t, chain2),
1248 
1249  if not self.cbeta:
1250  if (p1 is not None and p2 is not None):
1251  ccl.add_contribution(p1, p2)
1252  d1 = IMP.core.XYZ(p1)
1253  d2 = IMP.core.XYZ(p2)
1254 
1255  print("Distance_" + str(resid1) + "_" + chain1 + ":" +
1256  str(resid2) + "_" + chain2,
1257  IMP.core.get_distance(d1, d2))
1258 
1259  else:
1260  if (len(p1) == 3 and len(p2) == 3):
1261  p11n = p1[0].get_name()
1262  p12n = p1[1].get_name()
1263  p13n = p1[2].get_name()
1264  p21n = p2[0].get_name()
1265  p22n = p2[1].get_name()
1266  p23n = p2[2].get_name()
1267 
1268  print("CysteineCrossLink: generating CB cysteine "
1269  "cross-link restraint between")
1270  print("CysteineCrossLink: residue %d of chain %s and "
1271  "residue %d of chain %s"
1272  % (resid1, chain1, resid2, chain2))
1273  print("CysteineCrossLink: between particles %s %s %s and "
1274  "%s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
1275 
1276  ccl.add_contribution(p1, p2)
1277 
1278  if not failed:
1279  self.rs.add_restraint(ccl)
1280  ccl.set_name("CysteineCrossLink_" + str(resid1)
1281  + "_" + chain1 + ":" + str(resid2) + "_" + chain2)
1282 
1284  self.weight.get_particle()
1285  ).set_weights_are_optimized(weightissampled)
1286 
1287  def set_label(self, label):
1288  self.label = label
1289 
1290  def add_to_model(self):
1291  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs, add_to_rmf=True)
1292 
1293  def get_particles_to_sample(self):
1294  ps = {}
1295  if self.epsilonissampled:
1296  for eps in self.epsilons.keys():
1297  ps["Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps +
1298  "_" + self.label] = ([self.epsilons[eps]],
1299  self.epsilonmaxtrans)
1300  if self.betaissampled:
1301  ps["Nuisances_CysteineCrossLinkRestraint_beta_" +
1302  self.label] = ([self.beta], self.betamaxtrans)
1303  if self.weightissampled:
1304  ps["Weights_CysteineCrossLinkRestraint_" +
1305  self.label] = ([self.weight], self.weightmaxtrans)
1306  if self.sigmaissampled:
1307  ps["Nuisances_CysteineCrossLinkRestraint_" +
1308  self.label] = ([self.sigma], self.sigmamaxtrans)
1309  return ps
1310 
1311  def set_output_level(self, level="low"):
1312  # this might be "low" or "high"
1313  self.outputlevel = level
1314 
1315  def get_restraint(self):
1316  return self.rs
1317 
1318  def get_sigma(self):
1319  return self.sigma
1320 
1321  def get_output(self):
1322  output = {}
1323  score = self.rs.unprotected_evaluate(None)
1324  output["_TotalScore"] = str(score)
1325  output["CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
1326  output["CysteineCrossLinkRestraint_sigma_" +
1327  self.label] = str(self.sigma.get_scale())
1328  for eps in self.epsilons.keys():
1329  output["CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
1330  self.label] = str(self.epsilons[eps].get_scale())
1331  output["CysteineCrossLinkRestraint_beta_" +
1332  self.label] = str(self.beta.get_scale())
1333  for n in range(self.weight.get_number_of_weights()):
1334  output["CysteineCrossLinkRestraint_weights_" +
1335  str(n) + "_" + self.label] = str(self.weight.get_weight(n))
1336 
1337  if self.outputlevel == "high":
1338  for rst in self.rs.get_restraints():
1339  cclr = IMP.isd.CysteineCrossLinkRestraint.get_from(rst)
1340  output["CysteineCrossLinkRestraint_Total_Frequency_" +
1341  cclr.get_name() +
1342  "_" + self.label] = cclr.get_model_frequency()
1343  output["CysteineCrossLinkRestraint_Standard_Error_" +
1344  cclr.get_name() + "_" +
1345  self.label] = cclr.get_standard_error()
1346  return output
1347 
1348 
1349 class DisulfideCrossLinkRestraint(object):
1350  def __init__(self, representation_or_hier,
1351  selection_tuple1,
1352  selection_tuple2,
1353  length=6.5,
1354  resolution=1,
1355  slope=0.01,
1356  label="None"):
1357 
1358  self.m = representation_or_hier.get_model()
1359  ps1 = IMP.pmi.tools.select_by_tuple_2(representation_or_hier,
1360  selection_tuple1,
1361  resolution=resolution)
1362  ps2 = IMP.pmi.tools.select_by_tuple_2(representation_or_hier,
1363  selection_tuple2,
1364  resolution=resolution)
1365 
1366  self.rs = IMP.RestraintSet(self.m, 'likelihood')
1367  self.rslin = IMP.RestraintSet(self.m, 'linear_dummy_restraints')
1368 
1369  # dummy linear restraint used for Chimera display
1370  self.linear = IMP.core.Linear(0, 0.0)
1371  self.linear.set_slope(0.0)
1372  dps2 = IMP.core.DistancePairScore(self.linear)
1373 
1374  self.label = label
1375  self.psi_dictionary = {}
1376  self.sigma_dictionary = {}
1377  self.psi_is_sampled = False
1378  self.sigma_is_sampled = False
1379  self.xl = {}
1380 
1381  if len(ps1) > 1 or len(ps1) == 0:
1382  raise ValueError(
1383  "DisulfideBondRestraint: ERROR> first selection pattern "
1384  "selects multiple particles or sero particles")
1385 
1386  if len(ps2) > 1 or len(ps2) == 0:
1387  raise ValueError(
1388  "DisulfideBondRestraint: ERROR> second selection pattern "
1389  "selects multiple particles or sero particles")
1390 
1391  p1 = ps1[0]
1392  p2 = ps2[0]
1393 
1394  sigma = self.create_sigma("SIGMA_DISULFIDE_BOND")
1395  psi = self.create_psi("PSI_DISULFIDE_BOND")
1396 
1397  p1i = p1.get_index()
1398  p2i = p2.get_index()
1399  si = sigma.get_particle().get_index()
1400  psii = psi.get_particle().get_index()
1401 
1402  dr = IMP.isd.CrossLinkMSRestraint(self.m, length, slope)
1403 
1404  dr.add_contribution((p1i, p2i), (si, si), psii)
1405 
1406  if p1i != p2i:
1407  pr = IMP.core.PairRestraint(self.m, dps2, (p1i, p2i))
1408  pr.set_name("DISULFIDE_BOND_"+self.label)
1409  self.rslin.add_restraint(pr)
1410 
1411  lw = IMP.isd.LogWrapper([dr], 1.0)
1412  self.rs.add_restraint(lw)
1413 
1414  self.xl["Particle1"] = p1
1415  self.xl["Particle2"] = p2
1416  self.xl["Sigma"] = sigma
1417  self.xl["Psi"] = psi
1418 
1419  def add_to_model(self):
1420  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs, add_to_rmf=True)
1422  self.m, self.rslin, add_to_rmf=True)
1423 
1424  def get_hierarchies(self):
1425  return self.prot
1426 
1427  def get_restraint_sets(self):
1428  return self.rs
1429 
1430  def get_restraint(self):
1431  return self.rs
1432 
1433  def get_restraint_for_rmf(self):
1434  return self.rslin
1435 
1436  def get_restraints(self):
1437  rlist = []
1438  for r in self.rs.get_restraints():
1439  rlist.append(IMP.core.PairRestraint.get_from(r))
1440  return rlist
1441 
1442  def set_psi_is_sampled(self, is_sampled=True):
1443  self.psi_is_sampled = is_sampled
1444 
1445  def set_sigma_is_sampled(self, is_sampled=True):
1446  self.sigma_is_sampled = is_sampled
1447 
1448  def create_sigma(self, name):
1449  ''' a nuisance on the structural uncertainty '''
1450  if name in self.sigma_dictionary:
1451  return self.sigma_dictionary[name][0]
1452 
1453  sigmainit = 1.0
1454  sigmaminnuis = 0.0000001
1455  sigmamaxnuis = 1000.0
1456  sigmatrans = 0.5
1457  sigma = IMP.pmi.tools.SetupNuisance(
1458  self.m, sigmainit, sigmaminnuis, sigmamaxnuis,
1459  self.sigma_is_sampled).get_particle()
1460  self.sigma_dictionary[name] = (
1461  sigma,
1462  sigmatrans,
1463  self.sigma_is_sampled)
1464 
1465  return sigma
1466 
1467  def create_psi(self, name):
1468  ''' a nuisance on the inconsistency '''
1469  if name in self.psi_dictionary:
1470  return self.psi_dictionary[name][0]
1471 
1472  psiinit = 0.001
1473  psiminnuis = 0.0000001
1474  psimaxnuis = 0.4999999
1475  psitrans = 0.1
1476  psi = IMP.pmi.tools.SetupNuisance(self.m, psiinit,
1477  psiminnuis, psimaxnuis,
1478  self.psi_is_sampled).get_particle()
1479  self.psi_dictionary[name] = (
1480  psi,
1481  psitrans,
1482  self.psi_is_sampled)
1483 
1484  return psi
1485 
1486  def get_output(self):
1487  output = {}
1488  score = self.rs.unprotected_evaluate(None)
1489  output["_TotalScore"] = str(score)
1490  output["DisulfideBondRestraint_Data_Score_" + self.label] = str(score)
1491  output["DisulfideBondRestraint_Linear_Score_" +
1492  self.label] = self.rslin.unprotected_evaluate(None)
1493  return output
1494 
1495  def get_particles_to_sample(self):
1496  raise NotImplementedError(" ")
def set_sigma_is_sampled
Switch on/off the sampling of sigma particles.
Setup cross-link distance restraints from mass spectrometry data.
Add weights to a particle.
Definition: Weight.h:29
A restraint for ambiguous cross-linking MS data and multiple state approach.
def get_best_stats
return the probability, best distance, two coords, and possibly the psi for each xl ...
Various classes to hold sets of particles.
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
Miscellaneous utilities.
Definition: tools.py:1
Restrain atom pairs based on a set of crosslinks.
Add scale parameter to particle.
Definition: Scale.h:24
The type of an atom.
log
Definition: log.py:1
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:749
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
Uniform distribution with harmonic boundaries.
Definition: UniformPrior.h:22
Object used to hold a set of restraints.
Definition: RestraintSet.h:41
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:86
Classes to handle different kinds of restraints.
A RestraintSet subclass to track cross-link metadata.
Warning related to handling of structures.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def select_by_tuple_2
New tuple format: molname OR (start,stop,molname,copynum,statenum) Copy and state are optional...
Definition: tools.py:438
A decorator for a particle representing an atom.
Definition: atom/Atom.h:238
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:92
def get_movers
Get all need data to construct a mover in IMP.pmi.dof class.
def get_particle_pairs
Get a list of tuples containing the particle pairs.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
int get_state_index(Hierarchy h)
Walk up the hierarchy to find the current state.
Linear function
Definition: Linear.h:22
def plot_violations
Create CMM files, one for each state, of all cross-links.
Modify a set of continuous variables using a normal distribution.
Definition: NormalMover.h:20
Classes for writing output files and processing them.
Definition: output.py:1
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
def create_restraints_for_rmf
create dummy harmonic restraints for each XL but don't add to model Makes it easy to see each contrib...
The general base class for IMP exceptions.
Definition: exception.h:48
Setup cross-link distance restraints at atomic level The "atomic" aspect is that it models the partic...
Residue get_residue(Atom d, bool nothrow=false)
Return the Residue containing this atom.
Calculate the -Log of a list of restraints.
Definition: LogWrapper.h:26
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
Definition: VectorD.h:213
Applies a PairScore to a Pair.
Definition: PairRestraint.h:31
def get_output
Get the output of the restraint to be used by the IMP.pmi.output object.
def get_restraint_for_rmf
get the dummy restraints to be displayed in the rmf file
Functionality for loading, creating, manipulating and scoring atomic structures.
std::string get_chain_id(Hierarchy h)
Walk up the hierarchy to determine the chain id.
def load_nuisances_from_stat_file
Read a stat file and load all the sigmas.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A restraint for cysteine cross-linking data.