IMP logo
IMP Reference Guide  2.21.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_likelihood(self):
459  """Get the unweighted likelihood of the restraint"""
460  likelihood = 1
461  for restraint in self.xl_restraints:
462  likelihood *= restraint.get_probability()
463 
464  return likelihood
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_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.
def get_likelihood
Get the unweighted likelihood of the restraint.
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
Score distance between two particle centers using a harmonic function.
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:23
Score a pair of particles based on the distance between their centers.
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:196
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.