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