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