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