IMP logo
IMP Reference Guide  2.10.0
The Integrative Modeling Platform
restraints/crosslinking.py
1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling crosslinking 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 
25  """Setup cross-link distance restraints from mass spectrometry data.
26  The noise in the data and the structural uncertainty of cross-linked amino-acids
27  is inferred using Bayes theory of probability
28  @note Wraps an IMP::isd::CrossLinkMSRestraint
29  """
30  def __init__(self, representation=None,
31  root_hier=None,
32  CrossLinkDataBase=None,
33  length=10.0,
34  resolution=None,
35  slope=0.02,
36  label=None,
37  filelabel="None",
38  attributes_for_label=None,
39  weight=1.):
40  """Constructor.
41  @param representation DEPRECATED The IMP.pmi.representation.Representation
42  object that contain the molecular system
43  @param root_hier The canonical hierarchy containing all the states
44  @param CrossLinkDataBase The IMP.pmi.io.crosslink.CrossLinkDataBase
45  object that contains the cross-link dataset
46  @param length maximal cross-linker length (including the residue sidechains)
47  @param resolution what representation resolution should the cross-link
48  restraint be applied to.
49  @param slope The slope of a distance-linear scoring function that
50  funnels the score when the particles are
51  too far away. Suggested value 0.02.
52  @param label the extra text to label the restraint so that it is
53  searchable in the output
54  @param filelabel automatically generated file containing missing/included/excluded
55  cross-links will be labeled using this text
56  @param attributes_for_label
57  @param weight Weight of restraint
58  """
59 
60  use_pmi2 = True
61  if representation is not None:
62  use_pmi2 = False
63  if type(representation) != list:
64  representations = [representation]
65  else:
66  representations = representation
67  model = representations[0].prot.get_model()
68  elif root_hier is not None:
69  representations = []
70  model = root_hier.get_model()
71  else:
72  raise Exception("You must pass either representation or root_hier")
73 
74  super(CrossLinkingMassSpectrometryRestraint, self).__init__(
75  model, weight=weight, label=label)
76 
77  if CrossLinkDataBase is None:
78  raise Exception("You must pass a CrossLinkDataBase")
79  if not isinstance(CrossLinkDataBase,IMP.pmi.io.crosslink.CrossLinkDataBase):
80  raise TypeError("CrossLinkingMassSpectrometryRestraint: CrossLinkDataBase should be an IMP.pmi.io.crosslink.CrossLinkDataBase object")
81  self.CrossLinkDataBase = CrossLinkDataBase
82 
83  if resolution==0 or resolution is None:
84  raise Exception("You must pass a resolution and it can't be zero")
85 
86  indb = open("included." + filelabel + ".xl.db", "w")
87  exdb = open("excluded." + filelabel + ".xl.db", "w")
88  midb = open("missing." + filelabel + ".xl.db", "w")
89 
90  self.rs.set_name(self.rs.get_name() + "_Data")
91  self.rspsi = self._create_restraint_set("PriorPsi")
92  self.rssig = self._create_restraint_set("PriorSig")
93  self.rslin = self._create_restraint_set("Linear")
94 
95  # dummy linear restraint used for Chimera display
96  self.linear = IMP.core.Linear(0, 0.0)
97  self.linear.set_slope(0.0)
98  dps2 = IMP.core.DistancePairScore(self.linear)
99 
100  self.psi_is_sampled = True
101  self.sigma_is_sampled = True
102  self.psi_dictionary={}
103  self.sigma_dictionary={}
104  self.xl_list=[]
105  self.outputlevel = "low"
106 
107  restraints = []
108 
109  xl_groups = [p.get_cross_link_group(self)
110  for p, state in IMP.pmi.tools._all_protocol_outputs(
111  representations, root_hier)]
112 
113  # if PMI2, first add all the molecule copies as clones to the database
114  if use_pmi2:
115  copies_to_add = defaultdict(int)
116  print('gathering copies')
117  for xlid in self.CrossLinkDataBase.xlid_iterator():
118  for xl in self.CrossLinkDataBase[xlid]:
119  r1 = xl[self.CrossLinkDataBase.residue1_key]
120  c1 = xl[self.CrossLinkDataBase.protein1_key]
121  r2 = xl[self.CrossLinkDataBase.residue2_key]
122  c2 = xl[self.CrossLinkDataBase.protein2_key]
123  for c,r in ((c1,r1),(c2,r2)):
124  if c in copies_to_add:
125  continue
126  sel = IMP.atom.Selection(root_hier,
127  state_index=0,
128  molecule=c,
129  residue_index=r,
130  resolution=resolution).get_selected_particles()
131  if len(sel)>0:
132  copies_to_add[c] = len(sel)-1
133  print(copies_to_add)
134  for molname in copies_to_add:
135  if copies_to_add[molname]==0:
136  continue
137  fo1 = IMP.pmi.io.crosslink.FilterOperator(self.CrossLinkDataBase.protein1_key,operator.eq,molname)
138  self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein1_key,molname+'.0',fo1)
139  fo2 = IMP.pmi.io.crosslink.FilterOperator(self.CrossLinkDataBase.protein2_key,operator.eq,molname)
140  self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein2_key,molname+'.0',fo2)
141  for ncopy in range(copies_to_add[molname]):
142  self.CrossLinkDataBase.clone_protein('%s.0'%molname,'%s.%i'%(molname,ncopy+1))
143  print('done pmi2 prelims')
144 
145  for xlid in self.CrossLinkDataBase.xlid_iterator():
146  new_contribution=True
147  for xl in self.CrossLinkDataBase[xlid]:
148 
149  r1 = xl[self.CrossLinkDataBase.residue1_key]
150  c1 = xl[self.CrossLinkDataBase.protein1_key]
151  r2 = xl[self.CrossLinkDataBase.residue2_key]
152  c2 = xl[self.CrossLinkDataBase.protein2_key]
153 
154  # todo: check that offset is handled correctly
155  ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
156  group), group)
157  for p, group in
158  zip(IMP.pmi.tools._all_protocol_outputs(
159  representations, root_hier),
160  xl_groups)]
161 
162  if use_pmi2:
163  iterlist = range(len(IMP.atom.get_by_type(root_hier,IMP.atom.STATE_TYPE)))
164  else:
165  iterlist = representations
166  for nstate, r in enumerate(iterlist):
167  # loop over every state
168  xl[self.CrossLinkDataBase.state_key]=nstate
169  xl[self.CrossLinkDataBase.data_set_name_key]=self.label
170 
171  if use_pmi2:
172  name1 = c1
173  name2 = c2
174  copy1 = 0
175  copy2 = 0
176  if '.' in c1:
177  name1,copy1 = c1.split('.')
178  if '.' in c2:
179  name2,copy2 = c2.split('.')
180  ps1 = IMP.atom.Selection(root_hier,
181  state_index=nstate,
182  molecule=name1,
183  copy_index=int(copy1),
184  residue_index=r1,
185  resolution=resolution).get_selected_particles()
186  ps2 = IMP.atom.Selection(root_hier,
187  state_index=nstate,
188  molecule=name2,
189  copy_index=int(copy2),
190  residue_index=r2,
191  resolution=resolution).get_selected_particles()
192 
193  ps1 = [IMP.atom.Hierarchy(p) for p in ps1]
194  ps2 = [IMP.atom.Hierarchy(p) for p in ps2]
195  else:
196  ps1 = IMP.pmi.tools.select(
197  r,
198  resolution=resolution,
199  name=c1,
200  name_is_ambiguous=False,
201  residue=r1)
202  ps2 = IMP.pmi.tools.select(
203  r,
204  resolution=resolution,
205  name=c2,
206  name_is_ambiguous=False,
207  residue=r2)
208 
209  if len(ps1) > 1:
210  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
211  elif len(ps1) == 0:
212  print("CrossLinkingMassSpectrometryRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
213  midb.write(str(xl) + "\n")
214  continue
215 
216  if len(ps2) > 1:
217  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
218  elif len(ps2) == 0:
219  print("CrossLinkingMassSpectrometryRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
220  midb.write(str(xl) + "\n")
221  continue
222 
223  p1 = ps1[0]
224  p2 = ps2[0]
225 
226  if p1 == p2 and r1 == r2:
227  print("CrossLinkingMassSpectrometryRestraint: WARNING> same particle and same residue, skipping cross-link")
228  continue
229 
230  if new_contribution:
231  print("generating a new crosslink restraint")
232  new_contribution=False
234  self.model,
235  length,
236  slope)
237  restraints.append(dr)
238 
239 
240  if self.CrossLinkDataBase.sigma1_key not in xl.keys():
241  sigma1name="SIGMA"
242  xl[self.CrossLinkDataBase.sigma1_key]=sigma1name
243  else:
244  sigma1name=xl[self.CrossLinkDataBase.sigma1_key]
245  sigma1=self.create_sigma(sigma1name)
246 
247  if self.CrossLinkDataBase.sigma2_key not in xl.keys():
248  sigma2name="SIGMA"
249  xl[self.CrossLinkDataBase.sigma2_key]=sigma2name
250  else:
251  sigma2name=xl[self.CrossLinkDataBase.sigma2_key]
252  sigma2=self.create_sigma(sigma2name)
253 
254  if self.CrossLinkDataBase.psi_key not in xl.keys():
255  psiname="PSI"
256  xl[self.CrossLinkDataBase.psi_key]=psiname
257  else:
258  psiname=xl[self.CrossLinkDataBase.psi_key]
259 
260  psi=self.create_psi(psiname)
261 
262 
263  p1i = p1.get_particle_index()
264  xl["Particle1"]=p1
265  p2i = p2.get_particle_index()
266  xl["Particle2"]=p2
267  s1i = sigma1.get_particle().get_index()
268  xl["Particle_sigma1"]=sigma1
269  s2i = sigma2.get_particle().get_index()
270  xl["Particle_sigma2"]=sigma2
271  psii = psi.get_particle().get_index()
272  xl["Particle_psi"]=psi
273 
274  dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
275  xl["Restraint"]=dr
276 
277  print("--------------")
278  print("CrossLinkingMassSpectrometryRestraint: generating cross-link restraint between")
279  print("CrossLinkingMassSpectrometryRestraint: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
280  print("CrossLinkingMassSpectrometryRestraint: with sigma1 %s sigma2 %s psi %s" % (sigma1name, sigma2name, psiname))
281  print("CrossLinkingMassSpectrometryRestraint: between particles %s and %s" % (p1.get_name(), p2.get_name()))
282  print("==========================================\n")
283  for p, ex_xl in zip(IMP.pmi.tools._all_protocol_outputs(
284  representations, root_hier),
285  ex_xls):
286  p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
287  sigma1, sigma2, psi, ex_xl[1])
288 
289  # check if the two residues belong to the same rigid body
292  IMP.core.RigidMember(p1).get_rigid_body() ==
293  IMP.core.RigidMember(p2).get_rigid_body()):
294  xl["IntraRigidBody"]=True
295  else:
296  xl["IntraRigidBody"]=False
297 
298  xl_label=self.CrossLinkDataBase.get_short_cross_link_string(xl)
299  xl["ShortLabel"]=xl_label
300  dr.set_name(xl_label)
301 
302  if p1i != p2i:
303  pr = IMP.core.PairRestraint(self.model, dps2, (p1i, p2i))
304  pr.set_name(xl_label)
305  self.rslin.add_restraint(pr)
306 
307  self.xl_list.append(xl)
308 
309  indb.write(str(xl) + "\n")
310 
311  if len(self.xl_list) == 0:
312  raise SystemError("CrossLinkingMassSpectrometryRestraint: no crosslink was constructed")
313  self.xl_restraints = restraints
314  lw = IMP.isd.LogWrapper(restraints,1.0)
315  self.rs.add_restraint(lw)
316 
317  def __set_dataset(self, ds):
318  self.CrossLinkDataBase.dataset = ds
319  dataset = property(lambda self: self.CrossLinkDataBase.dataset,
320  __set_dataset)
321 
322  def get_hierarchies(self):
323  """ get the hierarchy """
324  return self.prot
325 
327  """ get the restraint set """
328  return self.rs
329 
330  def get_restraints(self):
331  """ get the restraints in a list """
332  return self.xl_restraints
333 
335  """ get the dummy restraints to be displayed in the rmf file """
336  return self.rslin
337 
339  """ Get a list of tuples containing the particle pairs """
340  ppairs = []
341  for i in range(len(self.pairs)):
342  p0 = self.pairs[i][0]
343  p1 = self.pairs[i][1]
344  ppairs.append((p0, p1))
345  return ppairs
346 
347  def set_output_level(self, level="low"):
348  """ Set the output level of the output """
349  self.outputlevel = level
350 
351  def set_psi_is_sampled(self, is_sampled=True):
352  """ Switch on/off the sampling of psi particles """
353  self.psi_is_sampled = is_sampled
354 
355  def set_sigma_is_sampled(self, is_sampled=True):
356  """ Switch on/off the sampling of sigma particles """
357  self.sigma_is_sampled = is_sampled
358 
359 
360  def create_sigma(self, name):
361  """ This is called internally. Creates a nuisance
362  on the structural uncertainty """
363  if name in self.sigma_dictionary:
364  return self.sigma_dictionary[name][0]
365 
366  sigmainit = 2.0
367  sigmaminnuis = 0.0000001
368  sigmamaxnuis = 1000.0
369  sigmamin = 0.01
370  sigmamax = 100.0
371  sigmatrans = 0.5
372  sigma = IMP.pmi.tools.SetupNuisance(self.model, sigmainit,
373  sigmaminnuis, sigmamaxnuis, self.sigma_is_sampled).get_particle()
374  self.sigma_dictionary[name] = (
375  sigma,
376  sigmatrans,
377  self.sigma_is_sampled)
378  self.rssig.add_restraint(
380  self.model,
381  sigma,
382  1000000000.0,
383  sigmamax,
384  sigmamin))
385  return sigma
386 
387  def create_psi(self, name):
388  """ This is called internally. Creates a nuisance
389  on the data uncertainty """
390  if name in self.psi_dictionary:
391  return self.psi_dictionary[name][0]
392 
393  psiinit=0.25
394  psiminnuis = 0.0000001
395  psimaxnuis = 0.4999999
396  psimin = 0.01
397  psimax = 0.49
398  psitrans = 0.1
399  psi = IMP.pmi.tools.SetupNuisance(self.model, psiinit,
400  psiminnuis, psimaxnuis,
401  self.psi_is_sampled).get_particle()
402  self.psi_dictionary[name] = (
403  psi,
404  psitrans,
405  self.psi_is_sampled)
406 
407  self.rspsi.add_restraint(
409  self.model,
410  psi,
411  1000000000.0,
412  psimax,
413  psimin))
414 
415  self.rspsi.add_restraint(IMP.isd.JeffreysRestraint(self.model, psi))
416  return psi
417 
418  def get_output(self):
419  """ Get the output of the restraint to be used by the IMP.pmi.output object"""
420  output = super(CrossLinkingMassSpectrometryRestraint, self).get_output()
421 
422  for xl in self.xl_list:
423 
424  xl_label=xl["ShortLabel"]
425  ln = xl["Restraint"]
426  p0 = xl["Particle1"]
427  p1 = xl["Particle2"]
428  output["CrossLinkingMassSpectrometryRestraint_Score_" +
429  xl_label] = str(-log(ln.unprotected_evaluate(None)))
430 
431  d0 = IMP.core.XYZ(p0)
432  d1 = IMP.core.XYZ(p1)
433  output["CrossLinkingMassSpectrometryRestraint_Distance_" +
434  xl_label] = str(IMP.core.get_distance(d0, d1))
435 
436 
437  for psiname in self.psi_dictionary:
438  output["CrossLinkingMassSpectrometryRestraint_Psi_" +
439  str(psiname) + self._label_suffix] = str(
440  self.psi_dictionary[psiname][0].get_scale())
441 
442  for sigmaname in self.sigma_dictionary:
443  output["CrossLinkingMassSpectrometryRestraint_Sigma_" +
444  str(sigmaname) + self._label_suffix] = str(
445  self.sigma_dictionary[sigmaname][0].get_scale())
446 
447 
448  return output
449 
450  def get_movers(self):
451  """ Get all need data to construct a mover in IMP.pmi.dof class"""
452  movers=[]
453  if self.sigma_is_sampled:
454  for sigmaname in self.sigma_dictionary:
455  mover_name="Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" + str(sigmaname) + "_" + self.label
456  particle=self.sigma_dictionary[sigmaname][0]
457  maxstep=(self.sigma_dictionary[sigmaname][1])
458  mv=IMP.core.NormalMover([particle],
459  IMP.FloatKeys([IMP.FloatKey("nuisance")]),maxstep)
460  mv.set_name(mover_name)
461  movers.append(mv)
462 
463  if self.psi_is_sampled:
464  for psiname in self.psi_dictionary:
465  mover_name="Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" + str(psiname) + "_" + self.label
466  particle=self.psi_dictionary[psiname][0]
467  maxstep=(self.psi_dictionary[psiname][1])
468  mv=IMP.core.NormalMover([particle],
469  IMP.FloatKeys([IMP.FloatKey("nuisance")]),maxstep)
470  mv.set_name(mover_name)
471  movers.append(mv)
472 
473  return movers
474 
475 
477  """ Get the particles to be sampled by the IMP.pmi.sampler object """
478  ps = {}
479  if self.sigma_is_sampled:
480  for sigmaname in self.sigma_dictionary:
481  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" +
482  str(sigmaname) + self._label_suffix] =\
483  ([self.sigma_dictionary[sigmaname][0]],
484  self.sigma_dictionary[sigmaname][1])
485 
486  if self.psi_is_sampled:
487  for psiname in self.psi_dictionary:
488  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
489  str(psiname) + self._label_suffix] =\
490  ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
491 
492  return ps
493 
494 
496  """Setup cross-link distance restraints at atomic level
497  The "atomic" aspect is that it models the particle uncertainty with a Gaussian
498  The noise in the data and the structural uncertainty of cross-linked amino-acids
499  is inferred using Bayes' theory of probability
500  @note Wraps an IMP::isd::AtomicCrossLinkMSRestraint
501  """
502  def __init__(self,
503  root_hier,
504  xldb,
505  atom_type="NZ",
506  length=10.0,
507  slope=0.01,
508  nstates=None,
509  label=None,
510  nuisances_are_optimized=True,
511  sigma_init=5.0,
512  psi_init = 0.01,
513  one_psi=True,
514  filelabel=None,
515  weight=1.):
516  """Constructor.
517  Automatically creates one "sigma" per crosslinked residue and one "psis" per pair.
518  Other nuisance options are available.
519  @note Will return an error if the data+extra_sel don't specify two particles per XL pair.
520  @param root_hier The root hierarchy on which you'll do selection
521  @param xldb CrossLinkDataBase object
522  @param atom_type Can either be "NZ" or "CA"
523  @param length The XL linker length
524  @param slope Linear term to add to the restraint function to help when far away
525  @param nstates The number of states to model. Defaults to the number of states in root.
526  @param label The output label for the restraint
527  @param nuisances_are_optimized Whether to optimize nuisances
528  @param sigma_init The initial value for all the sigmas
529  @param psi_init The initial value for all the psis
530  @param one_psi Use a single psi for all restraints (if False, creates one per XL)
531  @param filelabel automatically generated file containing missing/included/excluded
532  cross-links will be labeled using this text
533  @param weight Weight of restraint
534 
535  """
536 
537  # basic params
538  self.root = root_hier
539  rname = "AtomicXLRestraint"
540  super(AtomicCrossLinkMSRestraint, self).__init__(
541  self.root.get_model(), name="AtomicXLRestraint", label=label,
542  weight=weight)
543  self.xldb = xldb
544  self.length = length
545  self.sigma_is_sampled = nuisances_are_optimized
546  self.psi_is_sampled = nuisances_are_optimized
547 
548  if atom_type not in("NZ","CA"):
549  raise Exception("AtomicXLRestraint: atom_type must be NZ or CA")
550  self.atom_type = atom_type
551  self.nstates = nstates
552  if nstates is None:
553  self.nstates = len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE))
554  elif nstates!=len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE)):
555  print("Warning: nstates is not the same as the number of states in root")
556 
557  self.rs_psi = self._create_restraint_set("psi")
558  self.rs_sig = self._create_restraint_set("sigma")
559  self.rs_lin = self._create_restraint_set("linear")
560 
561  self.psi_dictionary = {}
562  self.sigma_dictionary = {}
563 
564  self.particles=defaultdict(set)
565  self.one_psi = one_psi
566  self.bonded_pairs = []
567  if self.one_psi:
568  print('creating a single psi for all XLs')
569  else:
570  print('creating one psi for each XL')
571 
572  # output logging file
573  if filelabel is not None:
574  indb = open("included." + filelabel + ".xl.db", "w")
575  exdb = open("excluded." + filelabel + ".xl.db", "w")
576  midb = open("missing." + filelabel + ".xl.db", "w")
577 
578 
579 
580  ### Setup nuisances
581  '''
582  # read ahead to get the number of XL's per residue
583  num_xls_per_res=defaultdict(int)
584  for unique_id in data:
585  for nstate in range(self.nstates):
586  for xl in data[unique_id]:
587  num_xls_per_res[str(xl.r1)]+=1
588  num_xls_per_res[str(xl.r2)]+=1
589 
590  # Will setup two sigmas based on promiscuity of the residue
591  sig_threshold=4
592  self.sig_low = setup_nuisance(self.m,self.rs_nuis,init_val=sigma_init,min_val=1.0,
593  max_val=100.0,is_opt=self.nuis_opt)
594  self.sig_high = setup_nuisance(self.m,self.rs_nuis,init_val=sigma_init,min_val=1.0,
595  max_val=100.0,is_opt=self.nuis_opt)
596  '''
597  self._create_sigma('sigma',sigma_init)
598  if one_psi:
599  self._create_psi('psi',psi_init)
600  else:
601  for xlid in self.xldb.xlid_iterator():
602  self._create_psi(xlid,psi_init)
603 
604  ### create all the XLs
605  xlrs=[]
606  for xlid in self.xldb.xlid_iterator():
607  # create restraint for this data point
608  if one_psi:
609  psip = self.psi_dictionary['psi'][0].get_particle_index()
610  else:
611  psip = self.psi_dictionary[unique_id][0].get_particle_index()
612  r = IMP.isd.AtomicCrossLinkMSRestraint(self.model,
613  self.length,
614  psip,
615  slope,
616  True)
617  num_contributions=0
618 
619  # add a contribution for each XL ambiguity option within each state
620  for nstate in self.nstates:
621  for xl in self.xldb[xlid]:
622  r1 = xl[self.xldb.residue1_key]
623  c1 = xl[self.xldb.protein1_key].strip()
624  r2 = xl[self.xldb.residue2_key]
625  c2 = xl[self.xldb.protein2_key].strip()
626 
627  # perform selection. these may contain multiples if Copies are used
628  ps1 = IMP.atom.Selection(self.root,
629  state_index=nstate,
630  atom_type = IMP.atom.AtomType(self.atom_type),
631  molecule=c1,
632  residue_index=r1).get_selected_particles()
633  ps2 = IMP.atom.Selection(self.root,
634  state_index=nstate,
635  atom_type = IMP.atom.AtomType(self.atom_type),
636  molecule=c2,
637  residue_index=r2).get_selected_particles()
638  if len(ps1) == 0:
639  print("AtomicXLRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
640  if filelabel is not None:
641  midb.write(str(xl) + "\n")
642  continue
643 
644  if len(ps2) == 0:
645  print("AtomicXLRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
646  if filelabel is not None:
647  midb.write(str(xl) + "\n")
648  continue
649 
650 
651  # figure out sig1 and sig2 based on num XLs
652  '''
653  num1=num_xls_per_res[str(xl.r1)]
654  num2=num_xls_per_res[str(xl.r2)]
655  if num1<sig_threshold:
656  sig1=self.sig_low
657  else:
658  sig1=self.sig_high
659  if num2<sig_threshold:
660  sig2=self.sig_low
661  else:
662  sig2=self.sig_high
663  '''
664  sig1 = self.sigma_dictionary['sigma'][0]
665  sig2 = self.sigma_dictionary['sigma'][0]
666 
667  # add each copy contribution to restraint
668  for p1,p2 in itertools.product(ps1,ps2):
669  if p1==p2:
670  continue
671  self.particles[nstate]|=set((p1,p2))
672  r.add_contribution([p1.get_index(),p2.get_index()],
673  [sig1.get_particle_index(),sig2.get_particle_index()])
674  num_contributions+=1
675 
676  if num_contributions>0:
677  print('XL:',xlid,'num contributions:',num_contributions)
678  xlrs.append(r)
679  if len(xlrs)==0:
680  raise Exception("You didn't create any XL restraints")
681  print('created',len(xlrs),'XL restraints')
682  rname = self.rs.get_name()
683  self.rs=IMP.isd.LogWrapper(xlrs, self.weight)
684  self.rs.set_name(rname)
685  self.rs.set_weight(self.weight)
686  self.restraint_sets = [self.rs] + self.restraint_sets[1:]
687 
688  def get_hierarchy(self):
689  return self.prot
690 
691  def _create_sigma(self, name,sigmainit):
692  """ This is called internally. Creates a nuisance
693  on the structural uncertainty """
694  if name in self.sigma_dictionary:
695  return self.sigma_dictionary[name][0]
696 
697  sigmaminnuis = 0.0000001
698  sigmamaxnuis = 1000.0
699  sigmamin = 0.01
700  sigmamax = 100.0
701  sigmatrans = 0.5
702  sigma = IMP.pmi.tools.SetupNuisance(self.model,
703  sigmainit,
704  sigmaminnuis,
705  sigmamaxnuis,
706  self.sigma_is_sampled).get_particle()
707  self.sigma_dictionary[name] = (
708  sigma,
709  sigmatrans,
710  self.sigma_is_sampled)
711  self.rs_sig.add_restraint(
713  self.model,
714  sigma,
715  1000000000.0,
716  sigmamax,
717  sigmamin))
718  return sigma
719 
720  def _create_psi(self, name,psiinit):
721  """ This is called internally. Creates a nuisance
722  on the data uncertainty """
723  if name in self.psi_dictionary:
724  return self.psi_dictionary[name][0]
725 
726  psiminnuis = 0.0000001
727  psimaxnuis = 0.4999999
728  psimin = 0.01
729  psimax = 0.49
730  psitrans = 0.1
731  psi = IMP.pmi.tools.SetupNuisance(self.model,
732  psiinit,
733  psiminnuis,
734  psimaxnuis,
735  self.psi_is_sampled).get_particle()
736  self.psi_dictionary[name] = (
737  psi,
738  psitrans,
739  self.psi_is_sampled)
740 
741  self.rs_psi.add_restraint(
743  self.model,
744  psi,
745  1000000000.0,
746  psimax,
747  psimin))
748 
749  self.rs_psi.add_restraint(IMP.isd.JeffreysRestraint(self.model, psi))
750  return psi
751 
753  """ create dummy harmonic restraints for each XL but don't add to model
754  Makes it easy to see each contribution to each XL in RMF
755  """
756  class MyGetRestraint(object):
757  def __init__(self,rs):
758  self.rs=rs
759  def get_restraint_for_rmf(self):
760  return self.rs
761 
762  dummy_model=IMP.Model()
763  hps = IMP.core.HarmonicDistancePairScore(self.length,1.0)
764  dummy_rs=[]
765  for nxl in range(self.get_number_of_restraints()):
766  xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
767  rs = IMP.RestraintSet(dummy_model, 'atomic_xl_'+str(nxl))
768  for ncontr in range(xl.get_number_of_contributions()):
769  ps=xl.get_contribution(ncontr)
770  dr = IMP.core.PairRestraint(hps,[self.model.get_particle(p) for p in ps],
771  'xl%i_contr%i'%(nxl,ncontr))
772  rs.add_restraint(dr)
773  dummy_rs.append(MyGetRestraint(rs))
774  return dummy_rs
775 
776 
778  """ Get the particles to be sampled by the IMP.pmi.sampler object """
779  ps = {}
780  if self.sigma_is_sampled:
781  for sigmaname in self.sigma_dictionary:
782  ps["Nuisances_AtomicCrossLinkingMSRestraint_Sigma_" +
783  str(sigmaname) + self._label_suffix] = \
784  ([self.sigma_dictionary[sigmaname][0]],
785  self.sigma_dictionary[sigmaname][1])
786  if self.psi_is_sampled:
787  for psiname in self.psi_dictionary:
788  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
789  str(psiname) + self._label_suffix] =\
790  ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
791  return ps
792 
793  def get_bonded_pairs(self):
794  return self.bonded_pairs
795 
796  def get_number_of_restraints(self):
797  return self.rs.get_number_of_restraints()
798 
799  def __repr__(self):
800  return 'XL restraint with '+str(len(self.rs.get_restraint(0).get_number_of_restraints())) \
801  + ' data points'
802 
803  def load_nuisances_from_stat_file(self,in_fn,nframe):
804  """Read a stat file and load all the sigmas.
805  This is potentially quite stupid.
806  It's also a hack since the sigmas should be stored in the RMF file.
807  Also, requires one sigma and one psi for ALL XLs.
808  """
809  import subprocess
810  sig_val = float(subprocess.check_output(["process_output.py","-f",in_fn,
811  "-s","AtomicXLRestraint_sigma"]).split('\n>')[1+nframe])
812  psi_val = float(subprocess.check_output(["process_output.py","-f",in_fn,
813  "-s","AtomicXLRestraint_psi"]).split('\n>')[1+nframe])
814  for nxl in range(self.rs.get_number_of_restraints()):
815  xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
816  psip = xl.get_psi()
817  IMP.isd.Scale(self.model,psip).set_scale(psi_val)
818  for contr in range(xl.get_number_of_contributions()):
819  sig1,sig2=xl.get_contribution_sigmas(contr)
820  IMP.isd.Scale(self.model,sig1).set_scale(sig_val)
821 
822  print('loaded nuisances from file')
823 
824  def plot_violations(self,out_prefix,
825  max_prob_for_violation=0.1,
826  min_dist_for_violation=1e9,
827  coarsen=False,
828  limit_to_chains=None,
829  exclude_chains=''):
830  """Create CMM files, one for each state, of all crosslinks.
831  will draw in GREEN if non-violated in all states (or if only one state)
832  will draw in PURPLE if non-violated only in a subset of states (draws nothing elsewhere)
833  will draw in RED in ALL states if all violated
834  (if only one state, you'll only see green and red)
835 
836  @param out_prefix Output xlink files prefix
837  @param max_prob_for_violation It's a violation if the probability is below this
838  @param min_dist_for_violation It's a violation if the min dist is above this
839  @param coarsen Use CA positions
840  @param limit_to_chains Try to visualize just these chains
841  @param exclude_to_chains Try to NOT visualize these chains
842  """
843  print('going to calculate violations and plot CMM files')
844  all_stats = self.get_best_stats()
845  all_dists = [s["low_dist"] for s in all_stats]
846 
847  # prepare one output file per state
848  out_fns=[]
849  out_nvs=[]
850  state_info=[]
851  cmds = defaultdict(set)
852  for nstate in range(self.nstates):
853  outf=open(out_prefix+str(nstate)+'.cmm','w')
854  outf.write('<marker_set name="xlinks_state%i"> \n' % nstate)
855  out_fns.append(outf)
856  out_nvs.append(0)
857  print('will limit to',limit_to_chains)
858  print('will exclude',exclude_chains)
859  state_info.append(self.get_best_stats(nstate,
860  limit_to_chains,
861  exclude_chains))
862 
863  for nxl in range(self.rs.get_number_of_restraints()):
864  # for this XL, check which states passed
865  npass=[]
866  nviol=[]
867  for nstate in range(self.nstates):
868  prob = state_info[nstate][nxl]["prob"]
869  low_dist = state_info[nstate][nxl]["low_dist"]
870  if prob<max_prob_for_violation or low_dist>min_dist_for_violation:
871  nviol.append(nstate)
872  else:
873  npass.append(nstate)
874 
875  # special coloring when all pass or all fail
876  all_pass=False
877  all_viol=False
878  if len(npass)==self.nstates:
879  all_pass=True
880  elif len(nviol)==self.nstates:
881  all_viol=True
882 
883  # finally, color based on above info
884  print(nxl,'state dists:',[state_info[nstate][nxl]["low_dist"] for nstate in range(self.nstates)],
885  'viol states:',nviol,'all viol?',all_viol)
886  for nstate in range(self.nstates):
887  if all_pass:
888  r=0.365; g=0.933; b=0.365;
889  continue
890  elif all_viol:
891  r=0.980; g=0.302; b=0.247;
892  continue
893  else:
894  if nstate in nviol:
895  continue
896  else:
897  #r=0.9; g=0.34; b=0.9;
898  r=0.365; g=0.933; b=0.365;
899  # now only showing if UNIQUELY PASSING in this state
900  pp = state_info[nstate][nxl]["low_pp"]
901  c1=IMP.core.XYZ(self.model,pp[0]).get_coordinates()
902  c2=IMP.core.XYZ(self.model,pp[1]).get_coordinates()
903 
904  r1 = IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[0])).get_index()
905  ch1 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[0]))
906  r2 = IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[0])).get_index()
907  ch2 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[0]))
908 
909  cmds[nstate].add((ch1,r1))
910  cmds[nstate].add((ch2,r2))
911 
912  outf = out_fns[nstate]
913  nv = out_nvs[nstate]
914  outf.write('<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
915  'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,c1[0],c1[1],c1[2],r,g,b))
916  outf.write('<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
917  'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv+1,c2[0],c2[1],c2[2],r,g,b))
918  outf.write('<link id1= "%d" id2="%d" radius="0.8" '
919  'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,nv+1,r,g,b))
920  out_nvs[nstate]+=2
921 
922  for nstate in range(self.nstates):
923  out_fns[nstate].write('</marker_set>\n')
924  out_fns[nstate].close()
925  cmd = ''
926  for ch,r in cmds[nstate]:
927  cmd+='#%i:%i.%s '%(nstate,r,ch)
928  print(cmd)
929 
930  return all_dists
931  def _get_contribution_info(self,xl,ncontr,use_CA=False):
932  """Return the particles at that contribution. If requested will return CA's instead"""
933  idx1=xl.get_contribution(ncontr)[0]
934  idx2=xl.get_contribution(ncontr)[1]
935  if use_CA:
936  idx1 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,idx1)),
937  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
938  idx2 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,idx2)),
939  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
940  dist = IMP.algebra.get_distance(IMP.core.XYZ(self.model,idx1).get_coordinates(),
941  IMP.core.XYZ(self.model,idx2).get_coordinates())
942  return idx1,idx2,dist
943 
944  def get_best_stats(self,limit_to_state=None,limit_to_chains=None,exclude_chains='',use_CA=False):
945  ''' return the probability, best distance, two coords, and possibly the psi for each xl
946  @param limit_to_state Only examine contributions from one state
947  @param limit_to_chains Returns the particles for certain "easy to visualize" chains
948  @param exclude_chains Even if you limit, don't let one end be in this list.
949  Only works if you also limit chains
950  @param use_CA Limit to CA particles
951  '''
952  ret = []
953  for nxl in range(self.rs.get_number_of_restraints()):
954  this_info = {}
955  xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
956  low_dist=1e6
957  low_contr = None
958  low_pp = None
959  state_contrs=[]
960  low_pp_lim = None
961  low_dist_lim = 1e6
962  for contr in range(xl.get_number_of_contributions()):
963  pp = xl.get_contribution(contr)
964  if use_CA:
965  idx1 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[0])),
966  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
967  idx2 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[1])),
968  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
969  pp = [idx1,idx2]
970  if limit_to_state is not None:
971  nstate = IMP.atom.get_state_index(IMP.atom.Atom(self.model,pp[0]))
972  if nstate!=limit_to_state:
973  continue
974  state_contrs.append(contr)
975  dist = IMP.core.get_distance(IMP.core.XYZ(self.model,pp[0]),
976  IMP.core.XYZ(self.model,pp[1]))
977  if limit_to_chains is not None:
978  c1 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[0]))
979  c2 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[1]))
980  if (c1 in limit_to_chains or c2 in limit_to_chains) and (
981  c1 not in exclude_chains and c2 not in exclude_chains):
982  if dist<low_dist_lim:
983  low_dist_lim = dist
984  low_pp_lim = pp
985  if dist<low_dist:
986  low_dist = dist
987  low_contr = contr
988  low_pp = pp
989  if limit_to_state is not None:
990  this_info["prob"] = xl.evaluate_for_contributions(state_contrs,None)
991  else:
992  this_info["prob"] = xl.unprotected_evaluate(None)
993  if limit_to_chains is not None:
994  this_info["low_pp"] = low_pp_lim
995  else:
996  this_info["low_pp"] = low_pp
997 
998  this_info["low_dist"] = low_dist
999  if not self.one_psi:
1000  pval = IMP.isd.Scale(self.model,xl.get_psi()).get_scale()
1001  this_info["psi"] = pval
1002  ret.append(this_info)
1003  return ret
1004 
1005  def print_stats(self):
1006  #print("XL restraint statistics\n<num> <prob> <bestdist> <sig1> <sig2> <is_viol>")
1007  stats = self.get_best_stats()
1008  for nxl,s in enumerate(stats):
1009  #print('%i %.4f %.4f %.4f %.4f %i'%(nxl,prob,low_dist,sig1,sig2,is_viol))
1010  print(s["low_dist"])
1011 
1012 
1013  def get_output(self):
1014  output = super(AtomicCrossLinkMSRestraint, self).get_output()
1015 
1016  ### HACK to make it easier to see the few sigmas
1017  #output["AtomicXLRestraint_sigma"] = self.sigma.get_scale()
1018  #if self.one_psi:
1019  # output["AtomicXLRestraint_psi"] = self.psi.get_scale()
1020  ######
1021 
1022  # count distances above length
1023  bad_count=0
1024  stats = self.get_best_stats()
1025  for nxl,s in enumerate(stats):
1026  if s['low_dist']>20.0:
1027  bad_count+=1
1028  output["AtomicXLRestraint_%i_%s"%(nxl,"Prob")]=str(s['prob'])
1029  output["AtomicXLRestraint_%i_%s"%(nxl,"BestDist")]=str(s['low_dist'])
1030  if not self.one_psi:
1031  output["AtomicXLRestraint_%i_%s"%(nxl,"psi")]=str(s['psi'])
1032  output["AtomicXLRestraint_NumViol"] = str(bad_count)
1033  return output
1034 
1035 
1036 @IMP.deprecated_object("2.5", "Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1037 class ConnectivityCrossLinkMS(object):
1038  """This restraint allows ambiguous crosslinking between multiple copies
1039  it is a variant of the SimplifiedCrossLinkMS
1040  """
1041  def __init__(
1042  self,
1043  representation,
1044  restraints_file,
1045  expdistance,
1046  strength=None,
1047  resolution=None):
1048 
1049  self.m = representation.prot.get_model()
1050  self.rs = IMP.RestraintSet(self.m, 'data')
1051  self.weight = 1.0
1052  self.label = "None"
1053  self.pairs = []
1054 
1055  self.outputlevel = "low"
1056  self.expdistance = expdistance
1057  self.strength = strength
1058 
1059  fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1060 
1061  for line in fl:
1062 
1063  tokens = line.split()
1064  # skip character
1065  if (tokens[0] == "#"):
1066  continue
1067  r1 = int(tokens[2])
1068  c1 = tokens[0]
1069  r2 = int(tokens[3])
1070  c2 = tokens[1]
1071 
1072  ps1 = IMP.pmi.tools.select(
1073  representation,
1074  resolution=resolution,
1075  name=c1,
1076  name_is_ambiguous=True,
1077  residue=r1)
1078  hrc1 = [representation.hier_db.particle_to_name[p] for p in ps1]
1079  if len(ps1) == 0:
1080  print("ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1081  continue
1082 
1083  ps2 = IMP.pmi.tools.select(
1084  representation,
1085  resolution=resolution,
1086  name=c2,
1087  name_is_ambiguous=True,
1088  residue=r2)
1089  hrc2 = [representation.hier_db.particle_to_name[p] for p in ps2]
1090  if len(ps2) == 0:
1091  print("ConnectivityCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1092  continue
1093 
1094  s1 = IMP.atom.Selection(ps1)
1095  s2 = IMP.atom.Selection(ps2)
1096  s1_inds = s1.get_selected_particle_indexes()
1097  s2_inds = s2.get_selected_particle_indexes()
1098  if not set(s1_inds).isdisjoint(set(s2_inds)):
1099  print("ConnectivityCrossLinkMS: WARNING> %s %s and %s %s "
1100  "select the same bead(s) - ignoring" % (c1, r1, c2, r2))
1101  continue
1102 
1103  # calculate the radii to estimate the slope of the restraint
1104  if self.strength is None:
1105  rad1 = 0
1106  rad2 = 0
1107  for p in ps1:
1108  rad1 += IMP.pmi.Uncertainty(p).get_uncertainty()
1109 
1110  for p in ps2:
1111  rad2 += IMP.pmi.Uncertainty(p).get_uncertainty()
1112 
1113  rad1 = rad1 / len(ps1)
1114  rad2 = rad2 / len(ps2)
1115 
1116  self.strength = 1 / (rad1 ** 2 + rad2 ** 2)
1117 
1118  sels = [s1, s2]
1120  sels,
1121  self.expdistance,
1122  self.strength)
1123 
1124  self.rs.add_restraint(cr)
1125  self.pairs.append((ps1, hrc1, c1, r1, ps2, hrc2, c2, r2, cr))
1126 
1127  def plot_restraint(
1128  self,
1129  uncertainty1,
1130  uncertainty2,
1131  maxdist=50,
1132  npoints=10):
1133 
1134  p1 = IMP.Particle(self.m)
1135  p2 = IMP.Particle(self.m)
1138  d1.set_radius(uncertainty1)
1139  d2.set_radius(uncertainty2)
1140  s1 = IMP.atom.Selection(p1)
1141  s2 = IMP.atom.Selection(p2)
1142  sels = [s1, s2]
1143  strength = 1 / (uncertainty1 ** 2 + uncertainty2 ** 2)
1145  sels,
1146  self.expdistance,
1147  strength)
1148  dists = []
1149  scores = []
1150  for i in range(npoints):
1151  d2.set_coordinates(
1152  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
1153  dists.append(IMP.core.get_distance(d1, d2))
1154  scores.append(cr.unprotected_evaluate(None))
1155  IMP.pmi.output.plot_xy_data(dists, scores)
1156 
1157  def set_label(self, label):
1158  self.label = label
1159  self.rs.set_name(label)
1160  for r in self.rs.get_restraints():
1161  r.set_name(label)
1162 
1163  def add_to_model(self):
1164  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1165 
1166  def get_restraint_sets(self):
1167  return self.rs
1168 
1169  def get_restraint(self):
1170  return self.rs
1171 
1172  def set_output_level(self, level="low"):
1173  # this might be "low" or "high"
1174  self.outputlevel = level
1175 
1176  def set_weight(self, weight):
1177  self.weight = weight
1178  self.rs.set_weight(weight)
1179 
1180  def get_output(self):
1181  # content of the crosslink database pairs
1182  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
1183  output = {}
1184  score = self.weight * self.rs.unprotected_evaluate(None)
1185  output["_TotalScore"] = str(score)
1186  output["ConnectivityCrossLinkMS_Score_" + self.label] = str(score)
1187  for n, p in enumerate(self.pairs):
1188 
1189  ps1 = p[0]
1190  hrc1 = p[1]
1191  c1 = p[2]
1192  r1 = p[3]
1193  ps2 = p[4]
1194  hrc2 = p[5]
1195  c2 = p[6]
1196  r2 = p[7]
1197  cr = p[8]
1198  for n1, p1 in enumerate(ps1):
1199  name1 = hrc1[n1]
1200 
1201  for n2, p2 in enumerate(ps2):
1202  name2 = hrc2[n2]
1203  d1 = IMP.core.XYZR(p1)
1204  d2 = IMP.core.XYZR(p2)
1205  label = str(r1) + ":" + name1 + "_" + str(r2) + ":" + name2
1206  output["ConnectivityCrossLinkMS_Distance_" +
1207  label] = str(IMP.core.get_distance(d1, d2))
1208 
1209  label = str(r1) + ":" + c1 + "_" + str(r2) + ":" + c2
1210  output["ConnectivityCrossLinkMS_Score_" +
1211  label] = str(self.weight * cr.unprotected_evaluate(None))
1212 
1213  return output
1214 
1215 @IMP.deprecated_object("2.5", "Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1216 class SimplifiedCrossLinkMS(object):
1217  def __init__(
1218  self,
1219  representation,
1220  restraints_file,
1221  expdistance,
1222  strength,
1223  resolution=None,
1224  columnmapping=None,
1225  truncated=True,
1226  spheredistancepairscore=True):
1227  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2
1228  # by default column 0 = protein1; column 1 = protein2; column 2 =
1229  # residue1; column 3 = residue2
1230 
1231  if columnmapping is None:
1232  columnmapping = {}
1233  columnmapping["Protein1"] = 0
1234  columnmapping["Protein2"] = 1
1235  columnmapping["Residue1"] = 2
1236  columnmapping["Residue2"] = 3
1237 
1238  self.m = representation.prot.get_model()
1239  self.rs = IMP.RestraintSet(self.m, 'data')
1240  self.weight = 1.0
1241  self.label = "None"
1242  self.pairs = []
1243  self.already_added_pairs = {}
1244 
1245  self.outputlevel = "low"
1246  self.expdistance = expdistance
1247  self.strength = strength
1248  self.truncated = truncated
1249  self.spheredistancepairscore = spheredistancepairscore
1250 
1251  # fill the cross-linker pmfs
1252  # to accelerate the init the list listofxlinkertypes might contain only
1253  # yht needed crosslinks
1254  protein1 = columnmapping["Protein1"]
1255  protein2 = columnmapping["Protein2"]
1256  residue1 = columnmapping["Residue1"]
1257  residue2 = columnmapping["Residue2"]
1258 
1259  fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1260 
1261  for line in fl:
1262 
1263  tokens = line.split()
1264  # skip character
1265  if (tokens[0] == "#"):
1266  continue
1267  r1 = int(tokens[residue1])
1268  c1 = tokens[protein1]
1269  r2 = int(tokens[residue2])
1270  c2 = tokens[protein2]
1271 
1272  ps1 = IMP.pmi.tools.select(
1273  representation,
1274  resolution=resolution,
1275  name=c1,
1276  name_is_ambiguous=False,
1277  residue=r1)
1278  ps2 = IMP.pmi.tools.select(
1279  representation,
1280  resolution=resolution,
1281  name=c2,
1282  name_is_ambiguous=False,
1283  residue=r2)
1284 
1285  if len(ps1) > 1:
1286  raise ValueError(
1287  "residue %d of chain %s selects multiple particles; "
1288  "particles are: %s"
1289  % (r1, c1, "".join(p.get_name() for p in ps1)))
1290  elif len(ps1) == 0:
1291  print("SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1292  continue
1293 
1294  if len(ps2) > 1:
1295  raise ValueError(
1296  "residue %d of chain %s selects multiple particles; "
1297  "particles are: %s"
1298  % (r2, c2, "".join(p.get_name() for p in ps2)))
1299  elif len(ps2) == 0:
1300  print("SimplifiedCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1301  continue
1302 
1303  p1 = ps1[0]
1304  p2 = ps2[0]
1305 
1306  if (p1, p2) in self.already_added_pairs:
1307  dr = self.already_added_pairs[(p1, p2)]
1308  weight = dr.get_weight()
1309  dr.set_weight(weight + 1.0)
1310  print("SimplifiedCrossLinkMS> crosslink %d %s %d %s was already found, adding 1.0 to the weight, weight is now %d" % (r1, c1, r2, c2, weight + 1.0))
1311  continue
1312 
1313  else:
1314 
1315  if self.truncated:
1316  limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
1318  self.expdistance,
1319  self.strength,
1320  self.expdistance +
1321  15.,
1322  limit)
1323  else:
1325  self.expdistance,
1326  self.strength)
1327  if self.spheredistancepairscore:
1329  else:
1330  df = IMP.core.DistancePairScore(hub)
1331  dr = IMP.core.PairRestraint(df, (p1, p2))
1332  dr.set_name(c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2))
1333 
1334  self.rs.add_restraint(dr)
1335  self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
1336  self.already_added_pairs[(p1, p2)] = dr
1337  self.already_added_pairs[(p2, p1)] = dr
1338 
1339  def set_label(self, label):
1340  self.label = label
1341 
1342  def add_to_model(self):
1343  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1344 
1345  def get_restraint_sets(self):
1346  return self.rs
1347 
1348  def get_restraint(self):
1349  return self.rs
1350 
1351  def get_restraints(self):
1352  rlist = []
1353  for r in self.rs.get_restraints():
1354  rlist.append(IMP.core.PairRestraint.get_from(r))
1355  return rlist
1356 
1357  def get_particle_pairs(self):
1358  ppairs = []
1359  for i in range(len(self.pairs)):
1360  p0 = self.pairs[i][0]
1361  p1 = self.pairs[i][1]
1362  ppairs.append((p0, p1))
1363  return ppairs
1364 
1365  def set_output_level(self, level="low"):
1366  # this might be "low" or "high"
1367  self.outputlevel = level
1368 
1369  def set_weight(self, weight):
1370  self.weight = weight
1371  self.rs.set_weight(weight)
1372 
1373  def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
1374 
1375  p1 = IMP.Particle(self.m)
1376  p2 = IMP.Particle(self.m)
1379  d1.set_radius(radius1)
1380  d2.set_radius(radius2)
1381  if self.truncated:
1382  limit = self.strength * (self.expdistance + 15) ** 2 + 10.0
1384  self.expdistance,
1385  self.strength,
1386  self.expdistance +
1387  15.,
1388  limit)
1389  else:
1391  self.expdistance,
1392  self.strength)
1394  dr = IMP.core.PairRestraint(df, (p1, p2))
1395  dists = []
1396  scores = []
1397  for i in range(npoints):
1398  d2.set_coordinates(
1399  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
1400  dists.append(IMP.core.get_distance(d1, d2))
1401  scores.append(dr.unprotected_evaluate(None))
1402  IMP.pmi.output.plot_xy_data(dists, scores)
1403 
1404  def get_output(self):
1405  # content of the crosslink database pairs
1406  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
1407  output = {}
1408  score = self.weight * self.rs.unprotected_evaluate(None)
1409  output["_TotalScore"] = str(score)
1410  output["SimplifiedCrossLinkMS_Score_" + self.label] = str(score)
1411  for i in range(len(self.pairs)):
1412 
1413  p0 = self.pairs[i][0]
1414  p1 = self.pairs[i][1]
1415  crosslinker = 'standard'
1416  ln = self.pairs[i][2]
1417  resid1 = self.pairs[i][3]
1418  chain1 = self.pairs[i][4]
1419  resid2 = self.pairs[i][5]
1420  chain2 = self.pairs[i][6]
1421 
1422  label = str(resid1) + ":" + chain1 + "_" + \
1423  str(resid2) + ":" + chain2
1424  output["SimplifiedCrossLinkMS_Score_" + crosslinker + "_" +
1425  label] = str(self.weight * ln.unprotected_evaluate(None))
1426 
1427  if self.spheredistancepairscore:
1428  d0 = IMP.core.XYZR(p0)
1429  d1 = IMP.core.XYZR(p1)
1430  else:
1431  d0 = IMP.core.XYZ(p0)
1432  d1 = IMP.core.XYZ(p1)
1433  output["SimplifiedCrossLinkMS_Distance_" +
1434  label] = str(IMP.core.get_distance(d0, d1))
1435 
1436  return output
1437 
1438 #
1439 
1440 @IMP.deprecated_object("2.5", "Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1441 class SigmoidalCrossLinkMS(object):
1442  def __init__(
1443  self, representation, restraints_file, inflection, slope, amplitude,
1444  linear_slope, resolution=None, columnmapping=None, csvfile=False,
1445  filters=None, filelabel="None"):
1446  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2
1447  # by default column 0 = protein1; column 1 = protein2; column 2 = residue1; column 3 = residue2
1448  # the filters applies to the csvfile, the format is
1449  # filters=[("Field1",">|<|=|>=|<=",value),("Field2","is","String"),("Field2","in","String")]
1450 
1451  if columnmapping is None:
1452  columnmapping = {}
1453  columnmapping["Protein1"] = 0
1454  columnmapping["Protein2"] = 1
1455  columnmapping["Residue1"] = 2
1456  columnmapping["Residue2"] = 3
1457 
1458  if csvfile:
1459  # in case the file is a csv file
1460  # columnmapping will contain the field names
1461  # that compare in the first line of the csv file
1462  db = tools.get_db_from_csv(restraints_file)
1463  else:
1464  db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1465 
1466  self.m = representation.prot.get_model()
1467  self.rs = IMP.RestraintSet(self.m, 'data')
1468  self.weight = 1.0
1469 
1470  self.label = "None"
1471  self.pairs = []
1472  self.already_added_pairs = {}
1473  self.inflection = inflection
1474  self.slope = slope
1475  self.amplitude = amplitude
1476  self.linear_slope = linear_slope
1477 
1478  self.outputlevel = "low"
1479 
1480  # fill the cross-linker pmfs
1481  # to accelerate the init the list listofxlinkertypes might contain only
1482  # yht needed crosslinks
1483  protein1 = columnmapping["Protein1"]
1484  protein2 = columnmapping["Protein2"]
1485  residue1 = columnmapping["Residue1"]
1486  residue2 = columnmapping["Residue2"]
1487 
1488  indb = open("included." + filelabel + ".xl.db", "w")
1489  exdb = open("excluded." + filelabel + ".xl.db", "w")
1490  midb = open("missing." + filelabel + ".xl.db", "w")
1491 
1492  for entry in db:
1493  if not csvfile:
1494  tokens = entry.split()
1495  # skip character
1496  if (tokens[0] == "#"):
1497  continue
1498  r1 = int(tokens[residue1])
1499  c1 = tokens[protein1]
1500  r2 = int(tokens[residue2])
1501  c2 = tokens[protein2]
1502  else:
1503  if filters is not None:
1504  if eval(tools.cross_link_db_filter_parser(filters)) == False:
1505  exdb.write(str(entry) + "\n")
1506  continue
1507  indb.write(str(entry) + "\n")
1508  r1 = int(entry[residue1])
1509  c1 = entry[protein1]
1510  r2 = int(entry[residue2])
1511  c2 = entry[protein2]
1512 
1513  ps1 = IMP.pmi.tools.select(
1514  representation,
1515  resolution=resolution,
1516  name=c1,
1517  name_is_ambiguous=False,
1518  residue=r1)
1519  ps2 = IMP.pmi.tools.select(
1520  representation,
1521  resolution=resolution,
1522  name=c2,
1523  name_is_ambiguous=False,
1524  residue=r2)
1525 
1526  if len(ps1) > 1:
1527  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1528  elif len(ps1) == 0:
1529  print("SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1530  midb.write(str(entry) + "\n")
1531  continue
1532 
1533  if len(ps2) > 1:
1534  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1535  elif len(ps2) == 0:
1536  print("SigmoidalCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1537  midb.write(str(entry) + "\n")
1538  continue
1539 
1540  p1 = ps1[0]
1541  p2 = ps2[0]
1542 
1543  if (p1, p2) in self.already_added_pairs:
1544  dr = self.already_added_pairs[(p1, p2)]
1545  weight = dr.get_weight()
1546  dr.increment_amplitude(amplitude)
1547  print("SigmoidalCrossLinkMS> crosslink %d %s %d %s was already found, adding %d to the amplitude, amplitude is now %d" % (r1, c1, r2, c2, amplitude, dr.get_amplitude()))
1548  dr.set_name(c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2)
1549  + "-ampl:" + str(dr.get_amplitude()))
1550  continue
1551 
1552  else:
1553 
1555  self.m,
1556  p1,
1557  p2,
1558  self.inflection,
1559  self.slope,
1560  self.amplitude,
1561  self.linear_slope)
1562  dr.set_name(c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2)
1563  + "-ampl:" + str(dr.get_amplitude()))
1564 
1565  self.rs.add_restraint(dr)
1566 
1567  self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
1568  self.already_added_pairs[(p1, p2)] = dr
1569  self.already_added_pairs[(p2, p1)] = dr
1570 
1571  def set_label(self, label):
1572  self.label = label
1573 
1574  def add_to_model(self):
1575  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1576 
1577  def get_restraint_sets(self):
1578  return self.rs
1579 
1580  def get_restraint(self):
1581  return self.rs
1582 
1583  def get_restraints(self):
1584  rlist = []
1585  for r in self.rs.get_restraints():
1586  rlist.append(IMP.core.PairRestraint.get_from(r))
1587  return rlist
1588 
1589  def get_particle_pairs(self):
1590  ppairs = []
1591  for i in range(len(self.pairs)):
1592  p0 = self.pairs[i][0]
1593  p1 = self.pairs[i][1]
1594  ppairs.append((p0, p1))
1595  return ppairs
1596 
1597  def set_output_level(self, level="low"):
1598  # this might be "low" or "high"
1599  self.outputlevel = level
1600 
1601  def set_weight(self, weight):
1602  self.weight = weight
1603  self.rs.set_weight(weight)
1604 
1605  def plot_restraint(self, radius1, radius2, maxdist=50, npoints=10):
1606  p1 = IMP.Particle(self.m)
1607  p2 = IMP.Particle(self.m)
1610  d1.set_radius(radius1)
1611  d2.set_radius(radius2)
1613  self.m,
1614  p1,
1615  p2,
1616  self.inflection,
1617  self.slope,
1618  self.amplitude,
1619  self.linear_slope)
1620  dists = []
1621  scores = []
1622  for i in range(npoints):
1623  d2.set_coordinates(
1624  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
1625  dists.append(IMP.core.get_distance(d1, d2))
1626  scores.append(dr.unprotected_evaluate(None))
1627  IMP.pmi.output.plot_xy_data(dists, scores)
1628 
1629  def get_output(self):
1630  # content of the crosslink database pairs
1631  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
1632  output = {}
1633  score = self.weight * self.rs.unprotected_evaluate(None)
1634  output["_TotalScore"] = str(score)
1635  output["SigmoidalCrossLinkMS_Score_" + self.label] = str(score)
1636  for i in range(len(self.pairs)):
1637 
1638  p0 = self.pairs[i][0]
1639  p1 = self.pairs[i][1]
1640  crosslinker = 'standard'
1641  ln = self.pairs[i][2]
1642  resid1 = self.pairs[i][3]
1643  chain1 = self.pairs[i][4]
1644  resid2 = self.pairs[i][5]
1645  chain2 = self.pairs[i][6]
1646 
1647  label = str(resid1) + ":" + chain1 + "_" + \
1648  str(resid2) + ":" + chain2
1649  output["SigmoidalCrossLinkMS_Score_" + crosslinker + "_" +
1650  label + "_" + self.label] = str(self.weight * ln.unprotected_evaluate(None))
1651  d0 = IMP.core.XYZR(p0)
1652  d1 = IMP.core.XYZR(p1)
1653  output["SigmoidalCrossLinkMS_Distance_" +
1654  label + "_" + self.label] = str(IMP.core.get_distance(d0, d1))
1655 
1656  return output
1657 
1658 @IMP.deprecated_object("2.5", "Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1659 class ISDCrossLinkMS(IMP.pmi.restraints._NuisancesBase):
1660  def __init__(self, representation,
1661  restraints_file,
1662  length,
1663  jackknifing=None,
1664  # jackknifing (Float [0,1]) default=None; percentage of data to be
1665  # removed randomly
1666  resolution=None,
1667  # resolution (Non-negative Integer) default=None; percentage of data
1668  # to be removed randomly
1669  slope=0.0,
1670  inner_slope=0.0,
1671  columnmapping=None,
1672  rename_dict=None,
1673  offset_dict=None,
1674  csvfile=False,
1675  ids_map=None,
1676  radius_map=None,
1677  filters=None,
1678  label="None",
1679  filelabel="None",
1680  automatic_sigma_classification=False,
1681  attributes_for_label=None):
1682 
1683  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2,idscore, XL unique id
1684  # by default column 0 = protein1; column 1 = protein2; column 2 = residue1; column 3 = residue2;
1685  # column 4 = idscores
1686  # attributes_for_label: anything in the csv database that must be added to the label
1687  # slope is the slope defined on the linear function
1688  # inner_slope is the slope defined on the restraint directly
1689  # suggestion: do not use both!
1690 
1691  if type(representation) != list:
1692  representations = [representation]
1693  else:
1694  representations = representation
1695 
1696  if columnmapping is None:
1697  columnmapping = {}
1698  columnmapping["Protein1"] = 0
1699  columnmapping["Protein2"] = 1
1700  columnmapping["Residue1"] = 2
1701  columnmapping["Residue2"] = 3
1702  columnmapping["IDScore"] = 4
1703  columnmapping["XLUniqueID"] = 5
1704 
1705  if csvfile:
1706  # in case the file is a csv file
1707  # columnmapping will contain the field names
1708  # that compare in the first line of the csv file
1709  db = IMP.pmi.tools.get_db_from_csv(restraints_file)
1710  else:
1711  db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1712 
1713  indb = open("included." + filelabel + ".xl.db", "w")
1714  exdb = open("excluded." + filelabel + ".xl.db", "w")
1715  midb = open("missing." + filelabel + ".xl.db", "w")
1716 
1717  self.m = representations[0].prot.get_model()
1718  self.rs = IMP.RestraintSet(self.m, 'data')
1719  self.rspsi = IMP.RestraintSet(self.m, 'prior_psi')
1720  self.rssig = IMP.RestraintSet(self.m, 'prior_sigmas')
1721  self.rslin = IMP.RestraintSet(self.m, 'prior_linear')
1722  self.rslen = IMP.RestraintSet(self.m, 'prior_length')
1723 
1724  self.weight = 1.0
1725  self.label = label
1726  self.pairs = []
1727  self.sigma_dictionary = {}
1728  self.psi_dictionary = {}
1729  self.psi_is_sampled = True
1730  self.sigma_is_sampled = True
1731 
1732  self.dataset = representations[0].get_file_dataset(restraints_file)
1733  if not self.dataset and os.path.exists(restraints_file):
1734  l = ihm.location.InputFileLocation(restraints_file,
1735  details="Crosslinks")
1736  self.dataset = ihm.dataset.CXMSDataset(l)
1737 
1738  xl_groups = [p.get_cross_link_group(self)
1739  for p, state in representations[0]._protocol_output]
1740 
1741  # isd_map is a dictionary/map that is used to determine the psi
1742  # parameter from identity scores (such as ID-Score, or FDR)
1743  if ids_map is None:
1744  self.ids_map = IMP.pmi.tools.map()
1745  self.ids_map.set_map_element(20.0, 0.05)
1746  self.ids_map.set_map_element(65.0, 0.01)
1747  else:
1748  self.ids_map = ids_map
1749 
1750  if radius_map is None:
1751  self.radius_map = IMP.pmi.tools.map()
1752  if automatic_sigma_classification:
1753  self.radius_map.set_map_element(10, 10)
1754  self.radius_map.set_map_element(1, 1)
1755  else:
1756  self.radius_map = radius_map
1757 
1758  self.outputlevel = "low"
1759 
1760  # small linear contribution for long range
1761  self.linear = IMP.core.Linear(0, 0.0)
1762  self.linear.set_slope(slope)
1763  dps2 = IMP.core.DistancePairScore(self.linear)
1764 
1765  protein1 = columnmapping["Protein1"]
1766  protein2 = columnmapping["Protein2"]
1767  residue1 = columnmapping["Residue1"]
1768  residue2 = columnmapping["Residue2"]
1769  idscore = columnmapping["IDScore"]
1770  try:
1771  xluniqueid = columnmapping["XLUniqueID"]
1772  except:
1773  xluniqueid = None
1774 
1775  restraints = []
1776 
1777  # we need this dictionary to create ambiguity (i.e., multistate)
1778  # if one id is already present in the dictionary, add the term to the
1779  # corresponding already generated restraint
1780 
1781  uniqueid_restraints_map = {}
1782 
1783  for nxl, entry in enumerate(db):
1784 
1785  if not jackknifing is None:
1786 
1787  # to be implemented
1788  # the problem is that in the replica exchange
1789  # you have to broadcast the same restraints to every
1790  # replica
1791 
1792  raise NotImplementedError("jackknifing not yet implemented")
1793 
1794  if not csvfile:
1795  tokens = entry.split()
1796  if len(tokens)==0:
1797  continue
1798 
1799  # skip character
1800  if (tokens[0] == "#"):
1801  continue
1802  try:
1803  r1 = int(tokens[residue1])
1804  c1 = tokens[protein1]
1805  r2 = int(tokens[residue2])
1806  c2 = tokens[protein2]
1807 
1808  if offset_dict is not None:
1809  if c1 in offset_dict: r1+=offset_dict[c1]
1810  if c2 in offset_dict: r2+=offset_dict[c2]
1811 
1812  if rename_dict is not None:
1813  if c1 in rename_dict: c1=rename_dict[c1]
1814  if c2 in rename_dict: c2=rename_dict[c2]
1815 
1816  if idscore is None:
1817  ids = 1.0
1818  else:
1819  ids = float(tokens[idscore])
1820  if xluniqueid is None:
1821  xlid = str(nxl)
1822  else:
1823  xlid = tokens[xluniqueid]
1824  except:
1825  print("this line was not accessible " + str(entry))
1826  if residue1 not in entry: print(str(residue1)+" keyword not in database")
1827  if residue2 not in entry: print(str(residue2)+" keyword not in database")
1828  if protein1 not in entry: print(str(protein1)+" keyword not in database")
1829  if protein2 not in entry: print(str(protein2)+" keyword not in database")
1830  if idscore not in entry: print(str(idscore)+" keyword not in database")
1831  if xluniqueid not in entry: print(str(xluniqueid)+" keyword not in database")
1832  continue
1833 
1834  else:
1835  if filters is not None:
1836  if eval(IMP.pmi.tools.cross_link_db_filter_parser(filters)) == False:
1837  exdb.write(str(entry) + "\n")
1838  continue
1839 
1840  try:
1841  r1 = int(entry[residue1])
1842  c1 = entry[protein1]
1843  r2 = int(entry[residue2])
1844  c2 = entry[protein2]
1845 
1846  if offset_dict is not None:
1847  if c1 in offset_dict: r1+=offset_dict[c1]
1848  if c2 in offset_dict: r2+=offset_dict[c2]
1849 
1850  if rename_dict is not None:
1851  if c1 in rename_dict: c1=rename_dict[c1]
1852  if c2 in rename_dict: c2=rename_dict[c2]
1853 
1854  if idscore is None:
1855  ids = 1.0
1856  else:
1857  try:
1858  ids = float(entry[idscore])
1859  except ValueError:
1860  ids = entry[idscore]
1861  if xluniqueid is None:
1862  xlid = str(nxl)
1863  else:
1864  xlid = entry[xluniqueid]
1865 
1866  except:
1867  print("this line was not accessible " + str(entry))
1868  if residue1 not in entry: print(str(residue1)+" keyword not in database")
1869  if residue2 not in entry: print(str(residue2)+" keyword not in database")
1870  if protein1 not in entry: print(str(protein1)+" keyword not in database")
1871  if protein2 not in entry: print(str(protein2)+" keyword not in database")
1872  if idscore not in entry: print(str(idscore)+" keyword not in database")
1873  if xluniqueid not in entry: print(str(xluniqueid)+" keyword not in database")
1874  continue
1875 
1876  # todo: check that offset is handled correctly
1877  ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
1878  group), group)
1879  for p, group in zip(representations[0]._protocol_output,
1880  xl_groups)]
1881 
1882  for nstate, r in enumerate(representations):
1883  # loop over every state
1884 
1885  ps1 = IMP.pmi.tools.select(
1886  r,
1887  resolution=resolution,
1888  name=c1,
1889  name_is_ambiguous=False,
1890  residue=r1)
1891  ps2 = IMP.pmi.tools.select(
1892  r,
1893  resolution=resolution,
1894  name=c2,
1895  name_is_ambiguous=False,
1896  residue=r2)
1897 
1898  if len(ps1) > 1:
1899  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1900  elif len(ps1) == 0:
1901  print("ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1902  midb.write(str(entry) + "\n")
1903  continue
1904 
1905  if len(ps2) > 1:
1906  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1907  elif len(ps2) == 0:
1908  print("ISDCrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1909  midb.write(str(entry) + "\n")
1910  continue
1911 
1912  p1 = ps1[0]
1913  p2 = ps2[0]
1914 
1915  if (p1 == p2) and (r1 == r2) :
1916  print("ISDCrossLinkMS Restraint: WARNING> on the identical bead particles and the identical residues, thus skipping this cross-link.")
1917  continue
1918 
1919  if xlid in uniqueid_restraints_map:
1920  print("Appending a crosslink restraint into the uniqueID %s" % str(xlid))
1921  dr = uniqueid_restraints_map[xlid]
1922  else:
1923  print("Generating a NEW crosslink restraint with a uniqueID %s" % str(xlid))
1925  self.m,
1926  length,
1927  inner_slope)
1928  restraints.append(dr)
1929  uniqueid_restraints_map[xlid] = dr
1930 
1931  mappedr1 = self.radius_map.get_map_element(
1932  IMP.pmi.Uncertainty(p1).get_uncertainty())
1933  sigma1 = self.get_sigma(mappedr1)[0]
1934  mappedr2 = self.radius_map.get_map_element(
1935  IMP.pmi.Uncertainty(p2).get_uncertainty())
1936  sigma2 = self.get_sigma(mappedr2)[0]
1937 
1938  psival = self.ids_map.get_map_element(ids)
1939  psi = self.get_psi(psival)[0]
1940 
1941 
1942  p1i = p1.get_particle_index()
1943  p2i = p2.get_particle_index()
1944  s1i = sigma1.get_particle().get_index()
1945  s2i = sigma2.get_particle().get_index()
1946 
1947  #print nstate, p1i, p2i, p1.get_name(), p2.get_name()
1948 
1949  psii = psi.get_particle().get_index()
1950 
1951  dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
1952  print("--------------")
1953  print("ISDCrossLinkMS: generating cross-link restraint between")
1954  print("ISDCrossLinkMS: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
1955  print("ISDCrossLinkMS: with sigma1 %f sigma2 %f psi %s" % (mappedr1, mappedr2, psival))
1956  print("ISDCrossLinkMS: between particles %s and %s" % (p1.get_name(), p2.get_name()))
1957  print("==========================================\n")
1958  indb.write(str(entry) + "\n")
1959  for p, ex_xl in zip(representations[0]._protocol_output,
1960  ex_xls):
1961  p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
1962  sigma1, sigma2, psi, ex_xl[1])
1963 
1964  # check if the two residues belong to the same rigid body
1967  IMP.core.RigidMember(p1).get_rigid_body() ==
1968  IMP.core.RigidMember(p2).get_rigid_body()):
1969  xlattribute = "intrarb"
1970  else:
1971  xlattribute = "interrb"
1972 
1973  if csvfile:
1974  if not attributes_for_label is None:
1975  for a in attributes_for_label:
1976  xlattribute = xlattribute + "_" + str(entry[a])
1977 
1978  xlattribute = xlattribute + "-State:" + str(nstate)
1979 
1980  dr.set_name(
1981  xlattribute + "-" + c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2) + "_" + self.label)
1982 
1983  if p1i != p2i:
1984  pr = IMP.core.PairRestraint(self.m, dps2, (p1i, p2i))
1985  pr.set_name(
1986  xlattribute + "-" + c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2) + "_" + self.label)
1987  self.rslin.add_restraint(pr)
1988 
1989 
1990  self.pairs.append(
1991  (p1,
1992  p2,
1993  dr,
1994  r1,
1995  c1,
1996  r2,
1997  c2,
1998  xlattribute,
1999  mappedr1,
2000  mappedr2,
2001  psival,
2002  xlid,
2003  nstate,
2004  ids))
2005 
2006  lw = IMP.isd.LogWrapper(restraints, self.weight)
2007  self.rs.add_restraint(lw)
2008 
2009  # Provide self.model for compatibility with newer code
2010  model = property(lambda s: s.m)
2011 
2012  def set_weight(self, weight):
2013  self.weight = weight
2014  self.rs.set_weight(weight)
2015 
2016  def set_slope_linear_term(self, slope):
2017  self.linear.set_slope(slope)
2018 
2019  def set_label(self, label):
2020  self.label = label
2021 
2022  def add_to_model(self):
2023  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
2024  IMP.pmi.tools.add_restraint_to_model(self.m, self.rspsi)
2025  IMP.pmi.tools.add_restraint_to_model(self.m, self.rssig)
2026  IMP.pmi.tools.add_restraint_to_model(self.m, self.rslen)
2027  IMP.pmi.tools.add_restraint_to_model(self.m, self.rslin)
2028 
2029  def get_hierarchies(self):
2030  return self.prot
2031 
2032  def get_restraint_sets(self):
2033  return self.rs
2034 
2035  def get_restraint(self):
2036  return self.rs
2037 
2038  def get_restraint_for_rmf(self):
2039  return self.rslin
2040 
2041  def get_restraints(self):
2042  rlist = []
2043  for r in self.rs.get_restraints():
2044  rlist.append(IMP.core.PairRestraint.get_from(r))
2045  return rlist
2046 
2047  def get_particle_pairs(self):
2048  ppairs = []
2049  for i in range(len(self.pairs)):
2050  p0 = self.pairs[i][0]
2051  p1 = self.pairs[i][1]
2052  ppairs.append((p0, p1))
2053  return ppairs
2054 
2055  def set_output_level(self, level="low"):
2056  # this might be "low" or "high"
2057  self.outputlevel = level
2058 
2059  def set_psi_is_sampled(self, is_sampled=True):
2060  self.psi_is_sampled = is_sampled
2061 
2062  def set_sigma_is_sampled(self, is_sampled=True):
2063  self.sigma_is_sampled = is_sampled
2064 
2065  def get_label(self,pairs_index):
2066  resid1 = self.pairs[pairs_index][3]
2067  chain1 = self.pairs[pairs_index][4]
2068  resid2 = self.pairs[pairs_index][5]
2069  chain2 = self.pairs[pairs_index][6]
2070  attribute = self.pairs[pairs_index][7]
2071  rad1 = self.pairs[pairs_index][8]
2072  rad2 = self.pairs[pairs_index][9]
2073  psi = self.pairs[pairs_index][10]
2074  xlid= self.pairs[pairs_index][11]
2075  label = attribute + "-" + \
2076  str(resid1) + ":" + chain1 + "_" + str(resid2) + ":" + \
2077  chain2 + "-" + str(rad1) + "-" + str(rad2) + "-" + str(psi)
2078  return label
2079 
2080  def write_db(self,filename):
2081  cldb=IMP.pmi.output.CrossLinkIdentifierDatabase()
2082 
2083  for pairs_index in range(len(self.pairs)):
2084 
2085  resid1 = self.pairs[pairs_index][3]
2086  chain1 = self.pairs[pairs_index][4]
2087  resid2 = self.pairs[pairs_index][5]
2088  chain2 = self.pairs[pairs_index][6]
2089  attribute = self.pairs[pairs_index][7]
2090  rad1 = self.pairs[pairs_index][8]
2091  rad2 = self.pairs[pairs_index][9]
2092  psi = self.pairs[pairs_index][10]
2093  xlid= self.pairs[pairs_index][11]
2094  nstate=self.pairs[pairs_index][12]
2095  ids=self.pairs[pairs_index][13]
2096 
2097  label=self.get_label(pairs_index)
2098  cldb.set_unique_id(label,xlid)
2099  cldb.set_protein1(label,chain1)
2100  cldb.set_protein2(label,chain2)
2101  cldb.set_residue1(label,resid1)
2102  cldb.set_residue2(label,resid2)
2103  cldb.set_idscore(label,ids)
2104  cldb.set_state(label,nstate)
2105  cldb.set_sigma1(label,rad1)
2106  cldb.set_sigma2(label,rad2)
2107  cldb.set_psi(label,psi)
2108  cldb.write(filename)
2109 
2110 
2111  def get_output(self):
2112  # content of the crosslink database pairs
2113  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
2114  output = {}
2115  score = self.weight * self.rs.unprotected_evaluate(None)
2116  output["_TotalScore"] = str(score)
2117  output["ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
2118  output["ISDCrossLinkMS_PriorSig_Score_" +
2119  self.label] = self.rssig.unprotected_evaluate(None)
2120  output["ISDCrossLinkMS_PriorPsi_Score_" +
2121  self.label] = self.rspsi.unprotected_evaluate(None)
2122  output["ISDCrossLinkMS_Linear_Score_" +
2123  self.label] = self.rslin.unprotected_evaluate(None)
2124  for i in range(len(self.pairs)):
2125 
2126  label=self.get_label(i)
2127  ln = self.pairs[i][2]
2128  p0 = self.pairs[i][0]
2129  p1 = self.pairs[i][1]
2130  output["ISDCrossLinkMS_Score_" +
2131  label + "_" + self.label] = str(self.weight * -log(ln.unprotected_evaluate(None)))
2132 
2133  d0 = IMP.core.XYZ(p0)
2134  d1 = IMP.core.XYZ(p1)
2135  output["ISDCrossLinkMS_Distance_" +
2136  label + "_" + self.label] = str(IMP.core.get_distance(d0, d1))
2137 
2138 
2139  for psiindex in self.psi_dictionary:
2140  output["ISDCrossLinkMS_Psi_" +
2141  str(psiindex) + "_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2142 
2143  for resolution in self.sigma_dictionary:
2144  output["ISDCrossLinkMS_Sigma_" +
2145  str(resolution) + "_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
2146 
2147 
2148  return output
2149 
2150  def get_particles_to_sample(self):
2151  ps = {}
2152  for resolution in self.sigma_dictionary:
2153  if self.sigma_dictionary[resolution][2] and self.sigma_is_sampled:
2154  ps["Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) + "_" + self.label] =\
2155  ([self.sigma_dictionary[resolution][0]],
2156  self.sigma_dictionary[resolution][1])
2157  if self.psi_is_sampled:
2158  for psiindex in self.psi_dictionary:
2159  if self.psi_dictionary[psiindex][2]:
2160  ps["Nuisances_ISDCrossLinkMS_Psi_" +
2161  str(psiindex) + "_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2162  return ps
2163 
2164 #
2165 class CysteineCrossLinkRestraint(object):
2166  def __init__(self, representations, filename, cbeta=False,
2167  betatuple=(0.03, 0.1),
2168  disttuple=(0.0, 25.0, 1000),
2169  omegatuple=(1.0, 1000.0, 50),
2170  sigmatuple=(0.3, 0.3, 1),
2171  betaissampled=True,
2172  sigmaissampled=False,
2173  weightissampled=True,
2174  epsilonissampled=True
2175  ):
2176  # the file must have residue1 chain1 residue2 chain2 fractionvalue epsilonname
2177  # epsilonname is a name for the epsilon particle that must be used for that particular
2178  # residue pair, eg, "Epsilon-Intra-Solvent", or
2179  # "Epsilon-Solvent-Membrane", etc.
2180 
2181  self.representations = representations
2182  self.m = self.representations[0].prot.get_model()
2183  self.rs = IMP.RestraintSet(self.m, 'Cysteine_Crosslink')
2184  self.cbeta = cbeta
2185  self.epsilonmaxtrans = 0.01
2186  self.sigmamaxtrans = 0.1
2187  self.betamaxtrans = 0.01
2188  self.weightmaxtrans = 0.1
2189  self.label = "None"
2190  self.outputlevel = "low"
2191  self.betaissampled = betaissampled
2192  self.sigmaissampled = sigmaissampled
2193  self.weightissampled = weightissampled
2194  self.epsilonissampled = epsilonissampled
2195 
2196  betalower = betatuple[0]
2197  betaupper = betatuple[1]
2198  beta = 0.04
2199  sigma = 10.0
2200  betangrid = 30
2201  crossdataprior = 1
2202 
2203  # beta
2204  self.beta = IMP.pmi.tools.SetupNuisance(
2205  self.m,
2206  beta,
2207  betalower,
2208  betaupper,
2209  betaissampled).get_particle(
2210  )
2211  # sigma
2212  self.sigma = IMP.pmi.tools.SetupNuisance(
2213  self.m,
2214  sigma,
2215  sigmatuple[0],
2216  sigmatuple[1],
2217  sigmaissampled).get_particle()
2218  # population particle
2219  self.weight = IMP.pmi.tools.SetupWeight(
2220  self.m,
2221  weightissampled).get_particle(
2222  )
2223 
2224  # read the file
2225  fl = IMP.pmi.tools.open_file_or_inline_text(filename)
2226  # determine the upperlimit for epsilon
2227  # and also how many epsilon are needed
2228  self.epsilons = {}
2229  data = []
2230  for line in fl:
2231  t = line.split()
2232  if t[0][0] == "#":
2233  continue
2234  if t[5] in self.epsilons:
2235  if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
2236  self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
2237  else:
2238  self.epsilons[t[5]] = IMP.pmi.tools.SetupNuisance(self.m,
2239  0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
2240  up = self.epsilons[t[5]].get_upper()
2241  low = self.epsilons[t[5]].get_lower()
2242  if up < low:
2243  self.epsilons[t[5]].set_upper(low)
2244 
2245  data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
2246 
2247  # create CrossLinkData
2248  if not self.cbeta:
2249  crossdata = IMP.pmi.tools.get_cross_link_data(
2250  "cysteine", "cysteine_CA_FES.txt.standard",
2251  disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2252  else:
2253  crossdata = IMP.pmi.tools.get_cross_link_data(
2254  "cysteine", "cysteine_CB_FES.txt.standard",
2255  disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2256 
2257  # create grids needed by CysteineCrossLinkData
2258  fmod_grid = IMP.pmi.tools.get_grid(0.0, 1.0, 300, True)
2259  omega2_grid = IMP.pmi.tools.get_log_grid(0.001, 10000.0, 100)
2260  beta_grid = IMP.pmi.tools.get_log_grid(betalower, betaupper, betangrid)
2261 
2262  for d in data:
2263  print("--------------")
2264  print("CysteineCrossLink: attempting to create a restraint " + str(d))
2265  resid1 = d[0]
2266  chain1 = d[1]
2267  resid2 = d[2]
2268  chain2 = d[3]
2269  fexp = d[4]
2270  epslabel = d[5]
2271 
2272  # CysteineCrossLinkData
2273 
2275  fexp,
2276  fmod_grid,
2277  omega2_grid,
2278  beta_grid)
2279 
2281  self.m,
2282  self.beta,
2283  self.sigma,
2284  self.epsilons[epslabel],
2285  self.weight,
2286  crossdata,
2287  ccldata)
2288 
2289  failed = False
2290  for i, representation in enumerate(self.representations):
2291 
2292  if not self.cbeta:
2293  p1 = None
2294  p2 = None
2295 
2296  p1 = IMP.pmi.tools.select(representation,
2297  resolution=1, name=chain1,
2298  name_is_ambiguous=False, residue=resid1)[0]
2299 
2300  if p1 is None:
2301  failed = True
2302 
2303  p2 = IMP.pmi.tools.select(representation,
2304  resolution=1, name=chain2,
2305  name_is_ambiguous=False, residue=resid2)[0]
2306 
2307  if p2 is None:
2308  failed = True
2309 
2310  else:
2311  # use cbetas
2312  p1 = []
2313  p2 = []
2314  for t in range(-1, 2):
2315  p = IMP.pmi.tools.select(representation,
2316  resolution=1, name=chain1,
2317  name_is_ambiguous=False, residue=resid1 + t)
2318 
2319  if len(p) == 1:
2320  p1 += p
2321  else:
2322  failed = True
2323  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
2324 
2325  p = IMP.pmi.tools.select(representation,
2326  resolution=1, name=chain2,
2327  name_is_ambiguous=False, residue=resid2 + t)
2328 
2329  if len(p) == 1:
2330  p2 += p
2331  else:
2332  failed = True
2333  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
2334 
2335  if not self.cbeta:
2336  if (p1 is not None and p2 is not None):
2337  ccl.add_contribution(p1, p2)
2338  d1 = IMP.core.XYZ(p1)
2339  d2 = IMP.core.XYZ(p2)
2340 
2341  print("Distance_" + str(resid1) + "_" + chain1 + ":" + str(resid2) + "_" + chain2, IMP.core.get_distance(d1, d2))
2342 
2343  else:
2344  if (len(p1) == 3 and len(p2) == 3):
2345  p11n = p1[0].get_name()
2346  p12n = p1[1].get_name()
2347  p13n = p1[2].get_name()
2348  p21n = p2[0].get_name()
2349  p22n = p2[1].get_name()
2350  p23n = p2[2].get_name()
2351 
2352  print("CysteineCrossLink: generating CB cysteine cross-link restraint between")
2353  print("CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
2354  print("CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
2355 
2356  ccl.add_contribution(p1, p2)
2357 
2358  if not failed:
2359  self.rs.add_restraint(ccl)
2360  ccl.set_name("CysteineCrossLink_" + str(resid1)
2361  + "_" + chain1 + ":" + str(resid2) + "_" + chain2)
2362 
2363  def set_label(self, label):
2364  self.label = label
2365 
2366  def add_to_model(self):
2367  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
2368 
2369  def get_particles_to_sample(self):
2370  ps = {}
2371  if self.epsilonissampled:
2372  for eps in self.epsilons.keys():
2373  ps["Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
2374  self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
2375  if self.betaissampled:
2376  ps["Nuisances_CysteineCrossLinkRestraint_beta_" +
2377  self.label] = ([self.beta], self.betamaxtrans)
2378  if self.weightissampled:
2379  ps["Weights_CysteineCrossLinkRestraint_" +
2380  self.label] = ([self.weight], self.weightmaxtrans)
2381  if self.sigmaissampled:
2382  ps["Nuisances_CysteineCrossLinkRestraint_" +
2383  self.label] = ([self.sigma], self.sigmamaxtrans)
2384  return ps
2385 
2386  def set_output_level(self, level="low"):
2387  # this might be "low" or "high"
2388  self.outputlevel = level
2389 
2390  def get_restraint(self):
2391  return self.rs
2392 
2393  def get_sigma(self):
2394  return self.sigma
2395 
2396  def get_output(self):
2397  output = {}
2398  score = self.rs.unprotected_evaluate(None)
2399  output["_TotalScore"] = str(score)
2400  output["CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
2401  output["CysteineCrossLinkRestraint_sigma_" +
2402  self.label] = str(self.sigma.get_scale())
2403  for eps in self.epsilons.keys():
2404  output["CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
2405  self.label] = str(self.epsilons[eps].get_scale())
2406  output["CysteineCrossLinkRestraint_beta_" +
2407  self.label] = str(self.beta.get_scale())
2408  for n in range(self.weight.get_number_of_states()):
2409  output["CysteineCrossLinkRestraint_weights_" +
2410  str(n) + "_" + self.label] = str(self.weight.get_weight(n))
2411 
2412  if self.outputlevel == "high":
2413  for rst in self.rs.get_restraints():
2414  output["CysteineCrossLinkRestraint_Total_Frequency_" +
2415  IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2416  "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
2417  output["CysteineCrossLinkRestraint_Standard_Error_" +
2418  IMP.isd.CysteineCrossLinkRestraint.get_from(
2419  rst).get_name(
2420  ) + "_"
2421  + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
2422  if len(self.representations) > 1:
2423  for i in range(len(self.prots)):
2424  output["CysteineCrossLinkRestraint_Frequency_Contribution_" +
2425  IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2426  "_State_" + str(i) + "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_frequencies()[i]
2427 
2428  return output
def set_sigma_is_sampled
Switch on/off the sampling of sigma particles.
Setup cross-link distance restraints from mass spectrometry data.
A function that is harmonic over an interval.
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.
Upper bound harmonic function (non-zero when feature > mean)
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
Miscellaneous utilities.
Definition: tools.py:1
Restrain atom pairs based on a set of crosslinks.
Add scale parameter to particle.
Definition: Scale.h:24
The type of an atom.
Add uncertainty to a particle.
Definition: Uncertainty.h:24
A score on the distance between the surfaces of two spheres.
log
Definition: log.py:1
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:690
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.
def get_restraint_for_rmf
Get the restraint for visualization in an RMF file.
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.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
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:69
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.
Simple sigmoidal score calculated between sphere surfaces.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
This restraint allows ambiguous crosslinking between multiple copies it is a variant of the Simplifie...
int get_state_index(Hierarchy h)
Walk up the hierarchy to find the current state.
Linear function
Definition: Linear.h:19
def plot_violations
Create CMM files, one for each state, of all crosslinks.
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
def deprecated_object
Python decorator to mark a class as deprecated.
Definition: __init__.py:9687
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...
VectorD< 3 > Vector3D
Definition: VectorD.h:395
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
Class to handle individual particles of a Model object.
Definition: Particle.h:41
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
Definition: VectorD.h:209
Restraint * create_connectivity_restraint(const Selections &s, double x0, double k, std::string name="Connectivity%1%")
Create a restraint connecting the selections.
Applies a PairScore to a Pair.
Definition: PairRestraint.h:29
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:726
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
def cross_link_db_filter_parser
example '"{ID_Score}" > 28 AND "{Sample}" == "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample...
Definition: tools.py:501
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
def get_output
Get outputs to write to stat files.
Base class for PMI restraints, which wrap IMP.Restraint(s).
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A restraint for cysteine cross-linking data.
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27