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