IMP Reference Guide  develop.27926d84dc,2024/04/20 The Integrative Modeling Platform
nonmantained.py
1 """@namespace IMP.pmi1.nonmaintained
2  Nonmaintained code.
3 """
4
5 from __future__ import print_function
6 import IMP
7 import IMP.atom
8
9 IMP.deprecated_module('2.8', __name__,
10  "If you use any functionality from this module, please let the "
11  "IMP developers know. It will be removed in the next IMP release.")
12
13 class Rods(object):
14
15  def __init__(self, m):
16  self.m = m
18  self.rigid_bodies = []
19  self.floppy_bodies = []
20  self.maxtrans_rb = 2.0
21  self.maxtrans_fb = 2.0
22  self.maxrot_rb = 0.15
23
24  def add_protein(self, name, res):
25  (firstres, lastres) = res
26  from math import pi, cos, sin
28  h.set_name(name)
29  nres = lastres - firstres
30  radius = (nres) * 5 / 2 / pi
31
32  for res in range(firstres, lastres):
33  alpha = 2 * pi / nres * (res - firstres)
34  x = radius * cos(alpha)
35  y = radius * sin(alpha)
36  p = IMP.Particle(self.m)
37  r = IMP.atom.Residue.setup_particle(p, IMP.atom.ALA, res)
39  d.set_coordinates(IMP.algebra.Vector3D((x, y, 0)))
40  d.set_coordinates_are_optimized(True)
43
44  def get_hierarchy(self):
45  return self.hier
46
47  def set_rod(self, chainname, res):
48  (firstres, lastres) = res
49  prb = IMP.Particle(self.m)
50  sel = IMP.atom.Selection(
51  self.hier,
52  molecule=chainname,
53  residue_indexes=list(range(
54  firstres,
55  lastres + 1)))
56  ps = sel.get_selected_particles()
58  self.rigid_bodies.append(rb)
59
60  def get_particles_to_sample(self):
61  # get the list of samplable particles with their type
62  # and the mover displacement. Everything wrapped in a dictionary,
63  # to be used by samplers modules
64  ps = {}
65  ps["Floppy_Bodies_Rods"] = (self.floppy_bodies, self.maxtrans_fb)
66  ps["Rigid_Bodies_Rods"] = (
67  self.rigid_bodies,
68  self.maxtrans_rb,
69  self.maxrot_rb)
70  return ps
71
72
74
75  def __init__(self, m):
76  self.m = m
77  self.beads = []
78  self.nresidues = 0
80  self.floppy_bodies = []
81  self.maxtrans_fb = 0.2
82  self.particle_database = {}
83
85
86  p = IMP.Particle(self.m)
87  p.set_name(label)
88  self.particle_database[label] = p
89  self.floppy_bodies.append(p)
90  # set default coordinates 0,0,0
92  d = IMP.core.XYZR.setup_particle(p, radius)
93  d.set_coordinates_are_optimized(True)
94  # a=IMP.atom.Atom.setup_particle(p,IMP.atom.AT_CA)
95  # p=IMP.Particle(self.m)
96  self.nresidues += 1
97  # r=IMP.atom.Residue.setup_particle(p,IMP.atom.ALA,self.nresidues)
101  if not color is None:
102  self.set_color(label, color)
103  return self.particle_database[label]
104
105  def set_color(self, label, value):
106  p = self.particle_database[label]
107  clr = IMP.display.get_rgb_color(value)
109
110  def set_floppy_bodies_max_trans(self, maxtrans):
111  self.maxtrans_fb = maxtrans
112
113  def get_hierarchy(self):
114  return self.hier
115
116  def get_bead(self, label):
117  return self.particle_database[label]
118
119  def set_maxtrans_fb(self, maxtrans_fb):
120  self.maxtrans_fb = maxtrans_fb
121
122  def get_particles_to_sample(self):
123  # get the list of samplable particles with their type
124  # and the mover displacement. Everything wrapped in a dictionary,
125  # to be used by samplers modules
126  ps = {}
127  ps["Floppy_Bodies_Beads"] = (self.floppy_bodies, self.maxtrans_fb)
128  return ps
129
130
131 class MultipleStates(object):
132
133  def __init__(self, nstates, m):
134  global itertools, tools, restraints
135
136  import itertools
137  import IMP.pmi1.tools as tools
138  import IMP.pmi1.restraints as restraints
139
140  self.floppy_bodies = []
141  self.rigid_bodies = []
142  for ncopy in range(nstates):
143  self.floppy_bodies.append([])
144  self.rigid_bodies.append([])
145
146  self.rigid_bodies_are_sampled = True
147  self.floppy_bodies_are_sampled = True
148  self.prot = []
149  self.refprot = []
150  self.prot_lowres = {}
151  self.nstates = nstates
152  self.label = "None"
153
154  # model decorator list
155  self.xyzmodellist = []
156  self.xyzreflist = []
157  self.maxtrans_rb = 0.15
158  self.maxrot_rb = 0.03
159  self.maxtrans_fb = 0.15
160  self.m = m
161
162  def get_model(self):
163  return self.m
164
165  def set_label(self, label):
166  self.label = label
167
168  def set_rigid_bodies_are_sampled(self, input=True):
169  self.rigid_bodies_are_sampled = input
170
171  def set_floppy_bodies_are_sampled(self, input=True):
172  self.floppy_bodies_are_sampled = input
173
174  def get_rigid_bodies(self):
175  return self.rigid_bodies
176
177  def set_rigid_bodies_max_trans(self, maxtrans):
178  self.maxtrans_rb = maxtrans
179
180  def set_rigid_bodies_max_rot(self, maxrot):
181  self.maxrot_rb = maxrot
182
183  def set_floppy_bodies_max_trans(self, maxtrans):
184  self.maxtrans_fb = maxtrans
185
186  def get_hierarchies(self):
187  return self.prot
188
189  def destroy_residues(self, segments):
190  # segments are defined as a list of tuples ex [(res1,res2,chain),....]
191  # this function must be called before the rigid body definition!
192  for prot in self.prot:
193  for segment in segments:
194  # rinterval=[(segment[0],segment[1]+1)]
195  if (segment[0] == -1 or segment[1] == -1):
196  s = IMP.atom.Selection(prot, chains=segment[2])
197  else:
198  s = IMP.atom.Selection(
199  prot,
200  chains=segment[2],
201  residue_indexes=list(range(segment[0],
202  segment[1] + 1)))
203  for p in s.get_selected_particles():
205  print("MultipleStates: one particle was not destroyed because it was a RigidMember.")
206  else:
207  # destroy the residue and the associated atom
208  a = IMP.atom.Atom(p)
209  r = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
210  IMP.atom.destroy(r)
211  # IMP.atom.destroy(a)
212  # IMP.atom.destroy(p)
214
216  self,
217  residuechainlist,
218  residue_type=IMP.atom.LYS):
219  # add a list of residues to the corresponding list
220  # for instance residuechainlist=[(35,"A"),(100,"B")] will add
221  # residue 35 to chain A and residue 100 to chain B
222  for rc in residuechainlist:
223
224  s = IMP.atom.Selection(
225  self.prot[0],
226  chains=rc[1],
227  residue_index=rc[0],
228  atom_type=IMP.atom.AT_CA)
229
230  print(s.get_selected_particles())
231  if len(s.get_selected_particles()) == 0:
232  for prot in self.prot:
233  print("adding " + str(rc))
234  p = IMP.Particle(self.m)
235  # set default coordinates 0,0,0
238  d.set_coordinates_are_optimized(True)
239  a = IMP.atom.Atom.setup_particle(p, IMP.atom.AT_CA)
240  p = IMP.Particle(self.m)
241  r = IMP.atom.Residue.setup_particle(p, residue_type, rc[0])
243  p = IMP.Particle(self.m)
244  c = IMP.atom.Chain.setup_particle(p, rc[1])
247  print(tools.get_residue_index_and_chain_from_particle(a))
248
249  else:
250  p = s.get_selected_particles()[0]
251  # , tools.get_residue_index_and_chain_from_particle(s.get_selected_particles()[0])
252  print(rc, s.get_selected_particles()[0])
253
254  # test that that was indeed added:
255  s = IMP.atom.Selection(
256  self.prot[0],
257  chains=rc[1],
258  residue_index=rc[0],
259  atom_type=IMP.atom.AT_CA)
260
261  print(s.get_selected_particles())
262
264  '''
265  this method generate beads in missing portions.
266  The segments argument must be a list of selections
267  in the form [(firstres,lastres,chain)]
268  each selection will generate a bead
269  '''
270  if xyzs is None:
271  xyzs = []
272  if radii is None:
273  radii = []
274  if colors is None:
275  colors = []
276
277  from math import pi
278
279  for n, s in enumerate(segments):
280  firstres = s[0]
281  lastres = s[1]
282  chainid = s[2]
283  nres = s[1] - s[0]
284  for prot in self.prot:
285  for prot in self.prot:
286  cps = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
287  for c in cps:
288  chain = IMP.atom.Chain(c)
289  if chain.get_id() == chainid:
290  p = IMP.Particle(self.m)
292  rindexes = list(range(firstres, lasteres + 1))
293  f.set_residue_indexes(rindexes)
294  f.set_name(
295  "Fragment_" + '%i-%i' %
296  (firstres, lastres))
298  mass = len(rindexes) * 110.0
300  if n + 1 > len(radii):
301  mass = len(rindexes) * 110.0
303  radius = (3 * vol / math.pi) ** (1 / 3)
304  else:
306
307  if n + 1 > len(xyzs):
308  x = 0
309  y = 0
310  z = 0
311  else:
312  x = xyzs[n][0]
313  y = xyzs[n][1]
314  z = xyzs[n][2]
315
316  if n + 1 <= len(colors):
317  clr = IMP.display.get_rgb_color(colors[n])
319
320  d = IMP.atom.XYZR.setup_particle(
321  p, IMP.algebra.Sphere3D(x, y, z, radius))
322
323  def renumber_residues(self, chainid, newfirstresiduenumber):
324  for prot in self.prot:
325  cps = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
326  for c in cps:
327  if IMP.atom.Chain(c).get_id() == chainid:
328  ps = c.get_children()
329  r = IMP.atom.Residue(ps[0])
330  ri = r.get_index()
331  offs = newfirstresiduenumber - ri
332  for p in ps:
333  r = IMP.atom.Residue(p)
334  ri = r.get_index()
335  r.set_index(ri + offs)
336
337  def destroy_everything_but_the_residues(self, segments):
338  # segments are defined as a list of tuples ex [(res1,res2,chain),....]
339  for prot in self.prot:
340  pstokeep = []
341  for segment in segments:
342
343  # rinterval=[(segment[0],segment[1]+1)]
344  if (segment[0] == -1 or segment[1] == -1):
345  s = IMP.atom.Selection(prot, chains=segment[2])
346  else:
347  s = IMP.atom.Selection(
348  prot,
349  chains=segment[2],
350  residue_indexes=list(range(segment[0],
351  segment[1] + 1)))
352  pstokeep += s.get_selected_particles()
353
354  for p in IMP.atom.get_leaves(prot):
355  if p not in pstokeep:
357  print("MultipleStates: one particle was not destroyed because it was a RigidMember.")
358  else:
359  # destroy the residue and the associated atom
360  a = IMP.atom.Atom(p).get_parent()
361  r = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
362  # IMP.atom.destroy(a)
363  IMP.atom.destroy(r)
364  # self.m.remove_particle(p)
365
366  def generate_linkers_restraint_and_floppy_bodies(self, segment):
367  '''
368  this methods automatically links the particles consecutively
369  according to the sequence. The restraint applied is a harmonic upper bound,
370  with a distance that is proportional to the number of residues
371  in the gap.
372  '''
373  # this function will create floppy bodies where there are not
374  # rigid bodies and moreover create a linker restraint between them
375  linker_restraint_objects = []
376  for ncopy, prot in enumerate(self.prot):
377  if (segment[0] == -1 or segment[1] == -1):
378  s = IMP.atom.Selection(prot, chains=segment[2])
379  else:
380  s = IMP.atom.Selection(
381  prot,
382  chains=segment[2],
383  residue_indexes=list(range(segment[0],
384  segment[1] + 1)))
385  residue_indexes = []
386  for p in s.get_selected_particles():
387
388  (r, c) = tools.get_residue_index_and_chain_from_particle(p)
389
391  Floppy = False
392  else:
393  (r, c) = tools.get_residue_index_and_chain_from_particle(p)
394  p.set_name(str(r) + ":" + str(c))
395  tools.set_floppy_body(p)
396  self.floppy_bodies[ncopy].append(p)
397  Floppy = True
398  residue_indexes.append((r, Floppy, c, p))
399
400  residue_indexes.sort()
401
402  pruned_residue_list = []
403  r0 = residue_indexes[0]
404  pruned_residue_list.append(r0)
405
406  # generate the list of residues that define the intervals
407  # between rigid bodies and floppy bodies
408  for i in range(1, len(residue_indexes)):
409  r = residue_indexes[i]
410  if r[1] == r0[1] and r[1] == False and IMP.core.RigidMember(r[3]).get_rigid_body() == IMP.core.RigidMember(r0[3]).get_rigid_body():
411  r0 = r
412  elif r[1] == r0[1] and r[1] == False and IMP.core.RigidMember(r[3]).get_rigid_body() != IMP.core.RigidMember(r0[3]).get_rigid_body():
413  pruned_residue_list.append(r0)
414  pruned_residue_list.append(r)
415  r0 = r
416  elif r[1] != r0[1] and r0[1] == False:
417  pruned_residue_list.append(r0)
418  pruned_residue_list.append(r)
419  r0 = r
420  elif r[1] == r0[1] and r0[1]:
421  pruned_residue_list.append(r)
422  r0 = r
423  elif r[1] != r0[1] and r[1] == False:
424  pruned_residue_list.append(r)
425  r0 = r
426
427  r0 = pruned_residue_list[0]
428  linkdomaindef = []
429
430  for i in range(1, len(pruned_residue_list)):
431  r = pruned_residue_list[i]
432  if r[1] == r0[1] and r[1] == False and IMP.core.RigidMember(r[3]).get_rigid_body() == IMP.core.RigidMember(r0[3]).get_rigid_body():
433  r0 = r
434  else:
435  linkdomaindef.append((r0[0], r[0], r[2]))
436  r0 = r
437
438  print(" creating linker between atoms defined by: " + str(linkdomaindef))
439
440  ld = restraints.LinkDomains(prot, linkdomaindef, 1.0, 3.0)
441  ld.set_label(str(ncopy))
444  prs = ld.get_pairs()
445
447
448  def get_ref_hierarchies(self):
449  return self.refprot
450
451  def get_number_of_states(self):
452  return self.nstates
453
454  def get_rigid_bodies(self):
455  rblist = []
456  for rbl in self.rigid_bodies:
457  for rb in rbl:
458  rblist.append(rb)
459  return rblist
460
461  def get_floppy_bodies(self):
462  fblist = []
463  for fbl in self.floppy_bodies:
464  for fb in fbl:
465  fblist.append(fb)
466  return fblist
467
468  def set_rigid_bodies(self, rigid_body_list):
469  if len(self.prot) == 0:
470  print("MultipleStates.set_rigid_bodies: hierarchy was not initialized")
471  exit()
472  for ncopy, prot in enumerate(self.prot):
473  rbl = []
474  for element in rigid_body_list:
475  atoms = []
476  for interval in element:
477  # rinterval upper bound is incremented by one because the
478  # residue_indexes attribute cuts the upper edge
479  # rinterval=[(interval[0],interval[1]+1)]
480  if (interval[0] == -1 or interval[1] == -1):
481  s = IMP.atom.Selection(prot, chains=interval[2])
482  else:
483  s = IMP.atom.Selection(
484  prot,
485  chains=interval[2],
486  residue_indexes=list(range(interval[0],
487  interval[1] + 1)))
488  for p in s.get_selected_particles():
489  atoms.append(IMP.core.XYZR(p))
490
491  # add low resolution representation to the rigid bodies
492  for key in self.prot_lowres:
493  if (interval[0] == -1 or interval[1] == -1):
494  s = IMP.atom.Selection(
495  self.prot_lowres[key][ncopy],
496  chains=interval[2])
497  else:
498  s = IMP.atom.Selection(
499  self.prot_lowres[key][
500  ncopy], chains=interval[2],
501  residue_indexes=list(range(interval[0], interval[1] + 1)))
502  for p in s.get_selected_particles():
503  atoms.append(IMP.core.XYZR(p))
504
505  if len(atoms) > 0:
506  prb = IMP.Particle(self.m)
507  rb = IMP.core.RigidBody.setup_particle(prb, atoms)
508  rb.set_name(str(element))
509  rbl.append(rb)
510  else:
511  print("MultipleStates.set_rigid_bodies: selection " + str(interval) + " has zero elements")
512  self.rigid_bodies[ncopy] += rbl
513
514  def set_floppy_bodies(self, floppy_body_list):
515  # define flexible regions within rigid bodies
516
517  if len(self.prot) == 0:
518  print("MultipleStates: hierarchy was not initialized")
519  exit()
520
521  for ncopy, prot in enumerate(self.prot):
522  atoms = []
523  for element in floppy_body_list:
524
525  for interval in element:
526  # rinterval upper bound is incremented by one because the
527  # residue_indexes attribute cuts the upper edge
528  # rinterval=[(interval[0],interval[1]+1)]
529  if (interval[0] == -1 or interval[1] == -1):
530  s = IMP.atom.Selection(prot, chains=interval[2])
531  else:
532  s = IMP.atom.Selection(
533  prot,
534  chains=interval[2],
535  residue_indexes=list(range(interval[0],
536  interval[1] + 1)))
537  for p in s.get_selected_particles():
538  (r, c) = tools.get_residue_index_and_chain_from_particle(
539  p)
540  tools.set_floppy_body(p)
541  p.set_name(str(r) + ":" + str(c))
542  atoms.append(IMP.core.XYZR(p))
543  self.floppy_bodies[ncopy] += atoms
544
545  def get_particles_to_sample(self):
546  # get the list of samplable particles with their type
547  # and the mover displacement. Everything wrapped in a dictionary,
548  # to be used by samplers modules
549  ps = {}
550  rblist = self.get_rigid_bodies()
551  fblist = self.get_floppy_bodies()
552  if self.rigid_bodies_are_sampled:
553  ps["Rigid_Bodies_MultipleStates"] = (
554  rblist,
555  self.maxtrans_rb,
556  self.maxrot_rb)
557  if self.floppy_bodies_are_sampled:
558  ps["Floppy_Bodies_MultipleStates"] = (fblist, self.maxtrans_fb)
559  return ps
560
561  def set_hierarchy_from_pdb(self, pdblistoflist):
562  "the input is a list of list of pdbs"
563  "one list for each copy"
564  # eg [["pdb1_copy0","pdb2_copy0"],["pdb1_copy1","pdb2_copy1"]]"
565  for copy in range(0, self.nstates):
566  prot = self.read_pdbs(pdblistoflist[copy])
567  self.prot.append(prot)
569  self.xyzmodellist.append(xyz)
570
571  def set_ref_hierarchy_from_pdb(self, pdblistoflist):
572  "the input is a list of list of pdbs"
573  "one list for each copy"
574  # eg [["pdb1_copy0","pdb2_copy0"],["pdb1_copy1","pdb2_copy1"]]"
575  for copy in range(0, self.nstates):
576  prot = self.read_pdbs(pdblistoflist[copy])
577  self.refprot.append(prot)
579  self.xyzreflist.append(xyz)
580
581  def read_pdbs(self, list_pdb_file):
582  """read pdbs from an external list file
583  create a simplified representation
584  if the pdbs are given a individual strings, it will read the
585  pdbs and give the chain name as specified in the pdb
586  If it s a tuple like (filename,chainname) it will read
587  the pdb and assign a name chainname
588  to the chain"""
589
590  # create an empty hierarchy
592  for pdb in list_pdb_file:
593  if type(pdb) == str:
594  h = IMP.atom.read_pdb(
595  pdb, self.m, IMP.atom.AndPDBSelector(
598
599  '''
600  #destroy CA atoms, for the future
601  for p in IMP.atom.get_leaves(h):
602  coor=IMP.core.XYZ(p).get_coordinates()
603  r=IMP.atom.Hierarchy(p).get_parent()
604  IMP.core.XYZ.setup_particle(r,coor)
605  IMP.atom.destroy(p)
606  '''
607
608  cps = IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE)
609
610  '''
611  #consolidate the chains
612  for c in cps:
613  cid=c.get_id()
614  s0=IMP.atom.Selection(hier, chains=cid)
615  try:
616  p=s0.get_selected_particles()[0]
617  re=IMP.atom.Residue(IMP.atom.Atom(p).get_parent()
618  ch=IMP.atom.Chain(re).get_parent())
619
620  except:
621  continue
622  '''
623
625  if type(pdb) == tuple:
626  h = IMP.atom.read_pdb(
627  pdb[0], self.m, IMP.atom.AndPDBSelector(
630
631  '''
632  #destroy CA atoms, for the future
633  for p in IMP.atom.get_leaves(h):
634  coor=IMP.core.XYZ(p).get_coordinates()
635  r=IMP.atom.Hierarchy(p).get_parent()
636  IMP.core.XYZ.setup_particle(r,coor)
637  IMP.atom.destroy(p)
638  '''
639
640  cps = IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE)
641  for cp in cps:
642  IMP.atom.Chain(cp).set_id(pdb[1])
644
645  return hier
646
647  def recenter(self, prot):
648  "recenter the hierarchy"
649  ps = IMP.atom.get_leaves(prot)
650  center = IMP.algebra.get_zero_vector_3d()
651  for l in ps:
652  center += IMP.core.XYZ(l).get_coordinates()
653  center /= len(ps)
654  for l in ps:
655  d = IMP.core.XYZ(l)
656  d.set_coordinates(d.get_coordinates() - center)
657  d.set_coordinates_are_optimized(True)
658
659  '''
660  # bug generating code: keeping it for history
661
662  rb=IMP.atom.create_rigid_body(prot)
663  rbcoord=rb.get_coordinates()
664  rot=IMP.algebra.get_identity_rotation_3d()
665  tmptrans=IMP.algebra.Transformation3D(rot,rbcoord)
666  trans=tmptrans.get_inverse()
667  IMP.core.transform(rb,trans)
668  IMP.core.RigidBody.teardown_particle(rb)
669  self.m.remove_particle(rb)
670  '''
671
672  def shuffle_configuration(self, bounding_box_length):
673  "shuffle configuration, used to restart the optimization"
674  "it only works if rigid bodies were initialized"
675  if len(self.rigid_bodies) == 0:
676  print("MultipleStates: rigid bodies were not initialized")
677  hbbl = bounding_box_length / 2
678  for rbl in self.rigid_bodies:
679  for rb in rbl:
680  ub = IMP.algebra.Vector3D(-hbbl, -hbbl, -hbbl)
681  lb = IMP.algebra.Vector3D(hbbl, hbbl, hbbl)
682  bb = IMP.algebra.BoundingBox3D(ub, lb)
683  translation = IMP.algebra.get_random_vector_in(bb)
685  transformation = IMP.algebra.Transformation3D(
686  rotation,
687  translation)
688  rb.set_reference_frame(
689  IMP.algebra.ReferenceFrame3D(transformation))
690
691  def generate_simplified_hierarchy(self, nres):
692  # generate a new multistate hierarchy
693  self.prot_lowres[nres] = []
694  for prot in self.prot:
695  sh = IMP.atom.create_simplified_along_backbone(prot, nres, False)
696  print(IMP.atom.get_leaves(sh))
697  # for p in IMP.atom.get_leaves(sh):
698  # IMP.atom.Atom.setup_particle(p,IMP.atom.AT_CA)
699  # s=IMP.atom.Selection(sh, chains="A",
700  # residue_index=958)
701  # print s.get_selected_particles()[0]
702  self.prot_lowres[nres].append(sh)
703
704  def get_simplified_hierarchy(self, nres):
705  return self.prot_lowres[nres]
706
707  def calculate_drms(self):
708  # calculate DRMSD matrix
709
710  if len(self.xyzmodellist) == 0:
711  print("MultipleStates: hierarchies were not initialized")
712
713  if len(self.xyzreflist) == 0:
714  print("MultipleStates: reference hierarchies were not initialized")
715
716  drmsd = {}
717  for i in range(len(self.xyzreflist)):
718  for j in range(len(self.xyzmodellist)):
719  try:
720  drmsdval = IMP.atom.get_drmsd(
721  self.xyzmodellist[j],
722  self.xyzreflist[i])
723  except:
724  drmsdval = tools.get_drmsd(
725  self.xyzmodellist[j],
726  self.xyzreflist[i])
727  drmsd["MultipleStates_DRMSD_" +
728  str(i) + "-Model_" + str(j)] = drmsdval
729
730  # calculate model-template assignment that gives minimum total drmsd
731  min_drmsd = []
733  s = 0.
734  for i, j in enumerate(assign):
735  s += drmsd["MultipleStates_DRMSD_" +
736  str(j) + "-Model_" + str(i)]
737  min_drmsd.append(s)
738
739  drmsd["MultipleStates_Total_DRMSD"] = min(min_drmsd)
740  return drmsd
741
742  def get_output(self):
743  output = {}
744  if len(self.refprot) != 0:
745  drms = self.calculate_drms()
746  output.update(drms)
747  output["MultipleStates_Total_Score_" +
748  self.label] = str(self.m.evaluate(False))
749  return output
750
751
752 #
753
755
756  def __init__(self, prot, resrangelist, kappa, length=5.0):
757  # generate a linker between residues using HarmonicUpperBound
758  # restraints. Define a list of linked residues,
759  # e.g. [(3,5,"A"),(9,10,"B")]
760  # will link residues 3 and 5 belonging to chain A and
761  # residues 9 and 10 belonging to chain B
762  self.rs = IMP.RestraintSet('linker')
763  self.prot = prot
764  self.kappa = kappa
765  self.resrangelist = resrangelist
766  self.length = length
767  self.label = "None"
768
769  self.m = self.prot.get_model()
770  self.pairs = []
771
772  for pair in self.resrangelist:
773  c0 = pair[2]
774  r0 = pair[0]
775  c1 = pair[2]
776  r1 = pair[1]
777  try:
778  s0 = IMP.atom.Selection(
779  self.prot,
780  chains=c0,
781  residue_index=r0,
782  atom_type=IMP.atom.AT_CA)
783  p0 = s0.get_selected_particles()[0]
784  except:
786  continue
787  try:
788  s1 = IMP.atom.Selection(
789  self.prot,
790  chains=c1,
791  residue_index=r1,
792  atom_type=IMP.atom.AT_CA)
793  p1 = s1.get_selected_particles()[0]
794  except:
796  continue
797  # check this is the residue length (should be 4, but I use a larger
798  # length)
799  dist0 = float(pair[1] - pair[0]) * self.length
800  h = IMP.core.HarmonicUpperBound(dist0, self.kappa)
802  pr = IMP.core.PairRestraint(dps, IMP.ParticlePair(p0, p1))
803  pr.set_name(
804  "LinkDomains_" + str(pair[0]) + "-" + str(pair[1]) + "_" + str(pair[2]))
806  self.pairs.append((p0, p1, r0, c0, r1, c1, pr))
807
808  def set_label(self, label):
809  self.label = label
810
813
814  def get_pairs(self):
815  return self.pairs
816
817  def get_hierarchy(self):
818  return self.prot
819
820  def get_kappa(self):
821  return self.kappa
822
823  def get_residue_ranges(self):
824  return self.resrangelist
825
826  def get_length(self):
827  return self.length
828
829  def get_restraint(self):
830  return self.rs
831
832  def get_restraints(self):
833  rlist = []
834  for r in self.rs.get_restraints():
835  rlist.append(IMP.core.PairRestraint.get_from(r))
836  return rlist
837
838  def print_chimera_pseudobonds(self, filesuffix, model=0):
839  f = open(filesuffix + ".chimera", "w")
840  atype = "ca"
841  for p in self.pairs:
842  s = "#" + str(model) + ":" + str(p[2]) + "." + p[3] + "@" + atype + \
843  " #" + str(model) + ":" + str(p[4]) + "." + p[5] + "@" + atype
844  f.write(s + "\n")
845  f.close()
846
847  def get_output(self):
848  self.m.update()
849  output = {}
850  score = self.rs.unprotected_evaluate(None)
851  output["_TotalScore"] = str(score)
852  output["LinkDomains_" + self.label] = str(score)
853  for rst in self.rs.get_restraints():
855  # IMP.core.PairRestraint.get_from(rst).get_name()+
856  # "_"+self.label]=IMP.core.PairRestraint.get_from(rst).evaluate(False)
857  output["LinkDomains_" + rst.get_name() +
858  "_" + self.label] = rst.unprotected_evaluate(None)
859
860  for i in range(len(self.pairs)):
861
862  p0 = self.pairs[i][0]
863  p1 = self.pairs[i][1]
864  r0 = self.pairs[i][2]
865  c0 = self.pairs[i][3]
866  r1 = self.pairs[i][4]
867  c1 = self.pairs[i][5]
868
869  label = str(r0) + ":" + c0 + "_" + str(r1) + ":" + c1
870
871  d0 = IMP.core.XYZ(p0)
872  d1 = IMP.core.XYZ(p1)
873  output["LinkDomains_Distance_" + label + "_" +
874  self.label] = str(IMP.core.get_distance(d0, d1))
875
876  return output
877
878
879 #
880
881
882 class UpperBound(object):
883
884  def __init__(self, prot, respairs, kappa, length=5.0):
885  # generate a restraint between residues using HarmonicUpperBound
886  # restraints. Define a list of linked residues,
887  # e.g. [(3,"A",5,"B"),(9,"B",10,"C")]
888  # will link residues 3 and 5 belonging to chain A and B and
889  # residues 9 and 10 belonging to chain B and C
890  self.rs = IMP.RestraintSet('upperbound')
891
892  self.prot = prot
893  self.kappa = kappa
894  self.respairs = respairs
895  self.length = length
896  self.label = "None"
897
898  self.m = self.prot.get_model()
899
900  for pair in self.respairs:
901  try:
902  s0 = IMP.atom.Selection(
903  self.prot,
904  chains=pair[1],
905  residue_index=pair[0],
906  atom_type=IMP.atom.AT_CA)
907  p0 = s0.get_selected_particles()[0]
908  except:
909  "UpperBound: error"
910  continue
911  try:
912  s1 = IMP.atom.Selection(
913  self.prot,
914  chains=pair[3],
915  residue_index=pair[2],
916  atom_type=IMP.atom.AT_CA)
917  p1 = s1.get_selected_particles()[0]
918  except:
919  "UpperBound: error"
920  continue
921  # check this is the residue length (should be 4, but I use a larger
922  # length)
923
924  h = IMP.core.HarmonicUpperBound(self.length, self.kappa)
926  pr = IMP.core.PairRestraint(dps, IMP.ParticlePair(p0, p1))
927  pr.set_name(
928  "UpperBound_" + str(pair[0]) + "-" + str(pair[1]) + "_" + str(pair[2]))
930
931  def set_label(self, label):
932  self.label = label
933
936
937  def get_hierarchy(self):
938  return self.prot
939
940  def get_kappa(self):
941  return self.kappa
942
943  def get_residue_ranges(self):
944  return self.respairs
945
946  def get_length(self):
947  return self.length
948
949  def get_restraint(self):
950  return self.rs
951
952  def get_output(self):
953  output = {}
954  self.m.update()
955  score = self.rs.unprotected_evaluate(None)
956  output["_TotalScore"] = str(score)
957  output["UpperBound_" + self.label] = str(score)
958  return output
959
960 #
961
962
963 class ExcludedVolumeResidue(object):
964
965  def __init__(self, prot, kappa):
966  self.rs = IMP.RestraintSet('excluded_volume')
967  self.prot = prot
968  self.kappa = kappa
969  self.label = "None"
970  self.m = self.prot.get_model()
971
972  atoms = IMP.atom.get_by_type(self.prot, IMP.atom.ATOM_TYPE)
973  for atom in atoms:
974  restype = IMP.atom.Residue(
975  IMP.atom.Atom(atom).get_parent()).get_residue_type()
981
982  evr = IMP.core.ExcludedVolumeRestraint(lsa, self.kappa)
984
985  def set_label(self, label):
986  self.label = label
987
988  def add_excluded_particle_pairs(self, excluded_particle_pairs):
989  # add pairs to be filtered when calculating the score
990  lpc = IMP.container.ListPairContainer(self.m)
993  IMP.core.ExcludedVolumeRestraint.get_from(
995
998
999  def get_hierarchy(self):
1000  return self.prot
1001
1002  def get_kappa(self):
1003  return self.kappa
1004
1005  def get_restraint(self):
1006  return self.rs
1007
1008  def get_output(self):
1009  self.m.update()
1010  output = {}
1011  score = self.rs.unprotected_evaluate(None)
1012  output["_TotalScore"] = str(score)
1013  output["ExcludedVolumeResidue_" + self.label] = str(score)
1014  return output
1015
1016 #
1017
1018
1019 class BipartiteExcludedVolumeResidue(object):
1020
1021  def __init__(self, prot1, prot2, kappa):
1022  self.rs = IMP.RestraintSet('bipartite_excluded_volume')
1023  self.prot1 = prot1
1024  self.prot2 = prot2
1025  self.kappa = kappa
1026  self.label = "None"
1027  self.m = self.prot.get_model()
1028
1029  atoms1 = IMP.atom.get_by_type(prot1, IMP.atom.ATOM_TYPE)
1032  for atom in atoms1:
1033  restype = IMP.atom.Residue(
1034  IMP.atom.Atom(atom).get_parent()).get_residue_type()
1038
1039  atoms2 = IMP.atom.get_by_type(prot2, IMP.atom.ATOM_TYPE)
1042  for atom in atoms2:
1043  restype = IMP.atom.Residue(
1044  IMP.atom.Atom(atom).get_parent()).get_residue_type()
1048
1050  ls_ref,
1051  ls_symm,
1052  kappa,
1053  10.0)
1054  ssps = IMP.core.SoftSpherePairScore(kappa)
1055  evr3 = IMP.container.PairsRestraint(ssps, cbpc)
1058
1059  def set_label(self, label):
1060  self.label = label
1061
1064
1065  def get_hierarchy(self):
1066  return self.prot1, self.prot2
1067
1068  def get_kappa(self):
1069  return self.kappa
1070
1071  def get_restraint(self):
1072  return self.rs
1073
1074  def get_output(self):
1075  self.m.update()
1076  output = {}
1077  score = self.rs.unprotected_evaluate(None)
1078  output["_TotalScore"] = str(score)
1079  output["BipartiteExcludedVolumeResidue_" + self.label] = str(score)
1080  return output
1081
1082 #
1083
1084
1085 class TemplateRestraint(object):
1086
1087  def __init__(self, ps1, ps2, cutoff=6.5, kappa=1.0, forcerb=False):
1088  self.m = ps1[0].get_model()
1089  self.label = "None"
1090  self.cutoff = cutoff
1091  self.kappa = kappa
1092  # this parameter ovverides the rigid body filter below
1093  self.forcerb = forcerb
1094  self.rs = IMP.RestraintSet('template_restraint')
1095  for p1 in ps1:
1096  for p2 in ps2:
1097  # check that the two particles are not in the same rigid body
1099  IMP.core.RigidMember(p1).get_rigid_body() == IMP.core.RigidMember(p2).get_rigid_body()) and not self.forcerb:
1100  continue
1101  d0 = IMP.core.XYZ(p1)
1102  d1 = IMP.core.XYZ(p2)
1103  dist = IMP.core.get_distance(d0, d1)
1104  if dist <= self.cutoff:
1105  hf = IMP.core.Harmonic(dist, self.kappa)
1106  dps = IMP.core.DistancePairScore(hf)
1107  pr = IMP.core.PairRestraint(dps, IMP.ParticlePair(p1, p2))
1109
1110  def set_label(self, label):
1111  self.label = label
1112
1115
1116  def get_cutoff(self):
1117  return self.cutoff
1118
1119  def get_kappa(self):
1120  return self.kappa
1121
1122  def get_restraint(self):
1123  return self.rs
1124
1125  def get_output(self):
1126  self.m.update()
1127  output = {}
1128  score = self.rs.unprotected_evaluate(None)
1129  output["_TotalScore"] = str(score)
1130  output["TemplateRestraint_" + self.label] = str(score)
1131  return output
1132
1133 #
1134
1135
1136 class MarginalChi3Restraint(object):
1137
1138  def __init__(self, part1, part2):
1139  global impisd2, tools
1140  import IMP.isd2 as impisd2
1141  import IMP.pmi1.tools as tools
1142
1143  self.m = part1.get_model()
1144  self.label = "None"
1145  self.rs = IMP.RestraintSet('chi3_restraint')
1146  self.sigmamaxtrans = 0.1
1147
1148  self.ps1 = IMP.atom.get_leaves(part1)
1149  self.ps2 = IMP.atom.get_leaves(part2)
1150  self.sigma = tools.SetupNuisance(
1151  self.m,
1152  1.0,
1153  0.1,
1154  100.0,
1155  True).get_particle(
1156  )
1157
1158  for i in range(len(self.ps1)):
1159  mc = impisd2.MarginalChi3Restraint(
1160  self.ps1[i],
1161  self.ps2[i],
1162  self.sigma)
1164
1165  def set_label(self, label):
1166  self.label = label
1167
1170
1171  def get_restraint(self):
1172  return self.rs
1173
1174  def get_particles_to_sample(self):
1175  ps = {}
1176  ps["Nuisances_MarginalChi3Restraint_Sigma_" +
1177  self.label] = ([self.sigma], self.sigmamaxtrans)
1178  return ps
1179
1180  def get_output(self):
1181  self.m.update()
1182  output = {}
1183  score = self.rs.unprotected_evaluate(None)
1184  output["_TotalScore"] = str(score)
1185  output["MarginalChi3Restraint_" + self.label] = str(score)
1186  output["MarginalChi3Restraint_Sigma_" +
1187  self.label] = str(self.sigma.get_scale())
1188  return output
1189
1190 #
1191
1192
1194
1195  '''
1196  this class initializes a CrossLinkMS restraint and contains
1197  all useful information, such as the cross-link database, contained in self.pairs
1198  If restraint_file=None, it will proceed creating simulated data
1199  '''
1200
1201  def __init__(self, prots,
1202  listofxlinkertypes=["BS3", "BS2G", "EGS"], map_between_protein_names_and_chains=None,
1203  sigmamin=1.0, sigmamax=1.0, sigmagrid=1, sigmaissampled=False, typeofprofile="gofr"):
1204  global impisd2, tools
1205  import IMP.isd2 as impisd2
1206  import IMP.pmi1.tools as tools
1207
1208  if map_between_protein_names_and_chains is None:
1209  map_between_protein_names_and_chains = {}
1210
1211  self.rs = IMP.RestraintSet('data')
1212  self.rs2 = IMP.RestraintSet('prior')
1213  self.prots = prots
1214  self.label = "None"
1215  self.pairs = []
1216  self.m = self.prots[0].get_model()
1217  self.sigmamin = sigmamin
1218  self.sigmamax = sigmamax
1219  self.sigmagrid = sigmagrid
1220  self.sigmaissampled = sigmaissampled
1221
1222  self.sigmatrans = 0.1
1223  self.sigmaglobal = tools.SetupNuisance(self.m, self.sigmamin,
1224  self.sigmamin, self.sigmamax, self.sigmaissampled).get_particle()
1225  self.outputlevel = "low"
1227  # simulated reaction rates
1228  self.reaction_rates = None
1229  self.allpairs_database = None
1230  self.residue_list = None
1231
1233  # this map is used in the new cross-link file reader
1234  #{"Nsp1":"A","Nup82":"B"} etc.
1235  self.mbpnc = map_between_protein_names_and_chains
1236  # check whether the file was initialized
1237
1238  #-------------------------------
1239  def get_crosslinker_dict(self, typeofprofile="gofr"):
1240  # fill the cross-linker pmfs
1241  # to accelerate the init the list listofxlinkertypes might contain only yht needed crosslinks
1242  # type of profile can be gofr or pfes
1243
1244  disttuple = (0.0, 200.0, 1000)
1245  omegatuple = (1.0, 1000.0, 30)
1246  sigmatuple = (self.sigmamin, self.sigmamax, self.sigmagrid)
1247
1248  crosslinker_dict = {}
1249  if "BS3" in self.listofxlinkertypes:
1251  "pmf_bs3l_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1252  don=None, doff=None, prior=0, type_of_profile=typeofprofile)
1253  if "BS2G" in self.listofxlinkertypes:
1255  "pmf_bs2gl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1256  don=None, doff=None, prior=0, type_of_profile=typeofprofile)
1257  if "EGS" in self.listofxlinkertypes:
1259  "pmf_egl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple,
1260  don=None, doff=None, prior=0, type_of_profile=typeofprofile)
1261  if "Short" in self.listofxlinkertypes:
1262  # setup a "short" xl with an half length of 10 Ang
1264  10.0,
1265  disttuple,
1266  omegatuple,
1267  sigmatuple)
1269
1270  #-------------------------------
1271
1272  def add_restraints(self, restraint_file=None, oldfile=True):
1273
1274  if restraint_file is None:
1275  # get the restraints from simulated data
1276  restraint_list = self.allpairs_database
1277
1278  else:
1279  # get the restraints from external file
1280  f = open(restraint_file)
1281  restraint_list = f.readlines()
1282
1283  self.index = 0
1284
1285  self.added_pairs_list = []
1286  self.missing_residues = []
1287  for line in restraint_list:
1288  # force_restraint=True makes all intra rigid body restraint to be
1289  # accepted
1290  force_restraint = False
1291
1292  if restraint_file is None:
1293  if line["Is_Detected"]:
1295  (r1, c1) = line["Identified_Pair1"]
1296  (r2, c2) = line["Identified_Pair2"]
1297  index += 1
1298  else:
1299  continue
1300
1301  elif oldfile:
1302  tokens = line.split()
1303  # skip character
1304  if (tokens[0] == "#"):
1305  continue
1306  r1 = int(tokens[0])
1307  c1 = tokens[1]
1308  r2 = int(tokens[2])
1309  c2 = tokens[3]
1310  crosslinker = tokens[4]
1311
1312  # two restraints with the same index will be ambiguous
1313  self.index = int(tokens[5])
1314
1315  # force restraint even if it belong to the same rigid body, use
1316  # it for ambiguous restraints
1317  if (tokens[len(tokens) - 1] == "F"):
1318  force_restraint = True
1319
1320  else:
1321  # read with the new file parser
1322  totallist = eval(line)
1324  # skip the rest
1325  continue
1326
1327  print('''CrossLinkMS: attempting to add restraint between
1328  residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
1329
1330  p1s = []
1331  p2s = []
1332
1333  try:
1334  s1 = IMP.atom.Selection(
1335  self.prots[0],
1336  residue_index=r1,
1337  chains=c1,
1338  atom_type=IMP.atom.AT_CA)
1339  p1 = (s1.get_selected_particles()[0])
1340  except:
1341  print("CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1342  continue
1343  try:
1344  s2 = IMP.atom.Selection(
1345  self.prots[0],
1346  residue_index=r2,
1347  chains=c2,
1348  atom_type=IMP.atom.AT_CA)
1349  p2 = (s2.get_selected_particles()[0])
1350  except:
1351  print("CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1352  continue
1353
1354  for copy in self.prots:
1355  s1 = IMP.atom.Selection(
1356  copy,
1357  residue_index=r1,
1358  chains=c1,
1359  atom_type=IMP.atom.AT_CA)
1360  p1s.append(s1.get_selected_particles()[0])
1361  s2 = IMP.atom.Selection(
1362  copy,
1363  residue_index=r2,
1364  chains=c2,
1365  atom_type=IMP.atom.AT_CA)
1366  p2s.append(s2.get_selected_particles()[0])
1367
1368  # check whether the atom pair belongs to the same rigid body
1369  if(IMP.core.RigidMember.get_is_setup(p1s[0]) and
1371  IMP.core.RigidMember(p1s[0]).get_rigid_body() ==
1372  IMP.core.RigidMember(p2s[0]).get_rigid_body() and not force_restraint):
1373  print('''CrossLinkMS: WARNING> residue %d of chain %s and
1374  residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
1375  continue
1376
1377  # this list contains the list of symmetric pairs to avoid
1378  # duplicates
1379  if (p1s[0], p2s[0], crosslinker) in self.added_pairs_list:
1380  print("CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1381  continue
1382  if (p2s[0], p1s[0], crosslinker) in self.added_pairs_list:
1383  print("CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1384  continue
1385
1386  print("CrossLinkMS: added pair %d %s %d %s" % (r1, c1, r2, c2))
1388
1389  rs_name = 'restraint_' + str(index)
1390
1391  ln = impisd2.CrossLinkMSRestraint(
1392  self.sigmaglobal,
1394  for i in range(len(p1s)):
1396
1397  for i in range(len(self.prots)):
1398  self.pairs.append((p1s[i], p2s[i],
1400  100, 100, (r1, c1, i), (r2, c2, i),
1401  crosslinker, i, ln))
1402
1404
1407  impisd2.UniformPrior(self.sigmaglobal,
1408  1000.0,
1409  self.sigmaglobal.get_upper() - 1.0,
1410  self.sigmaglobal.get_lower() + 0.1))
1411  print("CrossLinkMS: missing residues")
1412  for ms in self.missing_residues:
1413  print("CrossLinkMS:missing " + str(ms))
1414
1415 #---------------------------------
1417  force_restraint = False
1418  ambiguous_list = totallist[0]
1419  crosslinker = totallist[1]
1420  if (totallist[2] == "F"):
1421  force_restraint = True
1422
1423  p1s = []
1424  p2s = []
1425  r1s = []
1426  r2s = []
1427  c1s = []
1428  c2s = []
1429  self.index += 1
1430  for pair in ambiguous_list:
1431  error = False
1432
1433  try:
1434  c1 = self.mbpnc[pair[0][0]]
1435  except:
1436  "CrossLinkMS: WARNING> protein name " + \
1437  pair[0][0] + " was not defined"
1438  continue
1439  try:
1440  c2 = self.mbpnc[pair[1][0]]
1441  except:
1442  "CrossLinkMS: WARNING> protein name " + \
1443  pair[1][0] + " was not defined"
1444  continue
1445  r1 = int(pair[0][1])
1446  r2 = int(pair[1][1])
1447
1448  print('''CrossLinkMS: attempting to add restraint between
1449  residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
1450
1451  try:
1452  s1 = IMP.atom.Selection(
1453  self.prots[0],
1454  residue_index=r1,
1455  chains=c1,
1456  atom_type=IMP.atom.AT_CA)
1457  p1 = (s1.get_selected_particles()[0])
1458  except:
1459  print("CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r1, c1))
1460  error = True
1461  self.missing_residues.append((r1, c1))
1462  try:
1463  s2 = IMP.atom.Selection(
1464  self.prots[0],
1465  residue_index=r2,
1466  chains=c2,
1467  atom_type=IMP.atom.AT_CA)
1468  p2 = (s2.get_selected_particles()[0])
1469  except:
1470  print("CrossLinkMS: WARNING> residue %d of chain %s is not there" % (r2, c2))
1471  error = True
1472  self.missing_residues.append((r2, c2))
1473  if error:
1474  continue
1475
1476  s1 = IMP.atom.Selection(
1477  self.prots[0],
1478  residue_index=r1,
1479  chains=c1,
1480  atom_type=IMP.atom.AT_CA)
1481  p1 = s1.get_selected_particles()[0]
1482  s2 = IMP.atom.Selection(
1483  self.prots[0],
1484  residue_index=r2,
1485  chains=c2,
1486  atom_type=IMP.atom.AT_CA)
1487  p2 = s2.get_selected_particles()[0]
1488  # this list contains the list of symmetric pairs to avoid
1489  # duplicates
1490  if (p1, p2, crosslinker) in self.added_pairs_list:
1491  print("CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1492  continue
1493  if (p2, p1, crosslinker) in self.added_pairs_list:
1494  print("CrossLinkMS: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
1495  continue
1496
1497  # check whether the atom pair belongs to the same rigid body
1500  IMP.core.RigidMember(p1).get_rigid_body() ==
1501  IMP.core.RigidMember(p2).get_rigid_body() and not force_restraint):
1502  print('''CrossLinkMS: WARNING> residue %d of chain %s and
1503  residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
1504  continue
1505
1506  p1s.append(p1)
1507  p2s.append(p2)
1508  r1s.append(r1)
1509  r2s.append(r2)
1510  c1s.append(c1)
1511  c2s.append(c2)
1512
1513  print("CrossLinkMS: added pair %d %s %d %s" % (r1, c1, r2, c2))
1514
1516
1517  if len(p1s) > 0:
1518  rs_name = '{:05}'.format(self.index % 100000)
1519
1520  ln = impisd2.CrossLinkMSRestraint(
1521  self.sigmaglobal,
1523  for i in range(len(p1s)):
1524  print(rs_name, i)
1526  self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name,
1527  100, 100, (r1s[i], c1s[i], i), (r2s[i], c2s[i], i), crosslinker, i, ln))
1529
1530 #---------------------------------
1531
1532  def simulate_data(self, crosslinker, weights, sensitivity_threshold=0.1,
1533  false_positive_half=0.02, elapsed_time=0.01,
1534  ratemin=100, ratemax=100):
1535
1536  from random import choice
1537  from random import randrange
1538  from itertools import combinations
1539  from numpy.random import binomial
1540
1541  self.weights = weights
1542  self.sensitivity_threshold = sensitivity_threshold
1543  self.false_positive_half = false_positive_half
1544  self.elapsed_time = elapsed_time
1545  self.ratemin = ratemin
1546  self.ratemax = ratemax
1547
1548  # dictionary of random reaction rates
1549  # check if they were already initialized
1550  if self.reaction_rates is None:
1551  self.reaction_rates = {}
1552
1553  s0 = IMP.atom.Selection(self.prots[0],
1554  residue_type=IMP.atom.get_residue_type(
1555  "K"),
1556  atom_type=IMP.atom.AT_CA)
1557  self.residue_list = []
1558
1559  for p1 in s0.get_selected_particles():
1560  (r1, c1) = tools.get_residue_index_and_chain_from_particle(p1)
1561  self.residue_list.append((r1, c1))
1562  if self.ratemin != self.ratemax:
1563  self.reaction_rates[(
1564  r1,
1565  c1)] = randrange(self.ratemin,
1566  self.ratemax,
1567  1)
1568  else:
1569  self.reaction_rates[(r1, c1)] = self.ratemax
1570
1571  if self.allpairs_database is None:
1572  self.allpairs_database = []
1573
1574  # generate the restraints
1575  allcomb = list(combinations(self.residue_list, 2))
1576  for ((r1, c1), (r2, c2)) in allcomb:
1577
1578  p1s = []
1579  p2s = []
1580
1581  for copy in self.prots:
1582  s1 = IMP.atom.Selection(copy, residue_index=r1,
1583  chains=c1, atom_type=IMP.atom.AT_CA)
1584  p1s.append(s1.get_selected_particles()[0])
1585  s2 = IMP.atom.Selection(copy, residue_index=r2,
1586  chains=c2, atom_type=IMP.atom.AT_CA)
1587  p2s.append(s2.get_selected_particles()[0])
1588
1589  ln = impisd2.CrossLinkMSRestraint(
1590  self.sigmaglobal,
1592  for i in range(len(p1s)):
1594  d1 = IMP.core.XYZ(p1s[i])
1595  d2 = IMP.core.XYZ(p2s[i])
1596  dist = IMP.core.get_distance(d1, d2)
1597  reactionrate1 = self.reaction_rates[(r1, c1)]
1598  reactionrate2 = self.reaction_rates[(r2, c2)]
1599  prob = ln.get_marginal_probabilities()[i]
1600  effrate = float(reactionrate1 * reactionrate2) / \
1601  (reactionrate1 + reactionrate2)
1602  probt = self.weights[i] * \
1603  (1 - exp(-effrate * prob * elapsed_time))
1604  falsepositiveprob = exp(-probt / false_positive_half)
1605  falsepositivebool = False
1606  falsepositive = binomial(n=1, p=falsepositiveprob)
1607  if (falsepositive == 1):
1608  falsepositivebool = True
1609  randompair = choice(allcomb)
1610  randpair1 = randompair[0]
1611  randpair2 = randompair[1]
1612  else:
1613  randpair1 = (r1, c1)
1614  randpair2 = (r2, c2)
1615  if (probt > sensitivity_threshold):
1616  detectedbool = True
1617  else:
1618  detectedbool = False
1619
1620  self.allpairs_database.append({})
1621  self.allpairs_database[-1]["Particle1"] = p1s[i]
1622  self.allpairs_database[-1]["Particle2"] = p2s[i]
1623  self.allpairs_database[-1]["Distance"] = dist
1625  self.allpairs_database[-1]["IMPRestraint"] = ln
1626  self.allpairs_database[-1]["IMPRestraint_Probability"] = prob
1627  self.allpairs_database[-1]["Reaction_Rate1"] = reactionrate1
1628  self.allpairs_database[-1]["Reaction_Rate2"] = reactionrate2
1629  self.allpairs_database[-1]["Effective_Rate"] = effrate
1630  self.allpairs_database[-1]["CrossLink_Fraction"] = probt
1631  self.allpairs_database[
1632  -1]["Resid1_Chainid1_Copy1"] = (r1, c1, i)
1633  self.allpairs_database[
1634  -1]["Resid2_Chainid2_Copy2"] = (r2, c2, i)
1635  self.allpairs_database[
1636  -1]["Is_False_Positive"] = falsepositivebool
1637  self.allpairs_database[-1]["Identified_Pair1"] = randpair1
1638  self.allpairs_database[-1]["Identified_Pair2"] = randpair2
1639  self.allpairs_database[-1]["Is_Detected"] = detectedbool
1640
1641  def set_hierarchy(self, prots):
1642  # we use it to change the hierarchy
1643  self.prots = prots
1644
1645  def initialize_simulated_database(self):
1646  # we use it to restart the simulation
1647  self.allpairs_database = None
1648
1649  def get_number_detected_inter(self, xl_type):
1650  # we use it to see ho many xls of a give type (eg. BS3) were detected
1651  # as inter-chain
1652  ndetected = 0
1653  for el in self.allpairs_database:
1654  if el["Is_Detected"] and \
1655  ( el["Identified_Pair1"][1] != el["Identified_Pair2"][1] ) and \
1656  el["Crosslinker"] == xl_type:
1657  ndetected += 1
1658  return ndetected
1659
1660  def get_number_detected_inter_false_positive(self, xl_type):
1661  # we use it to see ho many xls of a give type (eg. BS3) were detected
1662  # as inter-chain
1663  ndetectedfp = 0
1664  for el in self.allpairs_database:
1665  if el["Is_Detected"] and \
1666  ( el["Identified_Pair1"][1] != el["Identified_Pair2"][1] ) and \
1667  el["Crosslinker"] == xl_type and el["Is_False_Positive"]:
1668  ndetectedfp += 1
1669  return ndetectedfp
1670
1671  def show_simulated_data(self, what="Inter"):
1672  #"what" can be "All", "Detected", "FalsePositive", "TruePositive", "Intra", "Inter"
1673  if not self.allpairs_database is None:
1674  detectedlist = []
1675  for el in self.allpairs_database:
1676  printbool = False
1677  if el["Is_Detected"]:
1678  p1 = el["Identified_Pair1"]
1679  p2 = el["Identified_Pair2"]
1680  isf = el["Is_False_Positive"]
1681  isinter = (
1682  el["Identified_Pair1"][1] != el["Identified_Pair2"][1])
1683  cl = el["Crosslinker"]
1684  detectedlist.append((p1, p2, isf, cl, isinter))
1685
1686  if el["Is_Detected"] and what == "Detected":
1687  printbool = True
1688  if el["Is_Detected"] and el["Is_False_Positive"] and what == "FalsePositive":
1689  printbool = True
1690  if el["Is_Detected"] and el["Is_False_Positive"] == False and what == "TruePositive":
1691  printbool = True
1692  if el["Is_Detected"] and what == "Intra" and \
1693  (el["Identified_Pair1"][1] == el["Identified_Pair2"][1]):
1694  printbool = True
1695  if el["Is_Detected"] and what == "Inter" and \
1696  (el["Identified_Pair1"][1] != el["Identified_Pair2"][1]):
1697  printbool = True
1698  if what == "All":
1699  printbool = True
1700
1701  if printbool:
1702  print("Residue1: %6s, chainid1: %6s, copy1: %6d" % el["Resid1_Chainid1_Copy1"])
1703  print("Residue2: %6s, chainid2: %6s, copy2: %6d" % el["Resid2_Chainid2_Copy2"])
1704  keylist = list(el.keys())
1705  keylist.sort()
1706  for k in keylist:
1707  print("----", k, el[k])
1708  print("r1 c1 r2 c2 FP XL Inter")
1709  for d in detectedlist:
1710  print(d[0][0], d[0][1], d[1][0], d[1][1], d[2], d[3], d[4])
1711  else:
1712  print("CrossLinkMS: Simulated data not initialized")
1713  exit()
1714
1715  def dump_simulated_data(self, filename="simulated_cross_link.dat"):
1716  # dump the whole simulated xl database on a file
1717  sclf = open(filename, "w")
1718  for el in self.allpairs_database:
1719  sclf.write(str(el))
1720  sclf.write("\n")
1721  sclf.close()
1722
1723  def write_simulated_data(self, filename="detected_cross_link.dat"):
1724  # dump the whole simulated xl database on a file
1725  sclf = open(filename, "w")
1726  index = 0
1727  for el in self.allpairs_database:
1728  if el["Is_Detected"]:
1729  index += 1
1730  p1 = el["Identified_Pair1"]
1731  p2 = el["Identified_Pair2"]
1732  isf = el["Is_False_Positive"]
1733  cl = el["Crosslinker"]
1734  s = " ".join(
1735  [str(p1[0]),
1736  p1[1],
1737  str(p2[0]),
1738  p2[1],
1739  cl,
1740  str(index),
1741  str(isf),
1742  "T"])
1743  sclf.write(s)
1744  sclf.write("\n")
1745  sclf.close()
1746
1747  def set_label(self, label):
1748  self.label = label
1749
1753
1754  def get_hierarchies(self):
1755  return self.prots
1756
1757  def get_particles(self):
1758  return self.sigmaglobal
1759
1760  def get_restraint_sets(self):
1761  return self.rs, self.rs2
1762
1763  def get_restraint(self):
1764  tmprs = IMP.RestraintSet('xlms')
1767  return tmprs
1768
1769  def set_output_level(self, level="low"):
1770  # this might be "low" or "high"
1771  self.outputlevel = level
1772
1773  def print_chimera_pseudobonds(self, filesuffix, model=0):
1774  f = open(filesuffix + ".chimera", "w")
1775  atype = "ca"
1776  for p in self.pairs:
1777  s = "#" + str(model) + ":" + str(p[6][0]) + "." + p[6][1] + "@" + atype + \
1778  " #" + str(model) + ":" + \
1779  str(p[7][0]) + "." + p[7][1] + "@" + atype
1780  f.write(s + "\n")
1781  f.close()
1782
1783  def get_particles_to_sample(self):
1784  ps = {}
1785  if self.sigmaissampled:
1787  self.label] = ([self.sigmaglobal], self.sigmatrans)
1788  return ps
1789
1790  def get_output(self):
1791  # content of the crosslink database pairs
1792  #self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name, 100, 100, (r1,c1,i), (r2,c2,i), crosslinker, i, ln))
1793  self.m.update()
1794  output = {}
1795  score = self.rs.unprotected_evaluate(None)
1796  score2 = self.rs2.unprotected_evaluate(None)
1797  output["_TotalScore"] = str(score + score2)
1798
1799  output["CrossLinkMS_Likelihood_" + self.label] = str(score)
1800  output["CrossLinkMS_Prior_" + self.label] = str(score2)
1801  output["CrossLinkMS_Sigma"] = str(self.sigmaglobal.get_scale())
1802
1803  if self.outputlevel == "high":
1804  for i in range(len(self.pairs)):
1805
1806  p0 = self.pairs[i][0]
1807  p1 = self.pairs[i][1]
1808  crosslinker = self.pairs[i][2]
1809  ln = self.pairs[i][10]
1810  index = self.pairs[i][9]
1811  rsname = self.pairs[i][3]
1812  resid1 = self.pairs[i][6][0]
1813  chain1 = self.pairs[i][6][1]
1814  copy1 = self.pairs[i][6][2]
1815  resid2 = self.pairs[i][7][0]
1816  chain2 = self.pairs[i][7][1]
1817  copy2 = self.pairs[i][7][2]
1818  label_copy = str(rsname) + ":" + str(index) + "-" + str(resid1) + \
1819  ":" + chain1 + "_" + "-" + \
1820  str(resid2) + ":" + chain2 + "_" + crosslinker
1822  label_copy] = str(ln.get_marginal_probabilities()[index])
1823
1824  if copy1 == 0:
1825  label = str(resid1) + ":" + chain1 + \
1826  "_" + str(resid2) + ":" + chain2
1829  str(rsname) + "_" + crosslinker + "_" + label] = str(ln.unprotected_evaluate(None))
1830
1831  d0 = IMP.core.XYZ(p0)
1832  d1 = IMP.core.XYZ(p1)
1834  label_copy] = str(IMP.core.get_distance(d0, d1))
1835
1836  return output
1837
1838 #
1839
1840
1841 class BinomialXLMSRestraint(object):
1842
1843  def __init__(self, m, prots,
1844  listofxlinkertypes=["BS3", "BS2G", "EGS"], map_between_protein_names_and_chains=None, typeofprofile='pfes'):
1845
1846  if map_between_protein_names_and_chains is None:
1847  map_between_protein_names_and_chains = {}
1848
1849  global impisd2, tools, exp
1850  import IMP.isd2 as impisd2
1851  import IMP.pmi1.tools as tools
1852
1853  self.setup = 0
1854
1855  self.label = "None"
1856  self.rs = IMP.RestraintSet('xlms')
1857  self.rs2 = IMP.RestraintSet('jeffreys')
1858  self.m = m
1859  self.prots = prots
1860  self.pairs = []
1861
1862  self.weightmaxtrans = 0.05
1863  self.weightissampled = False
1864
1865  self.sigmainit = 5.0
1866  self.sigmamin = 1.0
1867  self.sigmaminnuis = 0.0
1868  self.sigmamax = 10.0
1869  self.sigmamaxnuis = 11.0
1870  self.nsigma = 100
1871  self.sigmaissampled = True
1872  self.sigmamaxtrans = 0.1
1873
1874  self.betainit = 1.0
1875  self.betamin = 1.0
1876  self.betamax = 4.0
1877  if self.setup == 1:
1878  self.betaissampled = True
1879  print("BinomialXLMSRestraint: beta is sampled")
1880  if self.setup == 0:
1881  self.betaissampled = False
1882  print("BinomialXLMSRestraint: beta is NOT sampled")
1883  self.betamaxtrans = 0.01
1884
1885  '''
1886  self.deltainit=0.001
1887  self.deltamin=0.001
1888  self.deltamax=0.1
1889  self.deltaissampled=False
1890  self.deltamaxtrans=0.001
1891
1892  self.laminit=5.0
1893  self.lammin=0.01
1894  self.lamminnuis=0.00001
1895  self.lammax=10.0
1896  self.lammaxnuis=100.0
1897  self.lamissampled=False
1898  self.lammaxtrans=0.1
1899  '''
1900
1901  self.epsilon = 0.01
1902  self.psi_dictionary = {}
1903
1904  self.sigma = tools.SetupNuisance(self.m, self.sigmainit,
1905  self.sigmaminnuis, self.sigmamaxnuis, self.sigmaissampled).get_particle()
1906
1907  self.beta = tools.SetupNuisance(self.m, self.betainit,
1908  self.betamin, self.betamax, self.betaissampled).get_particle()
1909
1910  '''
1911  self.delta=tools.SetupNuisance(self.m,self.deltainit,
1912  self.deltamin,self.deltamax,self.deltaissampled).get_particle()
1913
1914  self.lam=tools.SetupNuisance(self.m,self.laminit,
1915  self.lamminnuis,self.lammaxnuis,self.lamissampled).get_particle()
1916
1917  self.weight=tools.SetupWeight(m,False).get_particle()
1918
1919  for n in range(len(self.prots)):
1921  '''
1922  self.outputlevel = "low"
1924  # simulated reaction rates
1925  self.reaction_rates = None
1926  self.allpairs_database = None
1927  self.residue_list = None
1928
1930  # this map is used in the new cross-link file reader
1931  #{"Nsp1":"A","Nup82":"B"} etc.
1932  self.mbpnc = map_between_protein_names_and_chains
1933  # check whether the file was initialized
1934
1935  def create_psi(self, index, value):
1936  if value is None:
1937  self.psiinit = 0.01
1938  self.psiissampled = True
1939  print("BinomialXLMSRestraint: psi " + str(index) + " is sampled")
1940  else:
1941  self.psiinit = value
1942  self.psiissampled = False
1943  print("BinomialXLMSRestraint: psi " + str(index) + " is NOT sampled")
1944  self.psiminnuis = 0.0000001
1945  self.psimaxnuis = 0.4999999
1946  self.psimin = 0.01
1947  self.psimax = 0.49
1948  self.psitrans = 0.01
1949  self.psi = tools.SetupNuisance(self.m, self.psiinit,
1950  self.psiminnuis, self.psimaxnuis, self.psiissampled).get_particle()
1951  self.psi_dictionary[index] = (
1952  self.psi,
1953  self.psitrans,
1954  self.psiissampled)
1955
1956  def get_psi(self, index, value):
1957  if not index in self.psi_dictionary:
1958  self.create_psi(index, value)
1959  return self.psi_dictionary[index]
1960
1961  def get_crosslinker_dict(self, typeofprofile):
1962  # fill the cross-linker pmfs
1963  # to accelerate the init the list listofxlinkertypes might contain only
1964  # yht needed crosslinks
1965
1966  disttuple = (0.0, 200.0, 500)
1967  omegatuple = (0.01, 1000.0, 30)
1968  sigmatuple = (self.sigmamin, self.sigmamax, self.nsigma)
1969
1970  crosslinker_dict = {}
1971  if "BS3" in self.listofxlinkertypes:
1973  "pmf_bs3l_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=None, doff=None, prior=1, type_of_profile=typeofprofile)
1974  if "BS2G" in self.listofxlinkertypes:
1976  "pmf_bs2gl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=None, doff=None, prior=1, type_of_profile=typeofprofile)
1977  if "EGS" in self.listofxlinkertypes:
1979  "pmf_egl_tip3p.txt.standard", disttuple, omegatuple, sigmatuple, don=None, doff=None, prior=1, type_of_profile=typeofprofile)
1980  if "Short" in self.listofxlinkertypes:
1981  # setup a "short" xl with an half length of 10 Ang
1983  10.0,
1984  disttuple,
1985  omegatuple,
1986  sigmatuple)
1988
1989  def add_restraints(self, restraint_file=None, oldfile=False):
1990
1991  # get the restraints from external file
1992  f = open(restraint_file)
1993  restraint_list = f.readlines()
1994
1995  self.index = 0
1996
1997  self.added_pairs_list = []
1998  self.missing_residues = []
1999  for line in restraint_list:
2000  # force_restraint=True makes all intra rigid body restraint to be
2001  # accepted
2002  force_restraint = False
2003
2004  # read with the new file parser
2005  totallist = eval(line)
2007
2009  impisd2.UniformPrior(
2010  self.sigma,
2011  1000000000.0,
2012  self.sigmamax,
2013  self.sigmamin))
2015
2016  for psiindex in self.psi_dictionary:
2017  if self.psi_dictionary[psiindex][2]:
2018  psip = self.psi_dictionary[psiindex][0]
2019
2020  if self.setup == 0:
2021  print("BinomialXLMSRestraint: setup 0, adding BinomialJeffreysPrior to psi particle " + str(psiindex))
2024
2026  impisd2.UniformPrior(
2027  psip,
2028  1000000000.0,
2029  self.psimax,
2030  self.psimin))
2031
2033  # the constructor variable specify what constroctor to use in
2034  # the restraint 0: use predetermined f.p.r. (psi)
2035  # 1: infer the f.p.r. defining a psi nuisance defined by a class
2036  force_restraint = False
2037  ambiguous_list = totallist[0]
2038  crosslinker = totallist[1]
2039  if (totallist[2] == "F"):
2040  force_restraint = True
2041
2042  p1s = []
2043  p2s = []
2044  r1s = []
2045  r2s = []
2046  c1s = []
2047  c2s = []
2048  psis = []
2049  psivalues = []
2050  self.index += 1
2051  for pair in ambiguous_list:
2052  error = False
2053
2054  try:
2055  c1 = self.mbpnc[pair[0][0]]
2056  except:
2057  print("BinomialXLMSRestraint: WARNING> protein name " + pair[0][0] + " was not defined")
2058  continue
2059  try:
2060  c2 = self.mbpnc[pair[1][0]]
2061  except:
2062  print("BinomialXLMSRestraint: WARNING> protein name " + pair[1][0] + " was not defined")
2063  continue
2064
2065  r1 = int(pair[0][1])
2066  r2 = int(pair[1][1])
2067  psi = float(pair[2])
2068  try:
2069  psivalue = float(pair[3])
2070  except:
2071  psivalue = None
2072
2073  print('''CrossLinkMS: attempting to add restraint between
2074  residue %d of chain %s and residue %d of chain %s''' % (r1, c1, r2, c2))
2075
2076  try:
2077  s1 = IMP.atom.Selection(
2078  self.prots[0],
2079  residue_index=r1,
2080  chains=c1,
2081  atom_type=IMP.atom.AT_CA)
2082  p1 = (s1.get_selected_particles()[0])
2083  except:
2084  print("BinomialXLMSRestraint: WARNING> residue %d of chain %s is not there" % (r1, c1))
2085  error = True
2086  self.missing_residues.append((r1, c1))
2087  try:
2088  s2 = IMP.atom.Selection(
2089  self.prots[0],
2090  residue_index=r2,
2091  chains=c2,
2092  atom_type=IMP.atom.AT_CA)
2093  p2 = (s2.get_selected_particles()[0])
2094  except:
2095  print("BinomialXLMSRestraint: WARNING> residue %d of chain %s is not there" % (r2, c2))
2096  error = True
2097  self.missing_residues.append((r2, c2))
2098  if error:
2099  continue
2100
2101  s1 = IMP.atom.Selection(
2102  self.prots[0],
2103  residue_index=r1,
2104  chains=c1,
2105  atom_type=IMP.atom.AT_CA)
2106  p1 = s1.get_selected_particles()[0]
2107  s2 = IMP.atom.Selection(
2108  self.prots[0],
2109  residue_index=r2,
2110  chains=c2,
2111  atom_type=IMP.atom.AT_CA)
2112  p2 = s2.get_selected_particles()[0]
2113  # this list contains the list of symmetric pairs to avoid
2114  # duplicates
2115  if (p1, p2, crosslinker) in self.added_pairs_list:
2116  print("BinomialXLMSRestraint: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2117  continue
2118  if (p2, p1, crosslinker) in self.added_pairs_list:
2119  print("BinomialXLMSRestraint: WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2120  continue
2121
2122  # check whether the atom pair belongs to the same rigid body
2125  IMP.core.RigidMember(p1).get_rigid_body() ==
2126  IMP.core.RigidMember(p2).get_rigid_body() and not force_restraint):
2127  print('''BinomialXLMSRestraint: WARNING> residue %d of chain %s and
2128  residue %d of chain %s belong to the same rigid body''' % (r1, c1, r2, c2))
2129  continue
2130
2131  p1s.append(p1)
2132  p2s.append(p2)
2133  r1s.append(r1)
2134  r2s.append(r2)
2135  c1s.append(c1)
2136  c2s.append(c2)
2137  psis.append(psi)
2138  psivalues.append(psivalue)
2139
2140  print("BinomialXLMSRestraint: added pair %d %s %d %s" % (r1, c1, r2, c2))
2141
2143
2144  if len(p1s) > 0:
2145  rs_name = '{:05}'.format(self.index % 100000)
2146
2147  if self.setup == 0:
2148  print("BinomialXLMSRestraint: constructor 0")
2149  ln = impisd2.BinomialCrossLinkMSRestraint(
2150  self.m,
2151  self.sigma,
2152  self.epsilon,
2154
2155  if self.setup == 1:
2156  print("BinomialXLMSRestraint: constructor 1")
2157  ln = impisd2.BinomialCrossLinkMSRestraint(
2158  self.m,
2159  self.sigma,
2160  self.beta,
2161  self.epsilon,
2163
2164  for i in range(len(p1s)):
2166
2167  psi = self.get_psi(psis[i], psivalues[i])
2168
2170  i,
2171  (p1s[i].get_index(),
2172  p2s[i].get_index()),
2173  psi[0].get_particle().get_index())
2174  self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name,
2175  self.index, 100, (r1s[i], c1s[i], i), (r2s[i], c2s[i], i), crosslinker, i, ln))
2176
2177  h = IMP.core.Linear(0, 0.03)
2180  dps, IMP.ParticlePair(p1s[i], p2s[i]))
2182
2184
2185  print("BinomialXLMSRestraint: missing residues")
2186  for ms in self.missing_residues:
2187  print("BinomialXLMSRestraint:missing " + str(ms))
2188
2191  # exit()
2193  def set_label(self, label):
2194  self.label = label
2195
2199
2200  def get_hierarchies(self):
2201  return self.prots
2202
2203  def get_particles(self):
2204  return self.sigmaglobal
2205
2206  def get_restraint_sets(self):
2207  return self.rs, self.rs2
2208
2209  def get_restraint(self):
2210  tmprs = IMP.RestraintSet('xlms')
2213  return tmprs
2214
2215  def set_output_level(self, level="low"):
2216  # this might be "low" or "high"
2217  self.outputlevel = level
2218
2219  def print_chimera_pseudobonds(self, filesuffix, model=0):
2220  f = open(filesuffix + ".chimera", "w")
2221  atype = "ca"
2222  for p in self.pairs:
2223  s = "#" + str(model) + ":" + str(p[6][0]) + "." + p[6][1] + "@" + atype + \
2224  " #" + str(model) + ":" + \
2225  str(p[7][0]) + "." + p[7][1] + "@" + atype
2226  f.write(s + "\n")
2227  f.close()
2228
2229  def print_chimera_pseudobonds_with_psiindexes(self, filesuffix, model=0):
2230
2231  f = open(filesuffix + ".chimera", "w")
2232  atype = "ca"
2233  for p in self.pairs:
2234  s = "#" + str(model) + ":" + str(p[6][0]) + "." + p[6][1] + "@" + atype + \
2235  " #" + str(model) + ":" + \
2236  str(p[7][0]) + "." + p[7][1] + "@" + atype
2237  f.write(s + "\n")
2238  f.close()
2239
2240  def get_output(self):
2241  # content of the crosslink database pairs
2242  #self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name, psi, 100, (r1,c1,i), (r2,c2,i), crosslinker, i, ln))
2243  self.m.update()
2244  output = {}
2245  score = self.rs.unprotected_evaluate(None)
2246  score2 = self.rs2.unprotected_evaluate(None)
2247  output["_TotalScore"] = str(score + score2)
2248  output["CrossLinkMS_Likelihood_" + self.label] = str(score)
2249  output["CrossLinkMS_Prior_" + self.label] = str(score2)
2250
2251  if self.outputlevel == "high":
2252  for i in range(len(self.pairs)):
2253
2254  p0 = self.pairs[i][0]
2255  p1 = self.pairs[i][1]
2256  crosslinker = self.pairs[i][2]
2257  psiindex = self.pairs[i][4]
2258  ln = self.pairs[i][10]
2259  index = self.pairs[i][9]
2260  rsname = self.pairs[i][3]
2261  resid1 = self.pairs[i][6][0]
2262  chain1 = self.pairs[i][6][1]
2263  copy1 = self.pairs[i][6][2]
2264  resid2 = self.pairs[i][7][0]
2265  chain2 = self.pairs[i][7][1]
2266  copy2 = self.pairs[i][7][2]
2267  label_copy = str(rsname) + ":" + str(index) + "-" + str(resid1) + \
2268  ":" + chain1 + "_" + "-" + \
2269  str(resid2) + ":" + chain2 + "_" + crosslinker
2270
2271  if copy1 == 0:
2272  label = str(resid1) + ":" + chain1 + \
2273  "_" + str(resid2) + ":" + chain2
2275  output["CrossLinkMS_Score_" + str(rsname) + "_" + str(index) + "_" +
2276  crosslinker + "_" + label] = str(ln.unprotected_evaluate(None))
2277
2278  d0 = IMP.core.XYZ(p0)
2279  d1 = IMP.core.XYZ(p1)
2281  str(index) + "_" + label_copy] = str(IMP.core.get_distance(d0, d1))
2282
2283  output["CrossLinkMS_Sigma_" + self.label] = str(self.sigma.get_scale())
2284  '''
2287  '''
2288  output["CrossLinkMS_Beta_" + self.label] = str(self.beta.get_scale())
2289  for psiindex in self.psi_dictionary:
2291  str(psiindex) + "_" + self.label] = str(self.psi_dictionary[psiindex][0].get_scale())
2292  '''
2293  for n in range(self.weight.get_number_of_weights()):
2295  '''
2296  return output
2297
2298  def get_particles_to_sample(self):
2299  ps = {}
2300  if self.sigmaissampled:
2302  self.label] = ([self.sigma], self.sigmamaxtrans)
2303  '''
2304  if self.deltaissampled:
2306  if self.lamissampled:
2308  '''
2309  if self.betaissampled:
2311  self.label] = ([self.beta], self.betamaxtrans)
2312
2313  for psiindex in self.psi_dictionary:
2314  if self.psi_dictionary[psiindex][2]:
2316  str(psiindex) + "_" + self.label] = ([self.psi_dictionary[psiindex][0]], self.psi_dictionary[psiindex][1])
2317
2318  '''
2319  if self.weightissampled:
2321  '''
2322  return ps
2323
2324 #
2325
2326
2328
2329  def __init__(self, prot, restraints_file, TruncatedHarmonic=True):
2330  """read crosslink restraints between two residue
2331  of different chains from an external text file
2332  syntax: part_name_1 part_name_2 distance error
2333  example: 0 1 1.0 0.1"""
2334  self.prot = prot
2335  self.rs = IMP.RestraintSet('xlms')
2336  self.pairs = []
2337  self.m = self.prot.get_model()
2338
2340  # this is an harmonic potential with mean 12 and sigma=25, therefore
2341  # k=1/sigma^2
2342
2343  if TruncatedHarmonic:
2344  # use truncated harmonic to account for outliers
2345  hf = IMP.core.TruncatedHarmonicBound(12.0, 1.0 / 25.0, 15.0, 5)
2346  else:
2347  hf = IMP.core.Harmonic(12.0, 1.0 / 25.0)
2348  dps = IMP.core.DistancePairScore(hf)
2349
2350  # small linear contribution for long range
2351  h = IMP.core.Linear(0, 0.03)
2352  dps2 = IMP.core.DistancePairScore(h)
2353
2354  index = 0
2355
2356  addedd_pairs_list = []
2357  for line in open(restraints_file):
2358
2359  # force_restraint=True makes all intra rigid body restraint to be
2360  # accepted
2361  force_restraint = True
2362  tokens = line.split()
2363
2364  # skip character
2365  if (tokens[0] == "#"):
2366  continue
2367  r1 = int(tokens[0])
2368  c1 = tokens[1]
2369  r2 = int(tokens[2])
2370  c2 = tokens[3]
2371  crosslinker = tokens[4]
2372
2373  # two restraints with the same index will be ambiguous
2374  index = int(tokens[5])
2375
2376  # force restraint even if it belong to the same rigid body, use it
2377  # for ambiguous restraints
2378  if (tokens[len(tokens) - 1] == "F"):
2379  force_restraint = True
2380
2381  print("attempting to add restraint between residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
2382
2383  p1s = []
2384  p2s = []
2385
2386  # apply the cross-link to the main copy
2387  try:
2388  s1 = IMP.atom.Selection(
2389  self.prot,
2390  residue_index=r1,
2391  chains=c1,
2392  atom_type=IMP.atom.AT_CA)
2393  p1 = (s1.get_selected_particles()[0])
2394  except:
2395  print("WARNING> residue %d of chain %s is not there" % (r1, c1))
2396  continue
2397  try:
2398  s2 = IMP.atom.Selection(
2399  self.prot,
2400  residue_index=r2,
2401  chains=c2,
2402  atom_type=IMP.atom.AT_CA)
2403  p2 = (s2.get_selected_particles()[0])
2404  except:
2405  print("WARNING> residue %d of chain %s is not there" % (r2, c2))
2406  continue
2407
2408  # part1=[]
2409  # part2=[]
2410  # this list contain the list of symmetric pairs to avoid
2411  # duplications
2412
2413  # this is the map between particle pairs and the restraints (there
2414  # might be more than one restraint per particle pair if you have
2415  # more than one cross-link type
2416
2417  print("attempting to add restraint between residue %d of chain %s and residue %d of chain %s" % (r1, c1, r2, c2))
2418
2419  # check whether the atom pair belongs to the same rigid body
2421  IMP.core.RigidMember(p1).get_rigid_body() == IMP.core.RigidMember(p2).get_rigid_body() and not force_restraint):
2422  print("WARNING> residue %d of chain %s and residue %d of chain %s belong to the same rigid body" % (r1, c1, r2, c2))
2423  continue
2424
2425  # this list contain the list of symmetric pairs to avoid
2426  # duplications
2427  if (p1, p2, crosslinker) in addedd_pairs_list:
2428  print("WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2429  continue
2430  if (p2, p1, crosslinker) in addedd_pairs_list:
2431  print("WARNING> pair %d %s %d %s already there" % (r1, c1, r2, c2))
2432  continue
2433
2434  print("added pair %d %s %d %s" % (r1, c1, r2, c2))
2435  index += 1
2437
2438  rs_name = 'restraint_' + str(index)
2439
2440  ln = IMP.core.PairRestraint(dps, IMP.ParticlePair(p1, p2))
2441  ln.set_name("CrossLinkMSSimple_" + str(r1)
2442  + ":" + str(c1) + "-" + str(r2) + ":" + str(c2))
2443  ln.set_weight(1.0)
2444
2446
2447  pr = IMP.core.PairRestraint(dps2, IMP.ParticlePair(p1, p2))
2448
2450
2451  self.pairs.append(
2452  (p1,
2453  p2,
2455  rs_name,
2456  100,
2457  100,
2458  (r1,
2459  c1),
2460  (r2,
2461  c2),
2463  ln,
2464  pr))
2465
2466  def set_label(self, label):
2467  self.label = label
2468
2471
2472  def get_restraint(self):
2473  return self.rs
2474
2475  def print_chimera_pseudobonds(self, filesuffix, model=0):
2476  f = open(filesuffix + ".chimera", "w")
2477  atype = "ca"
2478  for p in self.pairs:
2479  s = "#" + str(model) + ":" + str(p[6][0]) + "." + p[6][1] + "@" + atype + \
2480  " #" + str(model) + ":" + \
2481  str(p[7][0]) + "." + p[7][1] + "@" + atype
2482  f.write(s + "\n")
2483  f.close()
2484
2485  def get_output(self):
2486  # content of the crosslink database pairs
2487  #self.pairs.append((p1s[i], p2s[i], crosslinker, rs_name, 100, 100, (r1,c1), (r2,c2), crosslinker, ln,pr))
2488  self.m.update()
2489  output = {}
2490  score = self.rs.unprotected_evaluate(None)
2491  output["_TotalScore"] = str(score)
2492  output["CrossLinkMSSimple_Score_" + self.label] = str(score)
2493  for i in range(len(self.pairs)):
2494
2495  p0 = self.pairs[i][0]
2496  p1 = self.pairs[i][1]
2497  crosslinker = self.pairs[i][2]
2498  ln = self.pairs[i][9]
2499  pr = self.pairs[i][10]
2500  resid1 = self.pairs[i][6][0]
2501  chain1 = self.pairs[i][6][1]
2502  resid2 = self.pairs[i][7][0]
2503  chain2 = self.pairs[i][7][1]
2504
2505  label = str(resid1) + ":" + chain1 + "_" + \
2506  str(resid2) + ":" + chain2
2508  "_" + label] = str(ln.unprotected_evaluate(None))
2510  "_" + label] = str(pr.unprotected_evaluate(None))
2511  d0 = IMP.core.XYZ(p0)
2512  d1 = IMP.core.XYZ(p1)
2514  label] = str(IMP.core.get_distance(d0, d1))
2515
2516  return output
2517
2518
2519 def get_residue_index_and_chain_from_particle(p):
2520  rind = IMP.atom.Residue(IMP.atom.Atom(p).get_parent()).get_index()
2521  c = IMP.atom.Residue(IMP.atom.Atom(p).get_parent()).get_parent()
2522  cid = IMP.atom.Chain(c).get_id()
2523  return rind, cid
2524
2525
2526 def select_calpha_or_residue(
2527  prot,
2528  chain,
2529  resid,
2530  ObjectName="None:",
2531  SelectResidue=False):
2532  # use calphas
2533  p = None
2534  s = IMP.atom.Selection(prot, chains=chain,
2535  residue_index=resid, atom_type=IMP.atom.AT_CA)
2536
2537  ps = s.get_selected_particles()
2538  # check if the calpha selection is empty
2539  if ps:
2540  if len(ps) == 1:
2541  p = ps[0]
2542  else:
2543  print(ObjectName + " multiple residues selected for selection residue %s chain %s " % (resid, chain))
2544  else:
2545  # use the residue, in case of simplified representation
2546  s = IMP.atom.Selection(prot, chains=chain,
2547  residue_index=resid)
2548  ps = s.get_selected_particles()
2549  # check if the residue selection is empty
2550  if ps:
2551  if len(ps) == 1:
2552  p = ps[0]
2553  else:
2554  print(ObjectName + " multiple residues selected for selection residue %s chain %s " % (resid, chain))
2555
2556  else:
2557  print(ObjectName + " residue %s chain %s does not exist" % (resid, chain))
2558  return p
2559
2560
2561
2562
2563 class SetupMembranePoreRestraint(object):
2564  import math
2565
2566  def __init__(
2567  self,
2568  representation,
2569  selection_tuples_outside=None,
2570  selection_tuples_membrane=None,
2571  selection_tuples_inside=None,
2572  center=(0.0,0.0,0.0),
2573  z_tickness=100,
2575  membrane_tickness=40.0,
2576  resolution=1,
2577  label="None"):
2578
2579  self.weight = 1.0
2580  self.label = label
2581  self.m = representation.prot.get_model()
2582  self.rs = IMP.RestraintSet(self.m, label)
2583  self.representation=representation
2584  self.center=center
2586  self.z_tickness=z_tickness
2587  self.membrane_tickness=membrane_tickness
2588
2589
2590
2591
2592
2593  def create_representation(self):
2594  p=IMP.Particle(self.m)
2596  h.set_name("Membrane_Pore"+self.label)
2598
2599  inner=self.z_tickness/2-self.membrane_tickness/2
2600  outer=self.z_tickness/2+self.membrane_tickness/2
2601
2602  surface_1=[]
2603  surface_2=[]
2604  particles_surface_1=[]
2605  particles_surface_2=[]
2606
2607  for i in range(10):
2608  for j in range(10):
2609  v=self.math.pi/10.0*float(i)+self.math.pi/2.0
2610  u=2.0*self.math.pi/10.0*float(j)
2613  z=inner*self.math.sin(v)
2614  p2=IMP.Particle(self.m)
2617  d.set_coordinates((x,y,z))
2621  surface_1.append((x,y,z))
2622  particles_surface_1.append(p2)
2623
2626  z=outer*self.math.sin(v)
2627  p2=IMP.Particle(self.m)
2630  d.set_coordinates((x,y,z))
2634  surface_2.append((x,y,z))
2635  particles_surface_2.append(p2)
2636
2637  for i in range(10):
2638  for j in range(10):
2642  continue
2643  else:
2644  for n,z in enumerate([self.z_tickness/2+self.membrane_tickness/2,
2645  self.z_tickness/2-self.membrane_tickness/2,
2646  -self.z_tickness/2+self.membrane_tickness/2,
2647  -self.z_tickness/2-self.membrane_tickness/2]):
2648  p2=IMP.Particle(self.m)
2651  d.set_coordinates((x,y,z))
2655  if n == 0 or n == 3:
2656  surface_1.append((x,y,z))
2657  particles_surface_1.append(p2)
2658  if n == 1 or n == 2:
2659  surface_2.append((x,y,z))
2660  particles_surface_2.append(p2)
2661
2662  from scipy.spatial import Delaunay
2663  tri = Delaunay(surface_1)
2664
2665  '''
2666  edge_points = []
2667  edges = set()
2668
2669  def add_edge(i, j):
2670  """Add a line between the i-th and j-th points, if not in the list already"""
2671  if (i, j) in edges or (j, i) in edges:
2673  return
2674  edges.add( (i, j) )
2675  edge_points.append(points[ [i, j] ])
2676
2677  # loop over triangles:
2678  # ia, ib, ic = indices of corner points of the triangle
2679  for ia, ib, ic in tri.vertices:
2683  '''
2684
2685  '''
2686  edges=set()
2687
2688  for ia, ib, ic, id in tri.simplices:
2692
2693  for e in edges:
2694  p1=particles_surface_1[e[0]]
2695  p2=particles_surface_1[e[1]]
2696  print(p1,p2,e[0],e[1])
2697  IMP.atom.Bonded.setup_particle(p1)
2698  IMP.atom.Bonded.setup_particle(p2)
2699  IMP.atom.create_bond(IMP.atom.Bonded(p1),IMP.atom.Bonded(p2),1)
2700  '''
2701
2702  for i in range(len(particles_surface_1)-1):
2703  p1=particles_surface_1[i]
2704  p2=particles_surface_1[i+1]
2705  if not IMP.atom.Bonded.get_is_setup(p1):
2707  if not IMP.atom.Bonded.get_is_setup(p2):
A filter which returns true if a container contains the Pair.
double get_volume_from_residue_type(ResidueType rt)
Return an estimate for the volume of a given residue.
Simple 3D transformation class.
A function that is harmonic over an bounded interval.
void show_molecular_hierarchy(Hierarchy h)
Print out the molecular hierarchy.
static Atom setup_particle(Model *m, ParticleIndex pi, Atom other)
Definition: atom/Atom.h:246
static Fragment setup_particle(Model *m, ParticleIndex pi)
Definition: Fragment.h:67
Upper bound harmonic function (non-zero when feature > mean)
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
A decorator for a particle which has bonds.
Select atoms which are selected by both selectors.
Definition: pdb.h:417
A class to store a fixed array of same-typed values.
Definition: Array.h:40
def deprecated_module
Mark a Python module as deprecated.
Definition: __init__.py:9515
double get_drmsd(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1)
Calculate distance the root mean square deviation between two sets of 3D points.
Definition: atom/distance.h:49
Rotation3D get_random_rotation_3d(const Rotation3D &center, double distance)
Pick a rotation at random near the provided one.
Color get_rgb_color(double f)
Return the color for f from the RGB color map.
Return the radius of a sphere with a given volume.
Definition: Sphere3D.h:35
ParticlesTemp get_particles(Model *m, const ParticleIndexes &ps)
Get the particles from a list of indexes.
Classes to handle different kinds of restraints.
static Residue setup_particle(Model *m, ParticleIndex pi, ResidueType t, int index, int insertion_code)
Definition: Residue.h:160
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:749
Return all spatially-proximal pairs of particles (a,b) from the two SingletonContainers A and B...
static XYZ setup_particle(Model *m, ParticleIndex pi)
Definition: XYZ.h:51
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
void read_pdb(TextInput input, int model, Hierarchy h)
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
A reference frame in 3D.
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
Definition: RestraintSet.h:41
static Molecule setup_particle(Model *m, ParticleIndex pi)
Definition: Molecule.h:37
Store a list of ParticleIndexPairs.
Miscellaneous utilities.
Definition: /tools.py:1
Select all non-alternative ATOM records.
Definition: pdb.h:128
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Store a list of ParticleIndexes.
A decorator for a particle representing an atom.
Definition: atom/Atom.h:238
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Definition: Mass.h:48
static Bonded setup_particle(Model *m, ParticleIndex pi)
void add_restraints(RMF::FileHandle fh, const Restraints &hs)
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
Definition: Colored.h:62
Linear function
Definition: Linear.h:22
this class initializes a CrossLinkMS restraint and contains all useful information, such as the cross-link database, contained in self.pairs If restraint_file=None, it will proceed creating simulated data
Score a pair of particles based on the distance between their centers.
A decorator for a residue.
Definition: Residue.h:137
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
VectorD< 3 > Vector3D
Definition: VectorD.h:408
Prevent a set of particles and rigid bodies from inter-penetrating.
Class to handle individual particles of a Model object.
Definition: Particle.h:43
Store info for a chain of a protein.
Definition: Chain.h:61
Applies a PairScore to a Pair.
Definition: PairRestraint.h:31
Select all CA ATOM records.
Definition: pdb.h:142
Functionality for loading, creating, manipulating and scoring atomic structures.
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Definition: Chain.h:83
Hierarchies get_leaves(const Selection &h)
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
static RigidBody setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor ps)
Definition: rigid_bodies.h:180
Applies a PairScore to each Pair in a list.
ResidueType get_residue_type(char c)
Get the residue type from the 1-letter amino acid code.
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:27
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27