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