IMP logo
IMP Reference Guide  2.20.2
The Integrative Modeling Platform
/restraints/crosslinking.py
1 """@namespace IMP.pmi1.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.pmi1.tools
13 import IMP.pmi1.output
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.pmi1.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.pmi1.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.pmi1.io.crosslink.CrossLinkDataBase):
80  raise TypeError("CrossLinkingMassSpectrometryRestraint: CrossLinkDataBase should be an IMP.pmi1.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  if representations:
110  xl_groups = [p.get_cross_link_group(self)
111  for p, state in representations[0]._protocol_output]
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.pmi1.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.pmi1.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  if representations:
156  ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
157  group), group)
158  for p, group in
159  zip(representations[0]._protocol_output,
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.pmi1.tools.select(
197  r,
198  resolution=resolution,
199  name=c1,
200  name_is_ambiguous=False,
201  residue=r1)
202  ps2 = IMP.pmi1.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  if representations:
284  for p, ex_xl in zip(representations[0]._protocol_output,
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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.Uncertainty(p).get_uncertainty()
1109 
1110  for p in ps2:
1111  rad2 += IMP.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.tools.select(
1273  representation,
1274  resolution=resolution,
1275  name=c1,
1276  name_is_ambiguous=False,
1277  residue=r1)
1278  ps2 = IMP.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.tools.select(
1514  representation,
1515  resolution=resolution,
1516  name=c1,
1517  name_is_ambiguous=False,
1518  residue=r1)
1519  ps2 = IMP.pmi1.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.pmi1.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.pmi1.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.pmi1.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1659 class ISDCrossLinkMS(IMP.pmi1.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.pmi1.tools.get_db_from_csv(restraints_file)
1710  else:
1711  db = IMP.pmi1.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.pmi1.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.pmi1.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.pmi1.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.pmi1.tools.select(
1886  r,
1887  resolution=resolution,
1888  name=c1,
1889  name_is_ambiguous=False,
1890  residue=r1)
1891  ps2 = IMP.pmi1.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.pmi1.Uncertainty(p1).get_uncertainty())
1933  sigma1 = self.get_sigma(mappedr1)[0]
1934  mappedr2 = self.radius_map.get_map_element(
1935  IMP.pmi1.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  def set_weight(self, weight):
2010  self.weight = weight
2011  self.rs.set_weight(weight)
2012 
2013  def set_slope_linear_term(self, slope):
2014  self.linear.set_slope(slope)
2015 
2016  def set_label(self, label):
2017  self.label = label
2018 
2019  def add_to_model(self):
2020  IMP.pmi1.tools.add_restraint_to_model(self.m, self.rs)
2021  IMP.pmi1.tools.add_restraint_to_model(self.m, self.rspsi)
2022  IMP.pmi1.tools.add_restraint_to_model(self.m, self.rssig)
2023  IMP.pmi1.tools.add_restraint_to_model(self.m, self.rslen)
2024  IMP.pmi1.tools.add_restraint_to_model(self.m, self.rslin)
2025 
2026  def get_hierarchies(self):
2027  return self.prot
2028 
2029  def get_restraint_sets(self):
2030  return self.rs
2031 
2032  def get_restraint(self):
2033  return self.rs
2034 
2035  def get_restraint_for_rmf(self):
2036  return self.rslin
2037 
2038  def get_restraints(self):
2039  rlist = []
2040  for r in self.rs.get_restraints():
2041  rlist.append(IMP.core.PairRestraint.get_from(r))
2042  return rlist
2043 
2044  def get_particle_pairs(self):
2045  ppairs = []
2046  for i in range(len(self.pairs)):
2047  p0 = self.pairs[i][0]
2048  p1 = self.pairs[i][1]
2049  ppairs.append((p0, p1))
2050  return ppairs
2051 
2052  def set_output_level(self, level="low"):
2053  # this might be "low" or "high"
2054  self.outputlevel = level
2055 
2056  def set_psi_is_sampled(self, is_sampled=True):
2057  self.psi_is_sampled = is_sampled
2058 
2059  def set_sigma_is_sampled(self, is_sampled=True):
2060  self.sigma_is_sampled = is_sampled
2061 
2062  def get_label(self,pairs_index):
2063  resid1 = self.pairs[pairs_index][3]
2064  chain1 = self.pairs[pairs_index][4]
2065  resid2 = self.pairs[pairs_index][5]
2066  chain2 = self.pairs[pairs_index][6]
2067  attribute = self.pairs[pairs_index][7]
2068  rad1 = self.pairs[pairs_index][8]
2069  rad2 = self.pairs[pairs_index][9]
2070  psi = self.pairs[pairs_index][10]
2071  xlid= self.pairs[pairs_index][11]
2072  label = attribute + "-" + \
2073  str(resid1) + ":" + chain1 + "_" + str(resid2) + ":" + \
2074  chain2 + "-" + str(rad1) + "-" + str(rad2) + "-" + str(psi)
2075  return label
2076 
2077  def write_db(self,filename):
2078  cldb=IMP.pmi1.output.CrossLinkIdentifierDatabase()
2079 
2080  for pairs_index in range(len(self.pairs)):
2081 
2082  resid1 = self.pairs[pairs_index][3]
2083  chain1 = self.pairs[pairs_index][4]
2084  resid2 = self.pairs[pairs_index][5]
2085  chain2 = self.pairs[pairs_index][6]
2086  attribute = self.pairs[pairs_index][7]
2087  rad1 = self.pairs[pairs_index][8]
2088  rad2 = self.pairs[pairs_index][9]
2089  psi = self.pairs[pairs_index][10]
2090  xlid= self.pairs[pairs_index][11]
2091  nstate=self.pairs[pairs_index][12]
2092  ids=self.pairs[pairs_index][13]
2093 
2094  label=self.get_label(pairs_index)
2095  cldb.set_unique_id(label,xlid)
2096  cldb.set_protein1(label,chain1)
2097  cldb.set_protein2(label,chain2)
2098  cldb.set_residue1(label,resid1)
2099  cldb.set_residue2(label,resid2)
2100  cldb.set_idscore(label,ids)
2101  cldb.set_state(label,nstate)
2102  cldb.set_sigma1(label,rad1)
2103  cldb.set_sigma2(label,rad2)
2104  cldb.set_psi(label,psi)
2105  cldb.write(filename)
2106 
2107 
2108  def get_output(self):
2109  # content of the crosslink database pairs
2110  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
2111  output = {}
2112  score = self.weight * self.rs.unprotected_evaluate(None)
2113  output["_TotalScore"] = str(score)
2114  output["ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
2115  output["ISDCrossLinkMS_PriorSig_Score_" +
2116  self.label] = self.rssig.unprotected_evaluate(None)
2117  output["ISDCrossLinkMS_PriorPsi_Score_" +
2118  self.label] = self.rspsi.unprotected_evaluate(None)
2119  output["ISDCrossLinkMS_Linear_Score_" +
2120  self.label] = self.rslin.unprotected_evaluate(None)
2121  for i in range(len(self.pairs)):
2122 
2123  label=self.get_label(i)
2124  ln = self.pairs[i][2]
2125  p0 = self.pairs[i][0]
2126  p1 = self.pairs[i][1]
2127  output["ISDCrossLinkMS_Score_" +
2128  label + "_" + self.label] = str(self.weight * -log(ln.unprotected_evaluate(None)))
2129 
2130  d0 = IMP.core.XYZ(p0)
2131  d1 = IMP.core.XYZ(p1)
2132  output["ISDCrossLinkMS_Distance_" +
2133  label + "_" + self.label] = str(IMP.core.get_distance(d0, d1))
2134 
2135 
2136  for psiindex in self.psi_dictionary:
2137  output["ISDCrossLinkMS_Psi_" +
2138  str(psiindex) + "_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2139 
2140  for resolution in self.sigma_dictionary:
2141  output["ISDCrossLinkMS_Sigma_" +
2142  str(resolution) + "_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
2143 
2144 
2145  return output
2146 
2147  def get_particles_to_sample(self):
2148  ps = {}
2149  for resolution in self.sigma_dictionary:
2150  if self.sigma_dictionary[resolution][2] and self.sigma_is_sampled:
2151  ps["Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) + "_" + self.label] =\
2152  ([self.sigma_dictionary[resolution][0]],
2153  self.sigma_dictionary[resolution][1])
2154  if self.psi_is_sampled:
2155  for psiindex in self.psi_dictionary:
2156  if self.psi_dictionary[psiindex][2]:
2157  ps["Nuisances_ISDCrossLinkMS_Psi_" +
2158  str(psiindex) + "_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2159  return ps
2160 
2161 #
2162 class CysteineCrossLinkRestraint(object):
2163  def __init__(self, representations, filename, cbeta=False,
2164  betatuple=(0.03, 0.1),
2165  disttuple=(0.0, 25.0, 1000),
2166  omegatuple=(1.0, 1000.0, 50),
2167  sigmatuple=(0.3, 0.3, 1),
2168  betaissampled=True,
2169  sigmaissampled=False,
2170  weightissampled=True,
2171  epsilonissampled=True
2172  ):
2173  # the file must have residue1 chain1 residue2 chain2 fractionvalue epsilonname
2174  # epsilonname is a name for the epsilon particle that must be used for that particular
2175  # residue pair, eg, "Epsilon-Intra-Solvent", or
2176  # "Epsilon-Solvent-Membrane", etc.
2177 
2178  self.representations = representations
2179  self.m = self.representations[0].prot.get_model()
2180  self.rs = IMP.RestraintSet(self.m, 'Cysteine_Crosslink')
2181  self.cbeta = cbeta
2182  self.epsilonmaxtrans = 0.01
2183  self.sigmamaxtrans = 0.1
2184  self.betamaxtrans = 0.01
2185  self.weightmaxtrans = 0.1
2186  self.label = "None"
2187  self.outputlevel = "low"
2188  self.betaissampled = betaissampled
2189  self.sigmaissampled = sigmaissampled
2190  self.weightissampled = weightissampled
2191  self.epsilonissampled = epsilonissampled
2192 
2193  betalower = betatuple[0]
2194  betaupper = betatuple[1]
2195  beta = 0.04
2196  sigma = 10.0
2197  betangrid = 30
2198  crossdataprior = 1
2199 
2200  # beta
2201  self.beta = IMP.pmi1.tools.SetupNuisance(
2202  self.m,
2203  beta,
2204  betalower,
2205  betaupper,
2206  betaissampled).get_particle(
2207  )
2208  # sigma
2209  self.sigma = IMP.pmi1.tools.SetupNuisance(
2210  self.m,
2211  sigma,
2212  sigmatuple[0],
2213  sigmatuple[1],
2214  sigmaissampled).get_particle()
2215  # population particle
2216  self.weight = IMP.pmi1.tools.SetupWeight(
2217  self.m,
2218  isoptimized=False).get_particle(
2219  )
2220 
2221  # read the file
2222  fl = IMP.pmi1.tools.open_file_or_inline_text(filename)
2223  # determine the upperlimit for epsilon
2224  # and also how many epsilon are needed
2225  self.epsilons = {}
2226  data = []
2227  for line in fl:
2228  t = line.split()
2229  if t[0][0] == "#":
2230  continue
2231  if t[5] in self.epsilons:
2232  if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
2233  self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
2234  else:
2235  self.epsilons[t[5]] = IMP.pmi1.tools.SetupNuisance(self.m,
2236  0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
2237  up = self.epsilons[t[5]].get_upper()
2238  low = self.epsilons[t[5]].get_lower()
2239  if up < low:
2240  self.epsilons[t[5]].set_upper(low)
2241 
2242  data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
2243 
2244  # create CrossLinkData
2245  if not self.cbeta:
2246  crossdata = IMP.pmi1.tools.get_cross_link_data(
2247  "cysteine", "cysteine_CA_FES.txt.standard",
2248  disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2249  else:
2250  crossdata = IMP.pmi1.tools.get_cross_link_data(
2251  "cysteine", "cysteine_CB_FES.txt.standard",
2252  disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
2253 
2254  # create grids needed by CysteineCrossLinkData
2255  fmod_grid = IMP.pmi1.tools.get_grid(0.0, 1.0, 300, True)
2256  omega2_grid = IMP.pmi1.tools.get_log_grid(0.001, 10000.0, 100)
2257  beta_grid = IMP.pmi1.tools.get_log_grid(betalower, betaupper, betangrid)
2258 
2259  for d in data:
2260  print("--------------")
2261  print("CysteineCrossLink: attempting to create a restraint " + str(d))
2262  resid1 = d[0]
2263  chain1 = d[1]
2264  resid2 = d[2]
2265  chain2 = d[3]
2266  fexp = d[4]
2267  epslabel = d[5]
2268 
2269  # CysteineCrossLinkData
2270 
2272  fexp,
2273  fmod_grid,
2274  omega2_grid,
2275  beta_grid)
2276 
2278  self.m,
2279  self.beta,
2280  self.sigma,
2281  self.epsilons[epslabel],
2282  self.weight,
2283  crossdata,
2284  ccldata)
2285 
2286  failed = False
2287  for i, representation in enumerate(self.representations):
2288 
2289  if not self.cbeta:
2290  p1 = None
2291  p2 = None
2292 
2293  p1 = IMP.pmi1.tools.select(representation,
2294  resolution=1, name=chain1,
2295  name_is_ambiguous=False, residue=resid1)[0]
2296 
2297  if p1 is None:
2298  failed = True
2299 
2300  p2 = IMP.pmi1.tools.select(representation,
2301  resolution=1, name=chain2,
2302  name_is_ambiguous=False, residue=resid2)[0]
2303 
2304  if p2 is None:
2305  failed = True
2306 
2307  else:
2308  # use cbetas
2309  p1 = []
2310  p2 = []
2311  for t in range(-1, 2):
2312  p = IMP.pmi1.tools.select(representation,
2313  resolution=1, name=chain1,
2314  name_is_ambiguous=False, residue=resid1 + t)
2315 
2316  if len(p) == 1:
2317  p1 += p
2318  else:
2319  failed = True
2320  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
2321 
2322  p = IMP.pmi1.tools.select(representation,
2323  resolution=1, name=chain2,
2324  name_is_ambiguous=False, residue=resid2 + t)
2325 
2326  if len(p) == 1:
2327  p2 += p
2328  else:
2329  failed = True
2330  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
2331 
2332  if not self.cbeta:
2333  if (p1 is not None and p2 is not None):
2334  ccl.add_contribution(p1, p2)
2335  d1 = IMP.core.XYZ(p1)
2336  d2 = IMP.core.XYZ(p2)
2337 
2338  print("Distance_" + str(resid1) + "_" + chain1 + ":" + str(resid2) + "_" + chain2, IMP.core.get_distance(d1, d2))
2339 
2340  else:
2341  if (len(p1) == 3 and len(p2) == 3):
2342  p11n = p1[0].get_name()
2343  p12n = p1[1].get_name()
2344  p13n = p1[2].get_name()
2345  p21n = p2[0].get_name()
2346  p22n = p2[1].get_name()
2347  p23n = p2[2].get_name()
2348 
2349  print("CysteineCrossLink: generating CB cysteine cross-link restraint between")
2350  print("CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
2351  print("CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
2352 
2353  ccl.add_contribution(p1, p2)
2354 
2355  if not failed:
2356  self.rs.add_restraint(ccl)
2357  ccl.set_name("CysteineCrossLink_" + str(resid1)
2358  + "_" + chain1 + ":" + str(resid2) + "_" + chain2)
2359 
2361  self.weight.get_particle()
2362  ).set_weights_are_optimized(weightissampled)
2363 
2364  def set_label(self, label):
2365  self.label = label
2366 
2367  def add_to_model(self):
2368  IMP.pmi1.tools.add_restraint_to_model(self.m, self.rs)
2369 
2370  def get_particles_to_sample(self):
2371  ps = {}
2372  if self.epsilonissampled:
2373  for eps in self.epsilons.keys():
2374  ps["Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
2375  self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
2376  if self.betaissampled:
2377  ps["Nuisances_CysteineCrossLinkRestraint_beta_" +
2378  self.label] = ([self.beta], self.betamaxtrans)
2379  if self.weightissampled:
2380  ps["Weights_CysteineCrossLinkRestraint_" +
2381  self.label] = ([self.weight], self.weightmaxtrans)
2382  if self.sigmaissampled:
2383  ps["Nuisances_CysteineCrossLinkRestraint_" +
2384  self.label] = ([self.sigma], self.sigmamaxtrans)
2385  return ps
2386 
2387  def set_output_level(self, level="low"):
2388  # this might be "low" or "high"
2389  self.outputlevel = level
2390 
2391  def get_restraint(self):
2392  return self.rs
2393 
2394  def get_sigma(self):
2395  return self.sigma
2396 
2397  def get_output(self):
2398  output = {}
2399  score = self.rs.unprotected_evaluate(None)
2400  output["_TotalScore"] = str(score)
2401  output["CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
2402  output["CysteineCrossLinkRestraint_sigma_" +
2403  self.label] = str(self.sigma.get_scale())
2404  for eps in self.epsilons.keys():
2405  output["CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
2406  self.label] = str(self.epsilons[eps].get_scale())
2407  output["CysteineCrossLinkRestraint_beta_" +
2408  self.label] = str(self.beta.get_scale())
2409  for n in range(self.weight.get_number_of_weights()):
2410  output["CysteineCrossLinkRestraint_weights_" +
2411  str(n) + "_" + self.label] = str(self.weight.get_weight(n))
2412 
2413  if self.outputlevel == "high":
2414  for rst in self.rs.get_restraints():
2415  output["CysteineCrossLinkRestraint_Total_Frequency_" +
2416  IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2417  "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
2418  output["CysteineCrossLinkRestraint_Standard_Error_" +
2419  IMP.isd.CysteineCrossLinkRestraint.get_from(
2420  rst).get_name(
2421  ) + "_"
2422  + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
2423  if len(self.representations) > 1:
2424  for i in range(len(self.prots)):
2425  output["CysteineCrossLinkRestraint_Frequency_Contribution_" +
2426  IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
2427  "_State_" + str(i) + "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_frequencies()[i]
2428 
2429  return output
Setup cross-link distance restraints from mass spectrometry data.
Add weights to a particle.
Definition: Weight.h:29
A function that is harmonic over an bounded interval.
A restraint for ambiguous cross-linking MS data and multiple state approach.
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
Restrain atom pairs based on a set of crosslinks.
def get_best_stats
return the probability, best distance, two coords, and possibly the psi for each xl ...
Classes to handle different kinds of restraints.
Add scale parameter to particle.
Definition: Scale.h:24
The type of an atom.
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:749
Score distance between two particle centers using a harmonic function.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
def get_restraint_for_rmf
get the dummy restraints to be displayed in the rmf file
Uniform distribution with harmonic boundaries.
Definition: UniformPrior.h:22
Object used to hold a set of restraints.
Definition: RestraintSet.h:41
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:86
Setup cross-link distance restraints at atomic level The "atomic" aspect is that it models the part...
def get_output
Get outputs to write to stat files.
Add uncertainty to a particle.
Definition: /Uncertainty.h:24
Miscellaneous utilities.
Definition: /tools.py:1
def load_nuisances_from_stat_file
Read a stat file and load all the sigmas.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Simple sigmoidal score calculated between sphere surfaces.
A decorator for a particle representing an atom.
Definition: atom/Atom.h:238
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi1.sampler object.
def set_sigma_is_sampled
Switch on/off the sampling of sigma particles.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
int get_state_index(Hierarchy h)
Walk up the hierarchy to find the current state.
Linear function
Definition: Linear.h:22
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi1.sampler object.
Base class for PMI restraints, which wrap IMP.Restraint(s).
Modify a set of continuous variables using a normal distribution.
Definition: NormalMover.h:23
Score a pair of particles based on the distance between their centers.
def deprecated_object
Python decorator to mark a class as deprecated.
Definition: __init__.py:9522
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: /tools.py:727
def get_movers
Get all need data to construct a mover in IMP.pmi1.dof class.
def cross_link_db_filter_parser
example '"{ID_Score}" > 28 AND "{Sample}" == "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample...
Definition: /tools.py:504
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 plot_violations
Create CMM files, one for each state, of all crosslinks.
def get_restraint_for_rmf
Get the restraint for visualization in an RMF file.
The general base class for IMP exceptions.
Definition: exception.h:48
VectorD< 3 > Vector3D
Definition: VectorD.h:425
Residue get_residue(Atom d, bool nothrow=false)
Return the Residue containing this atom.
Calculate the -Log of a list of restraints.
Definition: LogWrapper.h:26
This restraint allows ambiguous crosslinking between multiple copies it is a variant of the Simplifie...
def get_particle_pairs
Get a list of tuples containing the particle pairs.
Class to handle individual particles of a Model object.
Definition: Particle.h:43
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
Definition: VectorD.h:213
def get_output
Get the output of the restraint to be used by the IMP.pmi1.output object.
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:31
Classes for writing output files and processing them.
Definition: /output.py:1
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: /tools.py:53
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.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
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...
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