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