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