IMP logo
IMP Reference Guide  2.11.0
The Integrative Modeling Platform
restraints/crosslinking.py
1 """@namespace IMP.pmi.restraints.crosslinking
2 Restraints for handling crosslinking data.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 import IMP.algebra
9 import IMP.atom
10 import IMP.isd
11 import IMP.container
12 import IMP.pmi.tools
13 import IMP.pmi.output
15 import IMP.pmi.restraints
16 from math import log
17 from collections import defaultdict
18 import itertools
19 import operator
20 import os
21 import ihm.location
22 import ihm.dataset
23 import warnings
24 
26  """Setup cross-link distance restraints from mass spectrometry data.
27  The noise in the data and the structural uncertainty of cross-linked amino-acids
28  is inferred using Bayes theory of probability
29  @note Wraps an IMP::isd::CrossLinkMSRestraint
30  """
31  def __init__(self, representation=None,
32  root_hier=None,
33  CrossLinkDataBase=None,
34  length=10.0,
35  resolution=None,
36  slope=0.02,
37  label=None,
38  filelabel="None",
39  attributes_for_label=None,
40  weight=1.):
41  """Constructor.
42  @param representation DEPRECATED The IMP.pmi.representation.Representation
43  object that contain the molecular system
44  @param root_hier The canonical hierarchy containing all the states
45  @param CrossLinkDataBase The IMP.pmi.io.crosslink.CrossLinkDataBase
46  object that contains the cross-link dataset
47  @param length maximal cross-linker length (including the residue sidechains)
48  @param resolution what representation resolution should the cross-link
49  restraint be applied to.
50  @param slope The slope of a distance-linear scoring function that
51  funnels the score when the particles are
52  too far away. Suggested value 0.02.
53  @param label the extra text to label the restraint so that it is
54  searchable in the output
55  @param filelabel automatically generated file containing missing/included/excluded
56  cross-links will be labeled using this text
57  @param attributes_for_label
58  @param weight Weight of restraint
59  """
60 
61  use_pmi2 = True
62  if representation is not None:
63  use_pmi2 = False
64  if type(representation) != list:
65  representations = [representation]
66  else:
67  representations = representation
68  model = representations[0].prot.get_model()
69  elif root_hier is not None:
70  representations = []
71  model = root_hier.get_model()
72  else:
73  raise Exception("You must pass either representation or root_hier")
74 
75  super(CrossLinkingMassSpectrometryRestraint, self).__init__(
76  model, weight=weight, label=label)
77 
78  if CrossLinkDataBase is None:
79  raise Exception("You must pass a CrossLinkDataBase")
80  if not isinstance(CrossLinkDataBase,IMP.pmi.io.crosslink.CrossLinkDataBase):
81  raise TypeError("CrossLinkingMassSpectrometryRestraint: CrossLinkDataBase should be an IMP.pmi.io.crosslink.CrossLinkDataBase object")
82  self.CrossLinkDataBase = CrossLinkDataBase
83 
84  if resolution==0 or resolution is None:
85  raise Exception("You must pass a resolution and it can't be zero")
86 
87  indb = open("included." + filelabel + ".xl.db", "w")
88  exdb = open("excluded." + filelabel + ".xl.db", "w")
89  midb = open("missing." + filelabel + ".xl.db", "w")
90 
91  self.rs.set_name(self.rs.get_name() + "_Data")
92  self.rspsi = self._create_restraint_set("PriorPsi")
93  self.rssig = self._create_restraint_set("PriorSig")
94  self.rslin = self._create_restraint_set("Linear")
95 
96  # dummy linear restraint used for Chimera display
97  self.linear = IMP.core.Linear(0, 0.0)
98  self.linear.set_slope(0.0)
99  dps2 = IMP.core.DistancePairScore(self.linear)
100 
101  self.psi_is_sampled = True
102  self.sigma_is_sampled = True
103  self.psi_dictionary={}
104  self.sigma_dictionary={}
105  self.xl_list=[]
106  self.outputlevel = "low"
107 
108  restraints = []
109 
110  xl_groups = [p.get_cross_link_group(self)
111  for p, state in IMP.pmi.tools._all_protocol_outputs(
112  representations, root_hier)]
113 
114  # if PMI2, first add all the molecule copies as clones to the database
115  if use_pmi2:
116  copies_to_add = defaultdict(int)
117  print('gathering copies')
118  for xlid in self.CrossLinkDataBase.xlid_iterator():
119  for xl in self.CrossLinkDataBase[xlid]:
120  r1 = xl[self.CrossLinkDataBase.residue1_key]
121  c1 = xl[self.CrossLinkDataBase.protein1_key]
122  r2 = xl[self.CrossLinkDataBase.residue2_key]
123  c2 = xl[self.CrossLinkDataBase.protein2_key]
124  for c,r in ((c1,r1),(c2,r2)):
125  if c in copies_to_add:
126  continue
127  sel = IMP.atom.Selection(root_hier,
128  state_index=0,
129  molecule=c,
130  residue_index=r,
131  resolution=resolution).get_selected_particles()
132  if len(sel)>0:
133  copies_to_add[c] = len(sel)-1
134  print(copies_to_add)
135  for molname in copies_to_add:
136  if copies_to_add[molname]==0:
137  continue
138  fo1 = IMP.pmi.io.crosslink.FilterOperator(self.CrossLinkDataBase.protein1_key,operator.eq,molname)
139  self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein1_key,molname+'.0',fo1)
140  fo2 = IMP.pmi.io.crosslink.FilterOperator(self.CrossLinkDataBase.protein2_key,operator.eq,molname)
141  self.CrossLinkDataBase.set_value(self.CrossLinkDataBase.protein2_key,molname+'.0',fo2)
142  for ncopy in range(copies_to_add[molname]):
143  self.CrossLinkDataBase.clone_protein('%s.0'%molname,'%s.%i'%(molname,ncopy+1))
144  print('done pmi2 prelims')
145 
146  for xlid in self.CrossLinkDataBase.xlid_iterator():
147  new_contribution=True
148  for xl in self.CrossLinkDataBase[xlid]:
149 
150  r1 = xl[self.CrossLinkDataBase.residue1_key]
151  c1 = xl[self.CrossLinkDataBase.protein1_key]
152  r2 = xl[self.CrossLinkDataBase.residue2_key]
153  c2 = xl[self.CrossLinkDataBase.protein2_key]
154 
155  # todo: check that offset is handled correctly
156  ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
157  group), group)
158  for p, group in
159  zip(IMP.pmi.tools._all_protocol_outputs(
160  representations, root_hier),
161  xl_groups)]
162 
163  if use_pmi2:
164  iterlist = range(len(IMP.atom.get_by_type(root_hier,IMP.atom.STATE_TYPE)))
165  else:
166  iterlist = representations
167  for nstate, r in enumerate(iterlist):
168  # loop over every state
169  xl[self.CrossLinkDataBase.state_key]=nstate
170  xl[self.CrossLinkDataBase.data_set_name_key]=self.label
171 
172  if use_pmi2:
173  name1 = c1
174  name2 = c2
175  copy1 = 0
176  copy2 = 0
177  if '.' in c1:
178  name1,copy1 = c1.split('.')
179  if '.' in c2:
180  name2,copy2 = c2.split('.')
181  ps1 = IMP.atom.Selection(root_hier,
182  state_index=nstate,
183  molecule=name1,
184  copy_index=int(copy1),
185  residue_index=r1,
186  resolution=resolution).get_selected_particles()
187  ps2 = IMP.atom.Selection(root_hier,
188  state_index=nstate,
189  molecule=name2,
190  copy_index=int(copy2),
191  residue_index=r2,
192  resolution=resolution).get_selected_particles()
193 
194  ps1 = [IMP.atom.Hierarchy(p) for p in ps1]
195  ps2 = [IMP.atom.Hierarchy(p) for p in ps2]
196  else:
197  ps1 = IMP.pmi.tools.select(
198  r,
199  resolution=resolution,
200  name=c1,
201  name_is_ambiguous=False,
202  residue=r1)
203  ps2 = IMP.pmi.tools.select(
204  r,
205  resolution=resolution,
206  name=c2,
207  name_is_ambiguous=False,
208  residue=r2)
209 
210  if len(ps1) > 1:
211  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
212  elif len(ps1) == 0:
213  warnings.warn(
214  "CrossLinkingMassSpectrometryRestraint: "
215  "residue %d of chain %s is not there" % (r1, c1),
217  midb.write(str(xl) + "\n")
218  continue
219 
220  if len(ps2) > 1:
221  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
222  elif len(ps2) == 0:
223  warnings.warn(
224  "CrossLinkingMassSpectrometryRestraint: "
225  "residue %d of chain %s is not there" % (r2, c2),
227  midb.write(str(xl) + "\n")
228  continue
229 
230  p1 = ps1[0]
231  p2 = ps2[0]
232 
233  if p1 == p2 and r1 == r2:
234  warnings.warn(
235  "CrossLinkingMassSpectrometryRestraint: "
236  "same particle and same residue, skipping "
237  "cross-link", IMP.pmi.StructureWarning)
238  continue
239 
240  if new_contribution:
241  print("generating a new crosslink restraint")
242  new_contribution=False
244  self.model,
245  length,
246  slope)
247  restraints.append(dr)
248 
249 
250  if self.CrossLinkDataBase.sigma1_key not in xl.keys():
251  sigma1name="SIGMA"
252  xl[self.CrossLinkDataBase.sigma1_key]=sigma1name
253  else:
254  sigma1name=xl[self.CrossLinkDataBase.sigma1_key]
255  sigma1=self.create_sigma(sigma1name)
256 
257  if self.CrossLinkDataBase.sigma2_key not in xl.keys():
258  sigma2name="SIGMA"
259  xl[self.CrossLinkDataBase.sigma2_key]=sigma2name
260  else:
261  sigma2name=xl[self.CrossLinkDataBase.sigma2_key]
262  sigma2=self.create_sigma(sigma2name)
263 
264  if self.CrossLinkDataBase.psi_key not in xl.keys():
265  psiname="PSI"
266  xl[self.CrossLinkDataBase.psi_key]=psiname
267  else:
268  psiname=xl[self.CrossLinkDataBase.psi_key]
269 
270  psi=self.create_psi(psiname)
271 
272 
273  p1i = p1.get_particle_index()
274  xl["Particle1"]=p1
275  p2i = p2.get_particle_index()
276  xl["Particle2"]=p2
277  s1i = sigma1.get_particle().get_index()
278  xl["Particle_sigma1"]=sigma1
279  s2i = sigma2.get_particle().get_index()
280  xl["Particle_sigma2"]=sigma2
281  psii = psi.get_particle().get_index()
282  xl["Particle_psi"]=psi
283 
284  dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
285  xl["Restraint"]=dr
286 
287  print("--------------")
288  print("CrossLinkingMassSpectrometryRestraint: generating cross-link restraint between")
289  print("CrossLinkingMassSpectrometryRestraint: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
290  print("CrossLinkingMassSpectrometryRestraint: with sigma1 %s sigma2 %s psi %s" % (sigma1name, sigma2name, psiname))
291  print("CrossLinkingMassSpectrometryRestraint: between particles %s and %s" % (p1.get_name(), p2.get_name()))
292  print("==========================================\n")
293  for p, ex_xl in zip(IMP.pmi.tools._all_protocol_outputs(
294  representations, root_hier),
295  ex_xls):
296  p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
297  sigma1, sigma2, psi, ex_xl[1])
298 
299  # check if the two residues belong to the same rigid body
302  IMP.core.RigidMember(p1).get_rigid_body() ==
303  IMP.core.RigidMember(p2).get_rigid_body()):
304  xl["IntraRigidBody"]=True
305  else:
306  xl["IntraRigidBody"]=False
307 
308  xl_label=self.CrossLinkDataBase.get_short_cross_link_string(xl)
309  xl["ShortLabel"]=xl_label
310  dr.set_name(xl_label)
311 
312  if p1i != p2i:
313  pr = IMP.core.PairRestraint(self.model, dps2, (p1i, p2i))
314  pr.set_name(xl_label)
315  self.rslin.add_restraint(pr)
316 
317  self.xl_list.append(xl)
318 
319  indb.write(str(xl) + "\n")
320 
321  if len(self.xl_list) == 0:
322  raise SystemError("CrossLinkingMassSpectrometryRestraint: no crosslink was constructed")
323  self.xl_restraints = restraints
324  lw = IMP.isd.LogWrapper(restraints,1.0)
325  self.rs.add_restraint(lw)
326 
327  def __set_dataset(self, ds):
328  self.CrossLinkDataBase.dataset = ds
329  dataset = property(lambda self: self.CrossLinkDataBase.dataset,
330  __set_dataset)
331 
332  def get_hierarchies(self):
333  """ get the hierarchy """
334  return self.prot
335 
337  """ get the restraint set """
338  return self.rs
339 
340  def get_restraints(self):
341  """ get the restraints in a list """
342  return self.xl_restraints
343 
345  """ get the dummy restraints to be displayed in the rmf file """
346  return self.rslin
347 
349  """ Get a list of tuples containing the particle pairs """
350  ppairs = []
351  for i in range(len(self.pairs)):
352  p0 = self.pairs[i][0]
353  p1 = self.pairs[i][1]
354  ppairs.append((p0, p1))
355  return ppairs
356 
357  def set_output_level(self, level="low"):
358  """ Set the output level of the output """
359  self.outputlevel = level
360 
361  def set_psi_is_sampled(self, is_sampled=True):
362  """ Switch on/off the sampling of psi particles """
363  self.psi_is_sampled = is_sampled
364 
365  def set_sigma_is_sampled(self, is_sampled=True):
366  """ Switch on/off the sampling of sigma particles """
367  self.sigma_is_sampled = is_sampled
368 
369 
370  def create_sigma(self, name):
371  """ This is called internally. Creates a nuisance
372  on the structural uncertainty """
373  if name in self.sigma_dictionary:
374  return self.sigma_dictionary[name][0]
375 
376  sigmainit = 2.0
377  sigmaminnuis = 0.0000001
378  sigmamaxnuis = 1000.0
379  sigmamin = 0.01
380  sigmamax = 100.0
381  sigmatrans = 0.5
382  sigma = IMP.pmi.tools.SetupNuisance(self.model, sigmainit,
383  sigmaminnuis, sigmamaxnuis, self.sigma_is_sampled).get_particle()
384  self.sigma_dictionary[name] = (
385  sigma,
386  sigmatrans,
387  self.sigma_is_sampled)
388  self.rssig.add_restraint(
390  self.model,
391  sigma,
392  1000000000.0,
393  sigmamax,
394  sigmamin))
395  return sigma
396 
397  def create_psi(self, name):
398  """ This is called internally. Creates a nuisance
399  on the data uncertainty """
400  if name in self.psi_dictionary:
401  return self.psi_dictionary[name][0]
402 
403  psiinit=0.25
404  psiminnuis = 0.0000001
405  psimaxnuis = 0.4999999
406  psimin = 0.01
407  psimax = 0.49
408  psitrans = 0.1
409  psi = IMP.pmi.tools.SetupNuisance(self.model, psiinit,
410  psiminnuis, psimaxnuis,
411  self.psi_is_sampled).get_particle()
412  self.psi_dictionary[name] = (
413  psi,
414  psitrans,
415  self.psi_is_sampled)
416 
417  self.rspsi.add_restraint(
419  self.model,
420  psi,
421  1000000000.0,
422  psimax,
423  psimin))
424 
425  self.rspsi.add_restraint(IMP.isd.JeffreysRestraint(self.model, psi))
426  return psi
427 
428  def get_output(self):
429  """ Get the output of the restraint to be used by the IMP.pmi.output object"""
430  output = super(CrossLinkingMassSpectrometryRestraint, self).get_output()
431 
432  for xl in self.xl_list:
433 
434  xl_label=xl["ShortLabel"]
435  ln = xl["Restraint"]
436  p0 = xl["Particle1"]
437  p1 = xl["Particle2"]
438  output["CrossLinkingMassSpectrometryRestraint_Score_" +
439  xl_label] = str(-log(ln.unprotected_evaluate(None)))
440 
441  d0 = IMP.core.XYZ(p0)
442  d1 = IMP.core.XYZ(p1)
443  output["CrossLinkingMassSpectrometryRestraint_Distance_" +
444  xl_label] = str(IMP.core.get_distance(d0, d1))
445 
446 
447  for psiname in self.psi_dictionary:
448  output["CrossLinkingMassSpectrometryRestraint_Psi_" +
449  str(psiname) + self._label_suffix] = str(
450  self.psi_dictionary[psiname][0].get_scale())
451 
452  for sigmaname in self.sigma_dictionary:
453  output["CrossLinkingMassSpectrometryRestraint_Sigma_" +
454  str(sigmaname) + self._label_suffix] = str(
455  self.sigma_dictionary[sigmaname][0].get_scale())
456 
457 
458  return output
459 
460  def get_movers(self):
461  """ Get all need data to construct a mover in IMP.pmi.dof class"""
462  movers=[]
463  if self.sigma_is_sampled:
464  for sigmaname in self.sigma_dictionary:
465  mover_name="Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" + str(sigmaname) + "_" + self.label
466  particle=self.sigma_dictionary[sigmaname][0]
467  maxstep=(self.sigma_dictionary[sigmaname][1])
468  mv=IMP.core.NormalMover([particle],
469  IMP.FloatKeys([IMP.FloatKey("nuisance")]),maxstep)
470  mv.set_name(mover_name)
471  movers.append(mv)
472 
473  if self.psi_is_sampled:
474  for psiname in self.psi_dictionary:
475  mover_name="Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" + str(psiname) + "_" + self.label
476  particle=self.psi_dictionary[psiname][0]
477  maxstep=(self.psi_dictionary[psiname][1])
478  mv=IMP.core.NormalMover([particle],
479  IMP.FloatKeys([IMP.FloatKey("nuisance")]),maxstep)
480  mv.set_name(mover_name)
481  movers.append(mv)
482 
483  return movers
484 
485 
487  """ Get the particles to be sampled by the IMP.pmi.sampler object """
488  ps = {}
489  if self.sigma_is_sampled:
490  for sigmaname in self.sigma_dictionary:
491  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Sigma_" +
492  str(sigmaname) + self._label_suffix] =\
493  ([self.sigma_dictionary[sigmaname][0]],
494  self.sigma_dictionary[sigmaname][1])
495 
496  if self.psi_is_sampled:
497  for psiname in self.psi_dictionary:
498  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
499  str(psiname) + self._label_suffix] =\
500  ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
501 
502  return ps
503 
504 
506  """Setup cross-link distance restraints at atomic level
507  The "atomic" aspect is that it models the particle uncertainty with a Gaussian
508  The noise in the data and the structural uncertainty of cross-linked amino-acids
509  is inferred using Bayes' theory of probability
510  @note Wraps an IMP::isd::AtomicCrossLinkMSRestraint
511  """
512  def __init__(self,
513  root_hier,
514  xldb,
515  atom_type="NZ",
516  length=10.0,
517  slope=0.01,
518  nstates=None,
519  label=None,
520  nuisances_are_optimized=True,
521  sigma_init=5.0,
522  psi_init = 0.01,
523  one_psi=True,
524  filelabel=None,
525  weight=1.):
526  """Constructor.
527  Automatically creates one "sigma" per crosslinked residue and one "psis" per pair.
528  Other nuisance options are available.
529  @note Will return an error if the data+extra_sel don't specify two particles per XL pair.
530  @param root_hier The root hierarchy on which you'll do selection
531  @param xldb CrossLinkDataBase object
532  @param atom_type Can either be "NZ" or "CA"
533  @param length The XL linker length
534  @param slope Linear term to add to the restraint function to help when far away
535  @param nstates The number of states to model. Defaults to the number of states in root.
536  @param label The output label for the restraint
537  @param nuisances_are_optimized Whether to optimize nuisances
538  @param sigma_init The initial value for all the sigmas
539  @param psi_init The initial value for all the psis
540  @param one_psi Use a single psi for all restraints (if False, creates one per XL)
541  @param filelabel automatically generated file containing missing/included/excluded
542  cross-links will be labeled using this text
543  @param weight Weight of restraint
544 
545  """
546 
547  # basic params
548  self.root = root_hier
549  rname = "AtomicXLRestraint"
550  super(AtomicCrossLinkMSRestraint, self).__init__(
551  self.root.get_model(), name="AtomicXLRestraint", label=label,
552  weight=weight)
553  self.xldb = xldb
554  self.length = length
555  self.sigma_is_sampled = nuisances_are_optimized
556  self.psi_is_sampled = nuisances_are_optimized
557 
558  if atom_type not in("NZ","CA"):
559  raise Exception("AtomicXLRestraint: atom_type must be NZ or CA")
560  self.atom_type = atom_type
561  self.nstates = nstates
562  if nstates is None:
563  self.nstates = len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE))
564  elif nstates!=len(IMP.atom.get_by_type(self.root,IMP.atom.STATE_TYPE)):
565  print("Warning: nstates is not the same as the number of states in root")
566 
567  self.rs_psi = self._create_restraint_set("psi")
568  self.rs_sig = self._create_restraint_set("sigma")
569  self.rs_lin = self._create_restraint_set("linear")
570 
571  self.psi_dictionary = {}
572  self.sigma_dictionary = {}
573 
574  self.particles=defaultdict(set)
575  self.one_psi = one_psi
576  self.bonded_pairs = []
577  if self.one_psi:
578  print('creating a single psi for all XLs')
579  else:
580  print('creating one psi for each XL')
581 
582  # output logging file
583  if filelabel is not None:
584  indb = open("included." + filelabel + ".xl.db", "w")
585  exdb = open("excluded." + filelabel + ".xl.db", "w")
586  midb = open("missing." + filelabel + ".xl.db", "w")
587 
588 
589 
590  ### Setup nuisances
591  '''
592  # read ahead to get the number of XL's per residue
593  num_xls_per_res=defaultdict(int)
594  for unique_id in data:
595  for nstate in range(self.nstates):
596  for xl in data[unique_id]:
597  num_xls_per_res[str(xl.r1)]+=1
598  num_xls_per_res[str(xl.r2)]+=1
599 
600  # Will setup two sigmas based on promiscuity of the residue
601  sig_threshold=4
602  self.sig_low = setup_nuisance(self.m,self.rs_nuis,init_val=sigma_init,min_val=1.0,
603  max_val=100.0,is_opt=self.nuis_opt)
604  self.sig_high = setup_nuisance(self.m,self.rs_nuis,init_val=sigma_init,min_val=1.0,
605  max_val=100.0,is_opt=self.nuis_opt)
606  '''
607  self._create_sigma('sigma',sigma_init)
608  if one_psi:
609  self._create_psi('psi',psi_init)
610  else:
611  for xlid in self.xldb.xlid_iterator():
612  self._create_psi(xlid,psi_init)
613 
614  ### create all the XLs
615  xlrs=[]
616  for xlid in self.xldb.xlid_iterator():
617  # create restraint for this data point
618  if one_psi:
619  psip = self.psi_dictionary['psi'][0].get_particle_index()
620  else:
621  psip = self.psi_dictionary[unique_id][0].get_particle_index()
622  r = IMP.isd.AtomicCrossLinkMSRestraint(self.model,
623  self.length,
624  psip,
625  slope,
626  True)
627  num_contributions=0
628 
629  # add a contribution for each XL ambiguity option within each state
630  for nstate in self.nstates:
631  for xl in self.xldb[xlid]:
632  r1 = xl[self.xldb.residue1_key]
633  c1 = xl[self.xldb.protein1_key].strip()
634  r2 = xl[self.xldb.residue2_key]
635  c2 = xl[self.xldb.protein2_key].strip()
636 
637  # perform selection. these may contain multiples if Copies are used
638  ps1 = IMP.atom.Selection(self.root,
639  state_index=nstate,
640  atom_type = IMP.atom.AtomType(self.atom_type),
641  molecule=c1,
642  residue_index=r1).get_selected_particles()
643  ps2 = IMP.atom.Selection(self.root,
644  state_index=nstate,
645  atom_type = IMP.atom.AtomType(self.atom_type),
646  molecule=c2,
647  residue_index=r2).get_selected_particles()
648  if len(ps1) == 0:
649  warnings.warn("AtomicXLRestraint: residue %d of "
650  "chain %s is not there" % (r1, c1),
652  if filelabel is not None:
653  midb.write(str(xl) + "\n")
654  continue
655 
656  if len(ps2) == 0:
657  warnings.warn("AtomicXLRestraint: residue %d of "
658  "chain %s is not there" % (r2, c2),
660  if filelabel is not None:
661  midb.write(str(xl) + "\n")
662  continue
663 
664 
665  # figure out sig1 and sig2 based on num XLs
666  '''
667  num1=num_xls_per_res[str(xl.r1)]
668  num2=num_xls_per_res[str(xl.r2)]
669  if num1<sig_threshold:
670  sig1=self.sig_low
671  else:
672  sig1=self.sig_high
673  if num2<sig_threshold:
674  sig2=self.sig_low
675  else:
676  sig2=self.sig_high
677  '''
678  sig1 = self.sigma_dictionary['sigma'][0]
679  sig2 = self.sigma_dictionary['sigma'][0]
680 
681  # add each copy contribution to restraint
682  for p1,p2 in itertools.product(ps1,ps2):
683  if p1==p2:
684  continue
685  self.particles[nstate]|=set((p1,p2))
686  r.add_contribution([p1.get_index(),p2.get_index()],
687  [sig1.get_particle_index(),sig2.get_particle_index()])
688  num_contributions+=1
689 
690  if num_contributions>0:
691  print('XL:',xlid,'num contributions:',num_contributions)
692  xlrs.append(r)
693  if len(xlrs)==0:
694  raise Exception("You didn't create any XL restraints")
695  print('created',len(xlrs),'XL restraints')
696  rname = self.rs.get_name()
697  self.rs=IMP.isd.LogWrapper(xlrs, self.weight)
698  self.rs.set_name(rname)
699  self.rs.set_weight(self.weight)
700  self.restraint_sets = [self.rs] + self.restraint_sets[1:]
701 
702  def get_hierarchy(self):
703  return self.prot
704 
705  def _create_sigma(self, name,sigmainit):
706  """ This is called internally. Creates a nuisance
707  on the structural uncertainty """
708  if name in self.sigma_dictionary:
709  return self.sigma_dictionary[name][0]
710 
711  sigmaminnuis = 0.0000001
712  sigmamaxnuis = 1000.0
713  sigmamin = 0.01
714  sigmamax = 100.0
715  sigmatrans = 0.5
716  sigma = IMP.pmi.tools.SetupNuisance(self.model,
717  sigmainit,
718  sigmaminnuis,
719  sigmamaxnuis,
720  self.sigma_is_sampled).get_particle()
721  self.sigma_dictionary[name] = (
722  sigma,
723  sigmatrans,
724  self.sigma_is_sampled)
725  self.rs_sig.add_restraint(
727  self.model,
728  sigma,
729  1000000000.0,
730  sigmamax,
731  sigmamin))
732  return sigma
733 
734  def _create_psi(self, name,psiinit):
735  """ This is called internally. Creates a nuisance
736  on the data uncertainty """
737  if name in self.psi_dictionary:
738  return self.psi_dictionary[name][0]
739 
740  psiminnuis = 0.0000001
741  psimaxnuis = 0.4999999
742  psimin = 0.01
743  psimax = 0.49
744  psitrans = 0.1
745  psi = IMP.pmi.tools.SetupNuisance(self.model,
746  psiinit,
747  psiminnuis,
748  psimaxnuis,
749  self.psi_is_sampled).get_particle()
750  self.psi_dictionary[name] = (
751  psi,
752  psitrans,
753  self.psi_is_sampled)
754 
755  self.rs_psi.add_restraint(
757  self.model,
758  psi,
759  1000000000.0,
760  psimax,
761  psimin))
762 
763  self.rs_psi.add_restraint(IMP.isd.JeffreysRestraint(self.model, psi))
764  return psi
765 
767  """ create dummy harmonic restraints for each XL but don't add to model
768  Makes it easy to see each contribution to each XL in RMF
769  """
770  class MyGetRestraint(object):
771  def __init__(self,rs):
772  self.rs=rs
773  def get_restraint_for_rmf(self):
774  return self.rs
775 
776  dummy_model=IMP.Model()
777  hps = IMP.core.HarmonicDistancePairScore(self.length,1.0)
778  dummy_rs=[]
779  for nxl in range(self.get_number_of_restraints()):
780  xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
781  rs = IMP.RestraintSet(dummy_model, 'atomic_xl_'+str(nxl))
782  for ncontr in range(xl.get_number_of_contributions()):
783  ps=xl.get_contribution(ncontr)
784  dr = IMP.core.PairRestraint(hps,[self.model.get_particle(p) for p in ps],
785  'xl%i_contr%i'%(nxl,ncontr))
786  rs.add_restraint(dr)
787  dummy_rs.append(MyGetRestraint(rs))
788  return dummy_rs
789 
790 
792  """ Get the particles to be sampled by the IMP.pmi.sampler object """
793  ps = {}
794  if self.sigma_is_sampled:
795  for sigmaname in self.sigma_dictionary:
796  ps["Nuisances_AtomicCrossLinkingMSRestraint_Sigma_" +
797  str(sigmaname) + self._label_suffix] = \
798  ([self.sigma_dictionary[sigmaname][0]],
799  self.sigma_dictionary[sigmaname][1])
800  if self.psi_is_sampled:
801  for psiname in self.psi_dictionary:
802  ps["Nuisances_CrossLinkingMassSpectrometryRestraint_Psi_" +
803  str(psiname) + self._label_suffix] =\
804  ([self.psi_dictionary[psiname][0]], self.psi_dictionary[psiname][1])
805  return ps
806 
807  def get_bonded_pairs(self):
808  return self.bonded_pairs
809 
810  def get_number_of_restraints(self):
811  return self.rs.get_number_of_restraints()
812 
813  def __repr__(self):
814  return 'XL restraint with '+str(len(self.rs.get_restraint(0).get_number_of_restraints())) \
815  + ' data points'
816 
817  def load_nuisances_from_stat_file(self,in_fn,nframe):
818  """Read a stat file and load all the sigmas.
819  This is potentially quite stupid.
820  It's also a hack since the sigmas should be stored in the RMF file.
821  Also, requires one sigma and one psi for ALL XLs.
822  """
823  import subprocess
824  sig_val = float(subprocess.check_output(["process_output.py","-f",in_fn,
825  "-s","AtomicXLRestraint_sigma"]).split('\n>')[1+nframe])
826  psi_val = float(subprocess.check_output(["process_output.py","-f",in_fn,
827  "-s","AtomicXLRestraint_psi"]).split('\n>')[1+nframe])
828  for nxl in range(self.rs.get_number_of_restraints()):
829  xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
830  psip = xl.get_psi()
831  IMP.isd.Scale(self.model,psip).set_scale(psi_val)
832  for contr in range(xl.get_number_of_contributions()):
833  sig1,sig2=xl.get_contribution_sigmas(contr)
834  IMP.isd.Scale(self.model,sig1).set_scale(sig_val)
835 
836  print('loaded nuisances from file')
837 
838  def plot_violations(self,out_prefix,
839  max_prob_for_violation=0.1,
840  min_dist_for_violation=1e9,
841  coarsen=False,
842  limit_to_chains=None,
843  exclude_chains=''):
844  """Create CMM files, one for each state, of all crosslinks.
845  will draw in GREEN if non-violated in all states (or if only one state)
846  will draw in PURPLE if non-violated only in a subset of states (draws nothing elsewhere)
847  will draw in RED in ALL states if all violated
848  (if only one state, you'll only see green and red)
849 
850  @param out_prefix Output xlink files prefix
851  @param max_prob_for_violation It's a violation if the probability is below this
852  @param min_dist_for_violation It's a violation if the min dist is above this
853  @param coarsen Use CA positions
854  @param limit_to_chains Try to visualize just these chains
855  @param exclude_to_chains Try to NOT visualize these chains
856  """
857  print('going to calculate violations and plot CMM files')
858  all_stats = self.get_best_stats()
859  all_dists = [s["low_dist"] for s in all_stats]
860 
861  # prepare one output file per state
862  out_fns=[]
863  out_nvs=[]
864  state_info=[]
865  cmds = defaultdict(set)
866  for nstate in range(self.nstates):
867  outf=open(out_prefix+str(nstate)+'.cmm','w')
868  outf.write('<marker_set name="xlinks_state%i"> \n' % nstate)
869  out_fns.append(outf)
870  out_nvs.append(0)
871  print('will limit to',limit_to_chains)
872  print('will exclude',exclude_chains)
873  state_info.append(self.get_best_stats(nstate,
874  limit_to_chains,
875  exclude_chains))
876 
877  for nxl in range(self.rs.get_number_of_restraints()):
878  # for this XL, check which states passed
879  npass=[]
880  nviol=[]
881  for nstate in range(self.nstates):
882  prob = state_info[nstate][nxl]["prob"]
883  low_dist = state_info[nstate][nxl]["low_dist"]
884  if prob<max_prob_for_violation or low_dist>min_dist_for_violation:
885  nviol.append(nstate)
886  else:
887  npass.append(nstate)
888 
889  # special coloring when all pass or all fail
890  all_pass=False
891  all_viol=False
892  if len(npass)==self.nstates:
893  all_pass=True
894  elif len(nviol)==self.nstates:
895  all_viol=True
896 
897  # finally, color based on above info
898  print(nxl,'state dists:',[state_info[nstate][nxl]["low_dist"] for nstate in range(self.nstates)],
899  'viol states:',nviol,'all viol?',all_viol)
900  for nstate in range(self.nstates):
901  if all_pass:
902  r=0.365; g=0.933; b=0.365;
903  continue
904  elif all_viol:
905  r=0.980; g=0.302; b=0.247;
906  continue
907  else:
908  if nstate in nviol:
909  continue
910  else:
911  #r=0.9; g=0.34; b=0.9;
912  r=0.365; g=0.933; b=0.365;
913  # now only showing if UNIQUELY PASSING in this state
914  pp = state_info[nstate][nxl]["low_pp"]
915  c1=IMP.core.XYZ(self.model,pp[0]).get_coordinates()
916  c2=IMP.core.XYZ(self.model,pp[1]).get_coordinates()
917 
918  r1 = IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[0])).get_index()
919  ch1 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[0]))
920  r2 = IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[0])).get_index()
921  ch2 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[0]))
922 
923  cmds[nstate].add((ch1,r1))
924  cmds[nstate].add((ch2,r2))
925 
926  outf = out_fns[nstate]
927  nv = out_nvs[nstate]
928  outf.write('<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
929  'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,c1[0],c1[1],c1[2],r,g,b))
930  outf.write('<marker id= "%d" x="%.3f" y="%.3f" z="%.3f" radius="0.8" '
931  'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv+1,c2[0],c2[1],c2[2],r,g,b))
932  outf.write('<link id1= "%d" id2="%d" radius="0.8" '
933  'r="%.2f" g="%.2f" b="%.2f"/> \n' % (nv,nv+1,r,g,b))
934  out_nvs[nstate]+=2
935 
936  for nstate in range(self.nstates):
937  out_fns[nstate].write('</marker_set>\n')
938  out_fns[nstate].close()
939  cmd = ''
940  for ch,r in cmds[nstate]:
941  cmd+='#%i:%i.%s '%(nstate,r,ch)
942  print(cmd)
943 
944  return all_dists
945  def _get_contribution_info(self,xl,ncontr,use_CA=False):
946  """Return the particles at that contribution. If requested will return CA's instead"""
947  idx1=xl.get_contribution(ncontr)[0]
948  idx2=xl.get_contribution(ncontr)[1]
949  if use_CA:
950  idx1 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,idx1)),
951  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
952  idx2 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,idx2)),
953  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
954  dist = IMP.algebra.get_distance(IMP.core.XYZ(self.model,idx1).get_coordinates(),
955  IMP.core.XYZ(self.model,idx2).get_coordinates())
956  return idx1,idx2,dist
957 
958  def get_best_stats(self,limit_to_state=None,limit_to_chains=None,exclude_chains='',use_CA=False):
959  ''' return the probability, best distance, two coords, and possibly the psi for each xl
960  @param limit_to_state Only examine contributions from one state
961  @param limit_to_chains Returns the particles for certain "easy to visualize" chains
962  @param exclude_chains Even if you limit, don't let one end be in this list.
963  Only works if you also limit chains
964  @param use_CA Limit to CA particles
965  '''
966  ret = []
967  for nxl in range(self.rs.get_number_of_restraints()):
968  this_info = {}
969  xl=IMP.isd.AtomicCrossLinkMSRestraint.get_from(self.rs.get_restraint(nxl))
970  low_dist=1e6
971  low_contr = None
972  low_pp = None
973  state_contrs=[]
974  low_pp_lim = None
975  low_dist_lim = 1e6
976  for contr in range(xl.get_number_of_contributions()):
977  pp = xl.get_contribution(contr)
978  if use_CA:
979  idx1 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[0])),
980  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
981  idx2 = IMP.atom.Selection(IMP.atom.get_residue(IMP.atom.Atom(self.model,pp[1])),
982  atom_type=IMP.atom.AtomType("CA")).get_selected_particle_indexes()[0]
983  pp = [idx1,idx2]
984  if limit_to_state is not None:
985  nstate = IMP.atom.get_state_index(IMP.atom.Atom(self.model,pp[0]))
986  if nstate!=limit_to_state:
987  continue
988  state_contrs.append(contr)
989  dist = IMP.core.get_distance(IMP.core.XYZ(self.model,pp[0]),
990  IMP.core.XYZ(self.model,pp[1]))
991  if limit_to_chains is not None:
992  c1 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[0]))
993  c2 = IMP.atom.get_chain_id(IMP.atom.Atom(self.model,pp[1]))
994  if (c1 in limit_to_chains or c2 in limit_to_chains) and (
995  c1 not in exclude_chains and c2 not in exclude_chains):
996  if dist<low_dist_lim:
997  low_dist_lim = dist
998  low_pp_lim = pp
999  if dist<low_dist:
1000  low_dist = dist
1001  low_contr = contr
1002  low_pp = pp
1003  if limit_to_state is not None:
1004  this_info["prob"] = xl.evaluate_for_contributions(state_contrs,None)
1005  else:
1006  this_info["prob"] = xl.unprotected_evaluate(None)
1007  if limit_to_chains is not None:
1008  this_info["low_pp"] = low_pp_lim
1009  else:
1010  this_info["low_pp"] = low_pp
1011 
1012  this_info["low_dist"] = low_dist
1013  if not self.one_psi:
1014  pval = IMP.isd.Scale(self.model,xl.get_psi()).get_scale()
1015  this_info["psi"] = pval
1016  ret.append(this_info)
1017  return ret
1018 
1019  def print_stats(self):
1020  #print("XL restraint statistics\n<num> <prob> <bestdist> <sig1> <sig2> <is_viol>")
1021  stats = self.get_best_stats()
1022  for nxl,s in enumerate(stats):
1023  #print('%i %.4f %.4f %.4f %.4f %i'%(nxl,prob,low_dist,sig1,sig2,is_viol))
1024  print(s["low_dist"])
1025 
1026 
1027  def get_output(self):
1028  output = super(AtomicCrossLinkMSRestraint, self).get_output()
1029 
1030  ### HACK to make it easier to see the few sigmas
1031  #output["AtomicXLRestraint_sigma"] = self.sigma.get_scale()
1032  #if self.one_psi:
1033  # output["AtomicXLRestraint_psi"] = self.psi.get_scale()
1034  ######
1035 
1036  # count distances above length
1037  bad_count=0
1038  stats = self.get_best_stats()
1039  for nxl,s in enumerate(stats):
1040  if s['low_dist']>20.0:
1041  bad_count+=1
1042  output["AtomicXLRestraint_%i_%s"%(nxl,"Prob")]=str(s['prob'])
1043  output["AtomicXLRestraint_%i_%s"%(nxl,"BestDist")]=str(s['low_dist'])
1044  if not self.one_psi:
1045  output["AtomicXLRestraint_%i_%s"%(nxl,"psi")]=str(s['psi'])
1046  output["AtomicXLRestraint_NumViol"] = str(bad_count)
1047  return output
1048 
1049 
1050 @IMP.deprecated_object("2.5", "Use IMP.pmi.restraints.crosslinking.CrossLinkingMassSpectrometryRestraint instead.")
1051 class ISDCrossLinkMS(IMP.pmi.restraints._NuisancesBase):
1052  def __init__(self, representation,
1053  restraints_file,
1054  length,
1055  jackknifing=None,
1056  # jackknifing (Float [0,1]) default=None; percentage of data to be
1057  # removed randomly
1058  resolution=None,
1059  # resolution (Non-negative Integer) default=None; percentage of data
1060  # to be removed randomly
1061  slope=0.0,
1062  inner_slope=0.0,
1063  columnmapping=None,
1064  rename_dict=None,
1065  offset_dict=None,
1066  csvfile=False,
1067  ids_map=None,
1068  radius_map=None,
1069  filters=None,
1070  label="None",
1071  filelabel="None",
1072  automatic_sigma_classification=False,
1073  attributes_for_label=None):
1074 
1075  # columnindexes is a list of column indexes for protein1, protein2, residue1, residue2,idscore, XL unique id
1076  # by default column 0 = protein1; column 1 = protein2; column 2 = residue1; column 3 = residue2;
1077  # column 4 = idscores
1078  # attributes_for_label: anything in the csv database that must be added to the label
1079  # slope is the slope defined on the linear function
1080  # inner_slope is the slope defined on the restraint directly
1081  # suggestion: do not use both!
1082 
1083  if type(representation) != list:
1084  representations = [representation]
1085  else:
1086  representations = representation
1087 
1088  if columnmapping is None:
1089  columnmapping = {}
1090  columnmapping["Protein1"] = 0
1091  columnmapping["Protein2"] = 1
1092  columnmapping["Residue1"] = 2
1093  columnmapping["Residue2"] = 3
1094  columnmapping["IDScore"] = 4
1095  columnmapping["XLUniqueID"] = 5
1096 
1097  if csvfile:
1098  # in case the file is a csv file
1099  # columnmapping will contain the field names
1100  # that compare in the first line of the csv file
1101  db = IMP.pmi.tools.get_db_from_csv(restraints_file)
1102  else:
1103  db = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
1104 
1105  indb = open("included." + filelabel + ".xl.db", "w")
1106  exdb = open("excluded." + filelabel + ".xl.db", "w")
1107  midb = open("missing." + filelabel + ".xl.db", "w")
1108 
1109  self.m = representations[0].prot.get_model()
1110  self.rs = IMP.RestraintSet(self.m, 'data')
1111  self.rspsi = IMP.RestraintSet(self.m, 'prior_psi')
1112  self.rssig = IMP.RestraintSet(self.m, 'prior_sigmas')
1113  self.rslin = IMP.RestraintSet(self.m, 'prior_linear')
1114  self.rslen = IMP.RestraintSet(self.m, 'prior_length')
1115 
1116  self.weight = 1.0
1117  self.label = label
1118  self.pairs = []
1119  self.sigma_dictionary = {}
1120  self.psi_dictionary = {}
1121  self.psi_is_sampled = True
1122  self.sigma_is_sampled = True
1123 
1124  if os.path.exists(restraints_file):
1125  l = ihm.location.InputFileLocation(restraints_file,
1126  details="Crosslinks")
1127  self.dataset = ihm.dataset.CXMSDataset(l)
1128  else:
1129  self.dataset = None
1130 
1131  xl_groups = [p.get_cross_link_group(self)
1132  for p, state in representations[0]._protocol_output]
1133 
1134  # isd_map is a dictionary/map that is used to determine the psi
1135  # parameter from identity scores (such as ID-Score, or FDR)
1136  if ids_map is None:
1137  self.ids_map = IMP.pmi.tools.map()
1138  self.ids_map.set_map_element(20.0, 0.05)
1139  self.ids_map.set_map_element(65.0, 0.01)
1140  else:
1141  self.ids_map = ids_map
1142 
1143  if radius_map is None:
1144  self.radius_map = IMP.pmi.tools.map()
1145  if automatic_sigma_classification:
1146  self.radius_map.set_map_element(10, 10)
1147  self.radius_map.set_map_element(1, 1)
1148  else:
1149  self.radius_map = radius_map
1150 
1151  self.outputlevel = "low"
1152 
1153  # small linear contribution for long range
1154  self.linear = IMP.core.Linear(0, 0.0)
1155  self.linear.set_slope(slope)
1156  dps2 = IMP.core.DistancePairScore(self.linear)
1157 
1158  protein1 = columnmapping["Protein1"]
1159  protein2 = columnmapping["Protein2"]
1160  residue1 = columnmapping["Residue1"]
1161  residue2 = columnmapping["Residue2"]
1162  idscore = columnmapping["IDScore"]
1163  try:
1164  xluniqueid = columnmapping["XLUniqueID"]
1165  except:
1166  xluniqueid = None
1167 
1168  restraints = []
1169 
1170  # we need this dictionary to create ambiguity (i.e., multistate)
1171  # if one id is already present in the dictionary, add the term to the
1172  # corresponding already generated restraint
1173 
1174  uniqueid_restraints_map = {}
1175 
1176  for nxl, entry in enumerate(db):
1177 
1178  if not jackknifing is None:
1179 
1180  # to be implemented
1181  # the problem is that in the replica exchange
1182  # you have to broadcast the same restraints to every
1183  # replica
1184 
1185  raise NotImplementedError("jackknifing not yet implemented")
1186 
1187  if not csvfile:
1188  tokens = entry.split()
1189  if len(tokens)==0:
1190  continue
1191 
1192  # skip character
1193  if (tokens[0] == "#"):
1194  continue
1195  try:
1196  r1 = int(tokens[residue1])
1197  c1 = tokens[protein1]
1198  r2 = int(tokens[residue2])
1199  c2 = tokens[protein2]
1200 
1201  if offset_dict is not None:
1202  if c1 in offset_dict: r1+=offset_dict[c1]
1203  if c2 in offset_dict: r2+=offset_dict[c2]
1204 
1205  if rename_dict is not None:
1206  if c1 in rename_dict: c1=rename_dict[c1]
1207  if c2 in rename_dict: c2=rename_dict[c2]
1208 
1209  if idscore is None:
1210  ids = 1.0
1211  else:
1212  ids = float(tokens[idscore])
1213  if xluniqueid is None:
1214  xlid = str(nxl)
1215  else:
1216  xlid = tokens[xluniqueid]
1217  except:
1218  print("this line was not accessible " + str(entry))
1219  if residue1 not in entry: print(str(residue1)+" keyword not in database")
1220  if residue2 not in entry: print(str(residue2)+" keyword not in database")
1221  if protein1 not in entry: print(str(protein1)+" keyword not in database")
1222  if protein2 not in entry: print(str(protein2)+" keyword not in database")
1223  if idscore not in entry: print(str(idscore)+" keyword not in database")
1224  if xluniqueid not in entry: print(str(xluniqueid)+" keyword not in database")
1225  continue
1226 
1227  else:
1228  if filters is not None:
1229  if eval(IMP.pmi.tools.cross_link_db_filter_parser(filters)) == False:
1230  exdb.write(str(entry) + "\n")
1231  continue
1232 
1233  try:
1234  r1 = int(entry[residue1])
1235  c1 = entry[protein1]
1236  r2 = int(entry[residue2])
1237  c2 = entry[protein2]
1238 
1239  if offset_dict is not None:
1240  if c1 in offset_dict: r1+=offset_dict[c1]
1241  if c2 in offset_dict: r2+=offset_dict[c2]
1242 
1243  if rename_dict is not None:
1244  if c1 in rename_dict: c1=rename_dict[c1]
1245  if c2 in rename_dict: c2=rename_dict[c2]
1246 
1247  if idscore is None:
1248  ids = 1.0
1249  else:
1250  try:
1251  ids = float(entry[idscore])
1252  except ValueError:
1253  ids = entry[idscore]
1254  if xluniqueid is None:
1255  xlid = str(nxl)
1256  else:
1257  xlid = entry[xluniqueid]
1258 
1259  except:
1260  print("this line was not accessible " + str(entry))
1261  if residue1 not in entry: print(str(residue1)+" keyword not in database")
1262  if residue2 not in entry: print(str(residue2)+" keyword not in database")
1263  if protein1 not in entry: print(str(protein1)+" keyword not in database")
1264  if protein2 not in entry: print(str(protein2)+" keyword not in database")
1265  if idscore not in entry: print(str(idscore)+" keyword not in database")
1266  if xluniqueid not in entry: print(str(xluniqueid)+" keyword not in database")
1267  continue
1268 
1269  # todo: check that offset is handled correctly
1270  ex_xls = [(p[0].add_experimental_cross_link(r1, c1, r2, c2,
1271  group), group)
1272  for p, group in zip(representations[0]._protocol_output,
1273  xl_groups)]
1274 
1275  for nstate, r in enumerate(representations):
1276  # loop over every state
1277 
1278  ps1 = IMP.pmi.tools.select(
1279  r,
1280  resolution=resolution,
1281  name=c1,
1282  name_is_ambiguous=False,
1283  residue=r1)
1284  ps2 = IMP.pmi.tools.select(
1285  r,
1286  resolution=resolution,
1287  name=c2,
1288  name_is_ambiguous=False,
1289  residue=r2)
1290 
1291  if len(ps1) > 1:
1292  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r1, c1, str(ps1)))
1293  elif len(ps1) == 0:
1294  warnings.warn("ISDCrossLinkMS: residue %d of chain %s "
1295  "is not there" % (r1, c1),
1297  midb.write(str(entry) + "\n")
1298  continue
1299 
1300  if len(ps2) > 1:
1301  raise ValueError("residue %d of chain %s selects multiple particles %s" % (r2, c2, str(ps2)))
1302  elif len(ps2) == 0:
1303  warnings.warn("ISDCrossLinkMS: residue %d of chain %s "
1304  "is not there" % (r2, c2),
1306  midb.write(str(entry) + "\n")
1307  continue
1308 
1309  p1 = ps1[0]
1310  p2 = ps2[0]
1311 
1312  if (p1 == p2) and (r1 == r2) :
1313  warnings.warn("ISDCrossLinkMS Restraint: on the identical "
1314  "bead particles and the identical residues, "
1315  "thus skipping this cross-link.",
1317  continue
1318 
1319  if xlid in uniqueid_restraints_map:
1320  print("Appending a crosslink restraint into the uniqueID %s" % str(xlid))
1321  dr = uniqueid_restraints_map[xlid]
1322  else:
1323  print("Generating a NEW crosslink restraint with a uniqueID %s" % str(xlid))
1325  self.m,
1326  length,
1327  inner_slope)
1328  restraints.append(dr)
1329  uniqueid_restraints_map[xlid] = dr
1330 
1331  mappedr1 = self.radius_map.get_map_element(
1332  IMP.pmi.Uncertainty(p1).get_uncertainty())
1333  sigma1 = self.get_sigma(mappedr1)[0]
1334  mappedr2 = self.radius_map.get_map_element(
1335  IMP.pmi.Uncertainty(p2).get_uncertainty())
1336  sigma2 = self.get_sigma(mappedr2)[0]
1337 
1338  psival = self.ids_map.get_map_element(ids)
1339  psi = self.get_psi(psival)[0]
1340 
1341 
1342  p1i = p1.get_particle_index()
1343  p2i = p2.get_particle_index()
1344  s1i = sigma1.get_particle().get_index()
1345  s2i = sigma2.get_particle().get_index()
1346 
1347  #print nstate, p1i, p2i, p1.get_name(), p2.get_name()
1348 
1349  psii = psi.get_particle().get_index()
1350 
1351  dr.add_contribution((p1i, p2i), (s1i, s2i), psii)
1352  print("--------------")
1353  print("ISDCrossLinkMS: generating cross-link restraint between")
1354  print("ISDCrossLinkMS: residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
1355  print("ISDCrossLinkMS: with sigma1 %f sigma2 %f psi %s" % (mappedr1, mappedr2, psival))
1356  print("ISDCrossLinkMS: between particles %s and %s" % (p1.get_name(), p2.get_name()))
1357  print("==========================================\n")
1358  indb.write(str(entry) + "\n")
1359  for p, ex_xl in zip(representations[0]._protocol_output,
1360  ex_xls):
1361  p[0].add_cross_link(p[1], ex_xl[0], p1, p2, length,
1362  sigma1, sigma2, psi, ex_xl[1])
1363 
1364  # check if the two residues belong to the same rigid body
1367  IMP.core.RigidMember(p1).get_rigid_body() ==
1368  IMP.core.RigidMember(p2).get_rigid_body()):
1369  xlattribute = "intrarb"
1370  else:
1371  xlattribute = "interrb"
1372 
1373  if csvfile:
1374  if not attributes_for_label is None:
1375  for a in attributes_for_label:
1376  xlattribute = xlattribute + "_" + str(entry[a])
1377 
1378  xlattribute = xlattribute + "-State:" + str(nstate)
1379 
1380  dr.set_name(
1381  xlattribute + "-" + c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2) + "_" + self.label)
1382 
1383  if p1i != p2i:
1384  pr = IMP.core.PairRestraint(self.m, dps2, (p1i, p2i))
1385  pr.set_name(
1386  xlattribute + "-" + c1 + ":" + str(r1) + "-" + c2 + ":" + str(r2) + "_" + self.label)
1387  self.rslin.add_restraint(pr)
1388 
1389 
1390  self.pairs.append(
1391  (p1,
1392  p2,
1393  dr,
1394  r1,
1395  c1,
1396  r2,
1397  c2,
1398  xlattribute,
1399  mappedr1,
1400  mappedr2,
1401  psival,
1402  xlid,
1403  nstate,
1404  ids))
1405 
1406  lw = IMP.isd.LogWrapper(restraints, self.weight)
1407  self.rs.add_restraint(lw)
1408 
1409  # Provide self.model for compatibility with newer code
1410  model = property(lambda s: s.m)
1411 
1412  def set_weight(self, weight):
1413  self.weight = weight
1414  self.rs.set_weight(weight)
1415 
1416  def set_slope_linear_term(self, slope):
1417  self.linear.set_slope(slope)
1418 
1419  def set_label(self, label):
1420  self.label = label
1421 
1422  def add_to_model(self):
1423  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1424  IMP.pmi.tools.add_restraint_to_model(self.m, self.rspsi)
1425  IMP.pmi.tools.add_restraint_to_model(self.m, self.rssig)
1426  IMP.pmi.tools.add_restraint_to_model(self.m, self.rslen)
1427  IMP.pmi.tools.add_restraint_to_model(self.m, self.rslin)
1428 
1429  def get_hierarchies(self):
1430  return self.prot
1431 
1432  def get_restraint_sets(self):
1433  return self.rs
1434 
1435  def get_restraint(self):
1436  return self.rs
1437 
1438  def get_restraint_for_rmf(self):
1439  return self.rslin
1440 
1441  def get_restraints(self):
1442  rlist = []
1443  for r in self.rs.get_restraints():
1444  rlist.append(IMP.core.PairRestraint.get_from(r))
1445  return rlist
1446 
1447  def get_particle_pairs(self):
1448  ppairs = []
1449  for i in range(len(self.pairs)):
1450  p0 = self.pairs[i][0]
1451  p1 = self.pairs[i][1]
1452  ppairs.append((p0, p1))
1453  return ppairs
1454 
1455  def set_output_level(self, level="low"):
1456  # this might be "low" or "high"
1457  self.outputlevel = level
1458 
1459  def set_psi_is_sampled(self, is_sampled=True):
1460  self.psi_is_sampled = is_sampled
1461 
1462  def set_sigma_is_sampled(self, is_sampled=True):
1463  self.sigma_is_sampled = is_sampled
1464 
1465  def get_label(self,pairs_index):
1466  resid1 = self.pairs[pairs_index][3]
1467  chain1 = self.pairs[pairs_index][4]
1468  resid2 = self.pairs[pairs_index][5]
1469  chain2 = self.pairs[pairs_index][6]
1470  attribute = self.pairs[pairs_index][7]
1471  rad1 = self.pairs[pairs_index][8]
1472  rad2 = self.pairs[pairs_index][9]
1473  psi = self.pairs[pairs_index][10]
1474  xlid= self.pairs[pairs_index][11]
1475  label = attribute + "-" + \
1476  str(resid1) + ":" + chain1 + "_" + str(resid2) + ":" + \
1477  chain2 + "-" + str(rad1) + "-" + str(rad2) + "-" + str(psi)
1478  return label
1479 
1480  def write_db(self,filename):
1481  cldb=IMP.pmi.output.CrossLinkIdentifierDatabase()
1482 
1483  for pairs_index in range(len(self.pairs)):
1484 
1485  resid1 = self.pairs[pairs_index][3]
1486  chain1 = self.pairs[pairs_index][4]
1487  resid2 = self.pairs[pairs_index][5]
1488  chain2 = self.pairs[pairs_index][6]
1489  attribute = self.pairs[pairs_index][7]
1490  rad1 = self.pairs[pairs_index][8]
1491  rad2 = self.pairs[pairs_index][9]
1492  psi = self.pairs[pairs_index][10]
1493  xlid= self.pairs[pairs_index][11]
1494  nstate=self.pairs[pairs_index][12]
1495  ids=self.pairs[pairs_index][13]
1496 
1497  label=self.get_label(pairs_index)
1498  cldb.set_unique_id(label,xlid)
1499  cldb.set_protein1(label,chain1)
1500  cldb.set_protein2(label,chain2)
1501  cldb.set_residue1(label,resid1)
1502  cldb.set_residue2(label,resid2)
1503  cldb.set_idscore(label,ids)
1504  cldb.set_state(label,nstate)
1505  cldb.set_sigma1(label,rad1)
1506  cldb.set_sigma2(label,rad2)
1507  cldb.set_psi(label,psi)
1508  cldb.write(filename)
1509 
1510 
1511  def get_output(self):
1512  # content of the crosslink database pairs
1513  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
1514  output = {}
1515  score = self.weight * self.rs.unprotected_evaluate(None)
1516  output["_TotalScore"] = str(score)
1517  output["ISDCrossLinkMS_Data_Score_" + self.label] = str(score)
1518  output["ISDCrossLinkMS_PriorSig_Score_" +
1519  self.label] = self.rssig.unprotected_evaluate(None)
1520  output["ISDCrossLinkMS_PriorPsi_Score_" +
1521  self.label] = self.rspsi.unprotected_evaluate(None)
1522  output["ISDCrossLinkMS_Linear_Score_" +
1523  self.label] = self.rslin.unprotected_evaluate(None)
1524  for i in range(len(self.pairs)):
1525 
1526  label=self.get_label(i)
1527  ln = self.pairs[i][2]
1528  p0 = self.pairs[i][0]
1529  p1 = self.pairs[i][1]
1530  output["ISDCrossLinkMS_Score_" +
1531  label + "_" + self.label] = str(self.weight * -log(ln.unprotected_evaluate(None)))
1532 
1533  d0 = IMP.core.XYZ(p0)
1534  d1 = IMP.core.XYZ(p1)
1535  output["ISDCrossLinkMS_Distance_" +
1536  label + "_" + self.label] = str(IMP.core.get_distance(d0, d1))
1537 
1538 
1539  for psiindex in self.psi_dictionary:
1540  output["ISDCrossLinkMS_Psi_" +
1541  str(psiindex) + "_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
1542 
1543  for resolution in self.sigma_dictionary:
1544  output["ISDCrossLinkMS_Sigma_" +
1545  str(resolution) + "_" + self.label] = str(self.sigma_dictionary[resolution][0].get_scale())
1546 
1547 
1548  return output
1549 
1550  def get_particles_to_sample(self):
1551  ps = {}
1552  for resolution in self.sigma_dictionary:
1553  if self.sigma_dictionary[resolution][2] and self.sigma_is_sampled:
1554  ps["Nuisances_ISDCrossLinkMS_Sigma_" + str(resolution) + "_" + self.label] =\
1555  ([self.sigma_dictionary[resolution][0]],
1556  self.sigma_dictionary[resolution][1])
1557  if self.psi_is_sampled:
1558  for psiindex in self.psi_dictionary:
1559  if self.psi_dictionary[psiindex][2]:
1560  ps["Nuisances_ISDCrossLinkMS_Psi_" +
1561  str(psiindex) + "_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
1562  return ps
1563 
1564 #
1565 class CysteineCrossLinkRestraint(object):
1566  def __init__(self, representations, filename, cbeta=False,
1567  betatuple=(0.03, 0.1),
1568  disttuple=(0.0, 25.0, 1000),
1569  omegatuple=(1.0, 1000.0, 50),
1570  sigmatuple=(0.3, 0.3, 1),
1571  betaissampled=True,
1572  sigmaissampled=False,
1573  weightissampled=True,
1574  epsilonissampled=True
1575  ):
1576  # the file must have residue1 chain1 residue2 chain2 fractionvalue epsilonname
1577  # epsilonname is a name for the epsilon particle that must be used for that particular
1578  # residue pair, eg, "Epsilon-Intra-Solvent", or
1579  # "Epsilon-Solvent-Membrane", etc.
1580 
1581  self.representations = representations
1582  self.m = self.representations[0].prot.get_model()
1583  self.rs = IMP.RestraintSet(self.m, 'Cysteine_Crosslink')
1584  self.cbeta = cbeta
1585  self.epsilonmaxtrans = 0.01
1586  self.sigmamaxtrans = 0.1
1587  self.betamaxtrans = 0.01
1588  self.weightmaxtrans = 0.1
1589  self.label = "None"
1590  self.outputlevel = "low"
1591  self.betaissampled = betaissampled
1592  self.sigmaissampled = sigmaissampled
1593  self.weightissampled = weightissampled
1594  self.epsilonissampled = epsilonissampled
1595 
1596  betalower = betatuple[0]
1597  betaupper = betatuple[1]
1598  beta = 0.04
1599  sigma = 10.0
1600  betangrid = 30
1601  crossdataprior = 1
1602 
1603  # beta
1604  self.beta = IMP.pmi.tools.SetupNuisance(
1605  self.m,
1606  beta,
1607  betalower,
1608  betaupper,
1609  betaissampled).get_particle(
1610  )
1611  # sigma
1612  self.sigma = IMP.pmi.tools.SetupNuisance(
1613  self.m,
1614  sigma,
1615  sigmatuple[0],
1616  sigmatuple[1],
1617  sigmaissampled).get_particle()
1618  # population particle
1619  self.weight = IMP.pmi.tools.SetupWeight(
1620  self.m,
1621  weightissampled).get_particle(
1622  )
1623 
1624  # read the file
1625  fl = IMP.pmi.tools.open_file_or_inline_text(filename)
1626  # determine the upperlimit for epsilon
1627  # and also how many epsilon are needed
1628  self.epsilons = {}
1629  data = []
1630  for line in fl:
1631  t = line.split()
1632  if t[0][0] == "#":
1633  continue
1634  if t[5] in self.epsilons:
1635  if 1.0 - float(t[4]) <= self.epsilons[t[5]].get_upper():
1636  self.epsilons[t[5]].set_upper(1.0 - float(t[4]))
1637  else:
1638  self.epsilons[t[5]] = IMP.pmi.tools.SetupNuisance(self.m,
1639  0.01, 0.01, 1.0 - float(t[4]), epsilonissampled).get_particle()
1640  up = self.epsilons[t[5]].get_upper()
1641  low = self.epsilons[t[5]].get_lower()
1642  if up < low:
1643  self.epsilons[t[5]].set_upper(low)
1644 
1645  data.append((int(t[0]), t[1], int(t[2]), t[3], float(t[4]), t[5]))
1646 
1647  # create CrossLinkData
1648  if not self.cbeta:
1649  crossdata = IMP.pmi.tools.get_cross_link_data(
1650  "cysteine", "cysteine_CA_FES.txt.standard",
1651  disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
1652  else:
1653  crossdata = IMP.pmi.tools.get_cross_link_data(
1654  "cysteine", "cysteine_CB_FES.txt.standard",
1655  disttuple, omegatuple, sigmatuple, disttuple[1], disttuple[1], 1)
1656 
1657  # create grids needed by CysteineCrossLinkData
1658  fmod_grid = IMP.pmi.tools.get_grid(0.0, 1.0, 300, True)
1659  omega2_grid = IMP.pmi.tools.get_log_grid(0.001, 10000.0, 100)
1660  beta_grid = IMP.pmi.tools.get_log_grid(betalower, betaupper, betangrid)
1661 
1662  for d in data:
1663  print("--------------")
1664  print("CysteineCrossLink: attempting to create a restraint " + str(d))
1665  resid1 = d[0]
1666  chain1 = d[1]
1667  resid2 = d[2]
1668  chain2 = d[3]
1669  fexp = d[4]
1670  epslabel = d[5]
1671 
1672  # CysteineCrossLinkData
1673 
1675  fexp,
1676  fmod_grid,
1677  omega2_grid,
1678  beta_grid)
1679 
1681  self.m,
1682  self.beta,
1683  self.sigma,
1684  self.epsilons[epslabel],
1685  self.weight,
1686  crossdata,
1687  ccldata)
1688 
1689  failed = False
1690  for i, representation in enumerate(self.representations):
1691 
1692  if not self.cbeta:
1693  p1 = None
1694  p2 = None
1695 
1696  p1 = IMP.pmi.tools.select(representation,
1697  resolution=1, name=chain1,
1698  name_is_ambiguous=False, residue=resid1)[0]
1699 
1700  if p1 is None:
1701  failed = True
1702 
1703  p2 = IMP.pmi.tools.select(representation,
1704  resolution=1, name=chain2,
1705  name_is_ambiguous=False, residue=resid2)[0]
1706 
1707  if p2 is None:
1708  failed = True
1709 
1710  else:
1711  # use cbetas
1712  p1 = []
1713  p2 = []
1714  for t in range(-1, 2):
1715  p = IMP.pmi.tools.select(representation,
1716  resolution=1, name=chain1,
1717  name_is_ambiguous=False, residue=resid1 + t)
1718 
1719  if len(p) == 1:
1720  p1 += p
1721  else:
1722  failed = True
1723  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid1 + t, chain1))
1724 
1725  p = IMP.pmi.tools.select(representation,
1726  resolution=1, name=chain2,
1727  name_is_ambiguous=False, residue=resid2 + t)
1728 
1729  if len(p) == 1:
1730  p2 += p
1731  else:
1732  failed = True
1733  print("\033[93m CysteineCrossLink: missing representation for residue %d of chain %s \033[0m" % (resid2 + t, chain2))
1734 
1735  if not self.cbeta:
1736  if (p1 is not None and p2 is not None):
1737  ccl.add_contribution(p1, p2)
1738  d1 = IMP.core.XYZ(p1)
1739  d2 = IMP.core.XYZ(p2)
1740 
1741  print("Distance_" + str(resid1) + "_" + chain1 + ":" + str(resid2) + "_" + chain2, IMP.core.get_distance(d1, d2))
1742 
1743  else:
1744  if (len(p1) == 3 and len(p2) == 3):
1745  p11n = p1[0].get_name()
1746  p12n = p1[1].get_name()
1747  p13n = p1[2].get_name()
1748  p21n = p2[0].get_name()
1749  p22n = p2[1].get_name()
1750  p23n = p2[2].get_name()
1751 
1752  print("CysteineCrossLink: generating CB cysteine cross-link restraint between")
1753  print("CysteineCrossLink: residue %d of chain %s and residue %d of chain %s" % (resid1, chain1, resid2, chain2))
1754  print("CysteineCrossLink: between particles %s %s %s and %s %s %s" % (p11n, p12n, p13n, p21n, p22n, p23n))
1755 
1756  ccl.add_contribution(p1, p2)
1757 
1758  if not failed:
1759  self.rs.add_restraint(ccl)
1760  ccl.set_name("CysteineCrossLink_" + str(resid1)
1761  + "_" + chain1 + ":" + str(resid2) + "_" + chain2)
1762 
1763  def set_label(self, label):
1764  self.label = label
1765 
1766  def add_to_model(self):
1767  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
1768 
1769  def get_particles_to_sample(self):
1770  ps = {}
1771  if self.epsilonissampled:
1772  for eps in self.epsilons.keys():
1773  ps["Nuisances_CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
1774  self.label] = ([self.epsilons[eps]], self.epsilonmaxtrans)
1775  if self.betaissampled:
1776  ps["Nuisances_CysteineCrossLinkRestraint_beta_" +
1777  self.label] = ([self.beta], self.betamaxtrans)
1778  if self.weightissampled:
1779  ps["Weights_CysteineCrossLinkRestraint_" +
1780  self.label] = ([self.weight], self.weightmaxtrans)
1781  if self.sigmaissampled:
1782  ps["Nuisances_CysteineCrossLinkRestraint_" +
1783  self.label] = ([self.sigma], self.sigmamaxtrans)
1784  return ps
1785 
1786  def set_output_level(self, level="low"):
1787  # this might be "low" or "high"
1788  self.outputlevel = level
1789 
1790  def get_restraint(self):
1791  return self.rs
1792 
1793  def get_sigma(self):
1794  return self.sigma
1795 
1796  def get_output(self):
1797  output = {}
1798  score = self.rs.unprotected_evaluate(None)
1799  output["_TotalScore"] = str(score)
1800  output["CysteineCrossLinkRestraint_Score_" + self.label] = str(score)
1801  output["CysteineCrossLinkRestraint_sigma_" +
1802  self.label] = str(self.sigma.get_scale())
1803  for eps in self.epsilons.keys():
1804  output["CysteineCrossLinkRestraint_epsilon_" + eps + "_" +
1805  self.label] = str(self.epsilons[eps].get_scale())
1806  output["CysteineCrossLinkRestraint_beta_" +
1807  self.label] = str(self.beta.get_scale())
1808  for n in range(self.weight.get_number_of_states()):
1809  output["CysteineCrossLinkRestraint_weights_" +
1810  str(n) + "_" + self.label] = str(self.weight.get_weight(n))
1811 
1812  if self.outputlevel == "high":
1813  for rst in self.rs.get_restraints():
1814  output["CysteineCrossLinkRestraint_Total_Frequency_" +
1815  IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
1816  "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_model_frequency()
1817  output["CysteineCrossLinkRestraint_Standard_Error_" +
1818  IMP.isd.CysteineCrossLinkRestraint.get_from(
1819  rst).get_name(
1820  ) + "_"
1821  + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_standard_error()
1822  if len(self.representations) > 1:
1823  for i in range(len(self.prots)):
1824  output["CysteineCrossLinkRestraint_Frequency_Contribution_" +
1825  IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_name() +
1826  "_State_" + str(i) + "_" + self.label] = IMP.isd.CysteineCrossLinkRestraint.get_from(rst).get_frequencies()[i]
1827 
1828  return output
def set_sigma_is_sampled
Switch on/off the sampling of sigma particles.
Setup cross-link distance restraints from mass spectrometry data.
A 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 ...
def add_to_model
Add the restraint to the model.
Various classes to hold sets of particles.
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
Miscellaneous utilities.
Definition: tools.py:1
Restrain atom pairs based on a set of crosslinks.
def set_weight
Set the weight to apply to all internal restraints.
Add scale parameter to particle.
Definition: Scale.h:24
The type of an atom.
Add uncertainty to a particle.
Definition: Uncertainty.h:24
log
Definition: log.py:1
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:690
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
def get_particles_to_sample
Get the particles to be sampled by the IMP.pmi.sampler object.
def get_restraint_for_rmf
Get the restraint for visualization in an RMF file.
Uniform distribution with harmonic boundaries.
Definition: UniformPrior.h:20
Object used to hold a set of restraints.
Definition: RestraintSet.h:36
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Classes to handle different kinds of restraints.
Warning related to handling of structures.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
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:86
def get_movers
Get all need data to construct a mover in IMP.pmi.dof class.
def get_particle_pairs
Get a list of tuples containing the particle pairs.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
def set_label
Set the unique label used in outputs and particle/restraint names.
int get_state_index(Hierarchy h)
Walk up the hierarchy to find the current state.
Linear function
Definition: Linear.h:19
def plot_violations
Create CMM files, one for each state, of all crosslinks.
Modify a set of continuous variables using a normal distribution.
Definition: NormalMover.h:20
Classes for writing output files and processing them.
Definition: output.py:1
def deprecated_object
Python decorator to mark a class as deprecated.
Definition: __init__.py:9704
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
def create_restraints_for_rmf
create dummy harmonic restraints for each XL but don't add to model Makes it easy to see each contrib...
The general base class for IMP exceptions.
Definition: exception.h:49
Setup cross-link distance restraints at atomic level The "atomic" aspect is that it models the part...
Residue get_residue(Atom d, bool nothrow=false)
Return the Residue containing this atom.
Calculate the -Log of a list of restraints.
Definition: LogWrapper.h:19
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
Definition: VectorD.h:209
Applies a PairScore to a Pair.
Definition: PairRestraint.h:29
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:545
def get_restraint
Get the primary restraint set.
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:337
Functionality for loading, creating, manipulating and scoring atomic structures.
std::string get_chain_id(Hierarchy h)
Walk up the hierarchy to determine the chain id.
def load_nuisances_from_stat_file
Read a stat file and load all the sigmas.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
def get_output
Get outputs to write to stat files.
Base class for PMI restraints, which wrap IMP.Restraint(s).
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A restraint for cysteine cross-linking data.