IMP logo
IMP Reference Guide  2.8.0
The Integrative Modeling Platform
pmi/restraints/proteomics.py
1 """@namespace IMP.pmi.restraints.proteomics
2 Restraints for handling various kinds of proteomics data.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 import IMP.algebra
9 import IMP.atom
10 import IMP.container
11 import IMP.pmi
12 import IMP.pmi.tools
13 import IMP.pmi.output
14 import numpy
15 import math
16 import math
17 import sys
18 
19 
20 class ConnectivityRestraint(object):
21 
22  '''
23  generate a connectivity restraint between domains
24  setting up the restraint
25  example:
26  cr=restraints.ConnectivityRestraint(simo,["CCC",(1,100,"TTT"),(100,150,"AAA")])
27  cr.add_to_model()
28  cr.set_label("CR1")
29  Multistate support =No
30  Selection type=selection tuple
31  Resolution=Yes
32  '''
33 
34  def __init__(
35  self,
36  representation,
37  selection_tuples,
38  kappa=10.0,
39  resolution=None,
40  label="None"):
41 
42  self.weight = 1.0
43  self.kappa = kappa
44  self.label = label
45  if self.label == "None":
46  self.label = str(selection_tuples)
47  self.m = representation.prot.get_model()
48  self.rs = IMP.RestraintSet(self.m, label)
49 
50  sels = []
51 
52  for s in selection_tuples:
53  particles = IMP.pmi.tools.select_by_tuple(representation, s,
54  resolution=resolution, name_is_ambiguous=True)
55  sel = IMP.atom.Selection(particles)
56  sels.append(sel)
57 
59  sels,
60  self.kappa,
61  self.label)
62  self.rs.add_restraint(cr)
63 
64  def set_label(self, label):
65  self.label = label
66  self.rs.set_name(label)
67  for r in self.rs.get_restraints():
68  r.set_name(label)
69 
70  def add_to_model(self):
72 
73  def get_restraint(self):
74  return self.rs
75 
76  def get_restraints(self):
77  rlist = []
78  for r in self.rs.get_restraints():
79  rlist.append(IMP.core.PairRestraint.get_from(r))
80  return rlist
81 
82  def set_weight(self, weight):
83  self.weight = weight
84  self.rs.set_weight(weight)
85 
86  def get_output(self):
87  self.m.update()
88  output = {}
89  score = self.weight * self.rs.unprotected_evaluate(None)
90  output["_TotalScore"] = str(score)
91  output["ConnectivityRestraint_" + self.label] = str(score)
92  return output
93 
94 #
95 
96 
97 class CompositeRestraint(object):
98 
99  '''
100  handleparticles is a selection tuple
101  compositeparticles is a list of selection tuples
102  '''
103 
104  def __init__(
105  self,
106  representation,
107  handleparticles_tuples,
108  compositeparticles_tuple_list,
109  cut_off=5.0,
110  lam=1.0,
111  plateau=0.0,
112  resolution=None,
113  label="None"):
114 
115  # composite particles: all particles beside the handle
116  self.label = label
117  self.m = representation.prot.get_model()
118  self.rs = IMP.RestraintSet(self.m, 'cr')
119 
120  self.handleparticles = []
121  for s in handleparticles_tuples:
122  self.handleparticles += IMP.pmi.tools.select_by_tuple(
123  representation, s,
124  resolution=resolution, name_is_ambiguous=True)
125  self.compositeparticles = []
126  compositeparticle_list = []
127  for list in compositeparticles_tuple_list:
128  tmplist = []
129  for s in list:
130  tmplist += IMP.pmi.tools.select_by_tuple(
131  representation, s,
132  resolution=resolution, name_is_ambiguous=True)
133  compositeparticle_list.append(tmplist)
134  self.compositeparticles += tmplist
135 
137  self.m,
138  self.handleparticles,
139  cut_off,
140  lam,
141  True,
142  plateau)
143 
144  for ps in compositeparticle_list:
145  # composite particles is a list of list of particles
146  ln.add_composite_particle(ps)
147 
148  self.rs.add_restraint(ln)
149 
150  def set_label(self, label):
151  self.label = label
152 
153  def get_handle_particles(self):
154  return self.handleparticles
155 
156  def get_composite_particles(self):
157  return self.compositeparticles
158 
159  def get_restraint(self):
160  return self.rs
161 
162  def add_to_model(self):
163  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
164 
165  def get_output(self):
166  self.m.update()
167  output = {}
168  score = self.rs.unprotected_evaluate(None)
169  output["_TotalScore"] = str(score)
170  output["CompositeRestraint_" + self.label] = str(score)
171  return output
172 
173 
174 #
176 
177  '''
178  this restraint allows ambiguous crosslinking between multiple copies
179  excluding between symmetric copies
180  It allows name ambiguity
181  '''
182 
183  def __init__(
184  self,
185  representation,
186  restraints_file,
187  cut_off=5.0,
188  lam=1.0,
189  plateau=0.01,
190  resolution=None,
191  label="None"):
192 
193  self.weight = 1.0
194  self.m = representation.prot.get_model()
195  self.rs = IMP.RestraintSet(self.m, 'data')
196  self.label = "None"
197  self.pairs = []
198 
199  self.outputlevel = "low"
200  self.cut_off = cut_off
201  self.lam = lam
202  self.plateau = plateau
203 
204  fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
205 
206  for line in fl:
207 
208  tokens = line.split()
209  # skip character
210  if (tokens[0] == "#"):
211  continue
212  r1 = int(tokens[2])
213  c1 = tokens[0]
214  r2 = int(tokens[3])
215  c2 = tokens[1]
216 
217  ps1 = IMP.pmi.tools.select(
218  representation,
219  resolution=resolution,
220  name=c1,
221  name_is_ambiguous=True,
222  residue=r1)
223  hrc1 = [representation.hier_db.particle_to_name[p] for p in ps1]
224  ps1nosym = [
225  p for p in ps1 if IMP.pmi.Symmetric(
226  p).get_symmetric(
227  ) == 0]
228  hrc1nosym = [representation.hier_db.particle_to_name[p]
229  for p in ps1nosym]
230 
231  if len(ps1) == 0:
232  print(
233  "AmbiguousCompositeRestraint: WARNING> residue %d of chain %s is not there" %
234  (r1, c1))
235  continue
236 
237  ps2 = IMP.pmi.tools.select(
238  representation,
239  resolution=resolution,
240  name=c2,
241  name_is_ambiguous=True,
242  residue=r2)
243  hrc2 = [representation.hier_db.particle_to_name[p] for p in ps2]
244  ps2nosym = [
245  p for p in ps2 if IMP.pmi.Symmetric(
246  p).get_symmetric(
247  ) == 0]
248  hrc2nosym = [representation.hier_db.particle_to_name[p]
249  for p in ps2nosym]
250 
251  if len(ps2) == 0:
252  print(
253  "AmbiguousCompositeRestraint: WARNING> residue %d of chain %s is not there" %
254  (r2, c2))
255  continue
256 
258  self.m,
259  ps1nosym,
260  self.cut_off,
261  self.lam,
262  True,
263  self.plateau)
264  cr.add_composite_particle(ps2)
265 
266  self.rs.add_restraint(cr)
267  self.pairs.append(
268  (ps1nosym,
269  hrc1nosym,
270  c1,
271  r1,
272  ps2,
273  hrc2,
274  c2,
275  r2,
276  cr))
277 
279  self.m,
280  ps1,
281  self.cut_off,
282  self.lam,
283  True,
284  self.plateau)
285  cr.add_composite_particle(ps2nosym)
286 
287  self.rs.add_restraint(cr)
288  self.pairs.append(
289  (ps1,
290  hrc1,
291  c1,
292  r1,
293  ps2nosym,
294  hrc2nosym,
295  c2,
296  r2,
297  cr))
298 
299  def plot_restraint(
300  self,
301  maxdist=100,
302  npoints=100):
303 
304  p1 = IMP.Particle(self.m)
305  p2 = IMP.Particle(self.m)
309  self.m,
310  [p1],
311  self.cut_off,
312  self.lam,
313  True,
314  self.plateau)
315  cr.add_composite_particle([p2])
316  dists = []
317  scores = []
318  for i in range(npoints):
319  d2.set_coordinates(
320  IMP.algebra.Vector3D(maxdist / npoints * float(i), 0, 0))
321  dists.append(IMP.core.get_distance(d1, d2))
322  scores.append(cr.unprotected_evaluate(None))
323  IMP.pmi.output.plot_xy_data(dists, scores)
324 
325  def set_label(self, label):
326  self.label = label
327  self.rs.set_name(label)
328  for r in self.rs.get_restraints():
329  r.set_name(label)
330 
331  def add_to_model(self):
332  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
333 
334  def get_hierarchies(self):
335  return self.prot
336 
337  def get_restraint_sets(self):
338  return self.rs
339 
340  def get_restraint(self):
341  return self.rs
342 
343  def set_output_level(self, level="low"):
344  # this might be "low" or "high"
345  self.outputlevel = level
346 
347  def set_weight(self, weight):
348  self.weight = weight
349  self.rs.set_weight(weight)
350 
351  def get_output(self):
352  # content of the crosslink database pairs
353  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
354  self.m.update()
355 
356  output = {}
357  score = self.weight * self.rs.unprotected_evaluate(None)
358  output["_TotalScore"] = str(score)
359  output["AmbiguousCompositeRestraint_Score_" + self.label] = str(score)
360  for n, p in enumerate(self.pairs):
361 
362  ps1 = p[0]
363  hrc1 = p[1]
364  c1 = p[2]
365  r1 = p[3]
366  ps2 = p[4]
367  hrc2 = p[5]
368  c2 = p[6]
369  r2 = p[7]
370  cr = p[8]
371  for n1, p1 in enumerate(ps1):
372  name1 = hrc1[n1]
373 
374  for n2, p2 in enumerate(ps2):
375  name2 = hrc2[n2]
376  d1 = IMP.core.XYZR(p1)
377  d2 = IMP.core.XYZR(p2)
378  label = str(r1) + ":" + name1 + "_" + str(r2) + ":" + name2
379  output["AmbiguousCompositeRestraint_Distance_" +
380  label] = str(IMP.core.get_distance(d1, d2))
381 
382  label = str(r1) + ":" + c1 + "_" + str(r2) + ":" + c2
383  output["AmbiguousCompositeRestraint_Score_" +
384  label] = str(self.weight * cr.unprotected_evaluate(None))
385 
386  return output
387 
388 
389 #
390 class SimplifiedPEMAP(object):
391 
392  def __init__(
393  self,
394  representation,
395  restraints_file,
396  expdistance,
397  strength,
398  resolution=None):
399 
400  self.m = representation.prot.get_model()
401  self.rs = IMP.RestraintSet(self.m, 'data')
402  self.label = "None"
403  self.pairs = []
404 
405  self.outputlevel = "low"
406  self.expdistance = expdistance
407  self.strength = strength
408 
409  fl = IMP.pmi.tools.open_file_or_inline_text(restraints_file)
410 
411  for line in fl:
412 
413  tokens = line.split()
414  # skip character
415  if (tokens[0] == "#"):
416  continue
417  r1 = int(tokens[2])
418  c1 = tokens[0]
419  r2 = int(tokens[3])
420  c2 = tokens[1]
421  pcc = float(tokens[4])
422 
423  ps1 = IMP.pmi.tools.select(
424  representation,
425  resolution=resolution,
426  name=c1,
427  name_is_ambiguous=False,
428  residue=r1)
429  if len(ps1) == 0:
430  print(
431  "SimplifiedPEMAP: WARNING> residue %d of chain %s is not there (w/ %d %s)" %
432  (r1, c1, r2, c2))
433  continue
434  if len(ps1) > 1:
435  print(
436  "SimplifiedPEMAP: WARNING> residue %d of chain %s selected multiple particles" %
437  (r1, c1))
438  continue
439 
440  ps2 = IMP.pmi.tools.select(
441  representation,
442  resolution=resolution,
443  name=c2,
444  name_is_ambiguous=False,
445  residue=r2)
446  if len(ps2) == 0:
447  print(
448  "SimplifiedPEMAP: WARNING> residue %d of chain %s is not there (w/ %d %s)" %
449  (r1, c1, r2, c2))
450  continue
451  if len(ps2) > 1:
452  print(
453  "SimplifiedPEMAP: WARNING> residue %d of chain %s selected multiple particles" %
454  (r2, c2))
455  continue
456 
457  p1 = ps1[0]
458  p2 = ps2[0]
459 
460  # This is harmonic potential for the pE-MAP data
461  upperdist = self.get_upper_bond(pcc)
462  limit = self.strength * (upperdist + 15) ** 2 + 10.0
464  upperdist,
465  self.strength,
466  upperdist +
467  15,
468  limit)
469 
470  # This is harmonic for the X-link
471  #hub= IMP.core.TruncatedHarmonicBound(17.0,self.strength,upperdist+15,limit)
472 
474  dr = IMP.core.PairRestraint(df, (p1, p2))
475  self.rs.add_restraint(dr)
476  self.pairs.append((p1, p2, dr, r1, c1, r2, c2))
477 
478  # Lower-bound restraint
479  lowerdist = self.get_lower_bond(pcc)
480  limit = self.strength * (lowerdist - 15) ** 2 + 10.0
482  lowerdist,
483  self.strength,
484  lowerdist +
485  15,
486  limit)
487 
488  # This is harmonic for the X-link
489  #hub2= IMP.core.TruncatedHarmonicBound(17.0,self.strength,upperdist+15,limit)
490 
492  dr2 = IMP.core.PairRestraint(df2, (p1, p2))
493  self.rs.add_restraint(dr2)
494  self.pairs.append((p1, p2, dr2, r1, c1, r2, c2))
495 
496  def get_upper_bond(self, pearsoncc):
497  # return (pearsoncc-1.)/-0.0075
498  return (pearsoncc - .5) / (-0.005415)
499 
500  def get_lower_bond(self, pearsoncc):
501  return (pearsoncc - 1.) / -0.0551
502 
503  def set_label(self, label):
504  self.label = label
505 
506  def add_to_model(self):
507  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
508 
509  def get_hierarchies(self):
510  return self.prot
511 
512  def get_restraint_sets(self):
513  return self.rs
514 
515  def set_output_level(self, level="low"):
516  # this might be "low" or "high"
517  self.outputlevel = level
518 
519  def get_output(self):
520  # content of the crosslink database pairs
521  # self.pairs.append((p1,p2,dr,r1,c1,r2,c2))
522  self.m.update()
523 
524  output = {}
525  score = self.rs.unprotected_evaluate(None)
526  output["_TotalScore"] = str(score)
527  output["SimplifiedPEMAP_Score_" + self.label] = str(score)
528  for i in range(len(self.pairs)):
529 
530  p0 = self.pairs[i][0]
531  p1 = self.pairs[i][1]
532  crosslinker = 'standard'
533  ln = self.pairs[i][2]
534  resid1 = self.pairs[i][3]
535  chain1 = self.pairs[i][4]
536  resid2 = self.pairs[i][5]
537  chain2 = self.pairs[i][6]
538 
539  label = str(resid1) + ":" + chain1 + "_" + \
540  str(resid2) + ":" + chain2
541  output["SimplifiedPEMAP_Score_" + crosslinker + "_" +
542  label] = str(ln.unprotected_evaluate(None))
543 
544  d0 = IMP.core.XYZ(p0)
545  d1 = IMP.core.XYZ(p1)
546  output["SimplifiedPEMAP_Distance_" +
547  label] = str(IMP.core.get_distance(d0, d1))
548 
549  return output
550 
551 
553 
554  '''
555  generates and wraps a IMP.pmi.ConnectivityRestraint between domains
556  example:
557  cr=restraints.ConnectivityNetworkRestraint(simo,["CCC",(1,100,"TTT"),(100,150,"AAA")])
558  cr.add_to_model()
559  cr.set_label("CR1")
560 
561  Multistate support =No
562  Selection type=selection tuple
563  Resolution=Yes
564  '''
565 
566  def __init__(
567  self,
568  representation=None,
569  objects=None,
570  kappa=10.0,
571  resolution=1.0,
572  label="None"):
573 
574  self.weight = 1.0
575  self.kappa = kappa
576  self.label = label
577  if self.label == "None":
578  self.label = str(selection_tuples)
579 
580  hiers=[]
581 
582  if representation is None:
583  print(objects)
584  for obj in objects:
585  hiers.append(IMP.pmi.tools.input_adaptor(obj,
586  resolution,
587  flatten=True))
588  self.m=hiers[0][0].get_model()
589  else:
590  self.m = representation.m
591  for s in objects:
592  hiers.append(IMP.pmi.tools.select_by_tuple(representation, s,
593  resolution=resolution,
594  name_is_ambiguous=False))
595 
596  #particles=[h.get_particle() for h in hiers]
597  cr = ConnectivityNetworkRestraint(self.m)
598  for hs in hiers:
599  cr.add_particles([h.get_particle() for h in hs])
600  self.rs = IMP.RestraintSet(self.m, label)
601  self.rs.add_restraint(cr)
602 
603  def set_label(self, label):
604  self.label = label
605  self.rs.set_name(label)
606  for r in self.rs.get_restraints():
607  r.set_name(label)
608 
609  def add_to_model(self):
610  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
611 
612  def get_restraint(self):
613  return self.rs
614 
615  def get_restraints(self):
616  rlist = []
617  for r in self.rs.get_restraints():
618  rlist.append(IMP.core.PairRestraint.get_from(r))
619  return rlist
620 
621  def set_weight(self, weight):
622  self.weight = weight
623  self.rs.set_weight(weight)
624 
625  def get_output(self):
626  self.m.update()
627  output = {}
628  score = self.weight * self.rs.unprotected_evaluate(None)
629  output["_TotalScore"] = str(score)
630  output["ConnectivityNetworkRestraint_" + self.label] = str(score)
631  return output
632 
633 
635 
636  '''
637  a python restraint that computes the score for a composite of proteins
638  Authors: G. Bouvier, R. Pellarin. Pasteur Institute.
639  '''
640 
641  def __init__(self, m, slope=1.0, theta=0.0, plateau=0.0000000001, linear_slope=0.015):
642  '''
643  input a list of particles, the slope and theta of the sigmoid potential
644  theta is the cutoff distance for a protein-protein contact
645  '''
646  # Import networkx here so that we don't introduce it as a dependency
647  # for *every* proteomics restraint, only this one
648  import networkx
649  self.networkx = networkx
650  IMP.Restraint.__init__(self, m, "ConnectivityNetworkRestraint %1%")
651  self.slope = slope
652  self.theta = theta
653  self.linear_slope = linear_slope
654  self.plateau = plateau
655  self.particles_blocks = []
656  self.particle_list = []
657 
658  def get_number_of_particle_blocks(self):
659  return len(self.particles_blocks)
660 
661  def get_number_of_particles_for_block(self, block_index):
662  return len(self.particles_blocks[block_index])
663 
664  def add_particles(self, particles):
665  self.particles_blocks.append(particles)
666  self.particle_list += particles
667 
668  def get_full_graph(self):
669  '''
670  get the full graph of distances between every particle pair
671  '''
672  import scipy.spatial
673  pdist_array = numpy.array(
675  pdist_mat = scipy.spatial.distance.squareform(pdist_array)
676  pdist_mat[pdist_mat < 0] = 0
677  graph = self.networkx.Graph(pdist_mat)
678  return graph
679 
681  """
682  return the minimum spanning tree
683  """
684  graph = self.get_full_graph()
685  graph = self.networkx.minimum_spanning_tree(graph)
686  return graph
687 
688  def sigmoid(self, x):
689  '''
690  a sigmoid function that scores the probability of a contact
691  between two proteins
692  '''
693  # return 1 - (x)**self.slope/ float(((x)**self.slope +
694  # self.theta**self.slope))
695  argvalue = (x - self.theta) / self.slope
696  return 1.0 - (1.0 - self.plateau) / (1.0 + math.exp(-argvalue))
697 
698  def unprotected_evaluate(self, da):
699  graph = self.get_minimum_spanning_tree()
700  score = 0.0
701  for e in graph.edges():
702  dist = graph.get_edge_data(*e)['weight']
703  prob = self.sigmoid(dist)
704  score += -numpy.log(prob)
705  score += self.linear_slope * dist
706  return score
707 
708  def do_get_inputs(self):
709  return self.particle_list
710 
711 
713 
714  '''
715 
716  '''
717 
718  def get_from_selection_tuple(self,tuples):
719  particles = []
720  for s in tuples:
721  ps = IMP.pmi.tools.select_by_tuple(
722  self.representation, s,
723  resolution=self.resolution, name_is_ambiguous=True)
724  particles += ps
725  return particles
726 
727  def __init__(
728  self,
729  representation=None,
730  objects_above=None,
731  objects_inside=None,
732  objects_below=None,
733  z_init=0.0,
734  z_min=0.0,
735  z_max=0.0,
736  thickness=30,
737  resolution=1,
738  label="None"):
739  self.resolution=resolution
740 
741  self.weight = 1.0
742  self.label = label
743  self.representation = representation
744  self.thickness = thickness
745  self.z_init=z_init
746 
747 
748  softness = 3.0
749  plateau = 1e-10
750  linear = 0.02
751 
752  mr=None
753  hierarchies_above = []
754  hierarchies_inside = []
755  hierarchies_below = []
756 
757  if representation is None:
758  hierarchies_above, hierarchies_inside, hierarchies_below = (
760  resolution,
761  flatten=True) for objects in [objects_above, objects_inside, objects_below])
762  for h in hierarchies_above, hierarchies_inside, hierarchies_below:
763  if h is not None:
764  self.m=hierarchies_above[0].get_model()
765  break
766  else:
767  self.m = representation.prot.get_model()
768  if objects_above is not None:
769  hierarchies_above = self.get_from_selection_tuple(objects_above)
770 
771  if objects_below is not None:
772  particles_below = self.get_from_selection_tuple(objects_below)
773 
774  if objects_inside is not None:
775  particles_inside = self.get_from_selection_tuple(objects_inside)
776 
777  self.z_center = IMP.pmi.tools.SetupNuisance(
778  self.m, self.z_init, z_min, z_max, isoptimized=True).get_particle()
780  self.m, self.z_center.get_particle_index(), self.thickness, softness, plateau, linear)
781  mr.add_particles_inside([h.get_particle().get_index()
782  for h in hierarchies_inside])
783  mr.add_particles_above([h.get_particle().get_index()
784  for h in hierarchies_above])
785  mr.add_particles_below([h.get_particle().get_index()
786  for h in hierarchies_below])
787 
788  self.rs = IMP.RestraintSet(self.m, label)
789  self.rs.add_restraint(mr)
790 
791  def create_box(self, x_center, y_center, hierarchy):
792 
793  z = self.z_center.get_nuisance()
794  p = IMP.Particle(self.m)
796  h.set_name("Membrane_" + self.label)
797  hierarchy.add_child(h)
798 
799  particles_box = []
800  p_origin = IMP.Particle(self.m)
801 
802  IMP.atom.Mass.setup_particle(p_origin, 100)
803  d = IMP.core.XYZR.setup_particle(p_origin)
804  d.set_coordinates((x_center, y_center, 0))
805  d.set_radius(1)
806  h_origin = IMP.atom.Hierarchy.setup_particle(p_origin)
807  h.add_child(h_origin)
808  particles_box.append(p_origin)
809 
810  p1 = IMP.Particle(self.m)
813  d.set_coordinates((x_center, y_center, z + 100))
814  d.set_radius(1)
816  h.add_child(h1)
817  particles_box.append(p1)
818 
819  p2 = IMP.Particle(self.m)
822  d.set_coordinates((x_center, y_center, z - 100))
823  d.set_radius(1)
825  h.add_child(h2)
826  particles_box.append(p2)
827 
830 
832 
833  p1.set_name("z_positive")
834  p2.set_name("z_negative")
835 
836  # to display the membrane in the rmf file
837  for offs in [self.thickness / 2, -self.thickness / 2]:
838  p1 = IMP.Particle(self.m)
841  d.set_coordinates((-100 + x_center, -100 + y_center, z + offs))
842  d.set_radius(1)
844  h.add_child(h1)
845  particles_box.append(p1)
846 
847  p2 = IMP.Particle(self.m)
850  d.set_coordinates((-100 + x_center, 100 + y_center, z + offs))
851  d.set_radius(1)
853  h.add_child(h2)
854  particles_box.append(p2)
855 
856  p3 = IMP.Particle(self.m)
859  d.set_coordinates((100 + x_center, -100 + y_center, z + offs))
860  d.set_radius(1)
862  h.add_child(h3)
863  particles_box.append(p3)
864 
865  p4 = IMP.Particle(self.m)
868  d.set_coordinates((100 + x_center, 100 + y_center, z + offs))
869  d.set_radius(1)
871  h.add_child(h4)
872  particles_box.append(p4)
873 
878 
883 
884  #sm = self._MembraneSingletonModifier(p_origin, self.z_center)
886  for p in particles_box:
887  smp = self._MembraneSingletonModifier(IMP.core.XYZ(p).get_z(),self.z_init, self.z_center)
889  IMP.core.XYZ(p).set_coordinates_are_optimized(True)
890  lcp.add(p.get_index())
891  c = IMP.container.SingletonsConstraint(smp, None, lcp)
892  self.m.add_score_state(c)
893  self.m.update()
894 
895  class _MembraneSingletonModifier(IMP.SingletonModifier):
896 
897  """A class that updates the membrane particles
898  """
899 
900  def __init__(self,z, z_init, z_nuisance):
901  IMP.SingletonModifier.__init__(
902  self, "MembraneSingletonModifier%1%")
903  self.z=z
904  self.z_init = z_init
905  self.z_nuisance = z_nuisance
906 
907  def apply_index(self, m, pi):
908  new_z = self.z_nuisance.get_nuisance()
909  d = IMP.core.XYZ(m, pi)
910  d.set_z(self.z + new_z-self.z_init)
911 
912 
913 
914  def do_get_inputs(self, m, pis):
915  return IMP.get_particles(m, pis)
916 
917  def do_get_outputs(self, m, pis):
918  return self.do_get_inputs(m, pis)
919 
920  def set_label(self, label):
921  self.label = label
922  self.rs.set_name(label)
923  for r in self.rs.get_restraints():
924  r.set_name(label)
925 
926  def add_to_model(self):
927  IMP.pmi.tools.add_restraint_to_model(self.m, self.rs)
928 
929  def get_restraint(self):
930  return self.rs
931 
932  def set_weight(self, weight):
933  self.weight = weight
934  self.rs.set_weight(weight)
935 
936  def get_particles_to_sample(self):
937  ps = {}
938 
939  ps["Nuisances_MembraneRestraint_Z_" +
940  self.label] = ([self.z_center], 2.0)
941  return ps
942 
943  def get_movers(self):
944  movers=[]
945  mover_name="Nuisances_MembraneRestraint_Z_" + self.label
946  particle=self.z_center
947  maxstep=2.0
948  mv=IMP.core.NormalMover([particle],
949  IMP.FloatKeys([IMP.FloatKey("nuisance")]),maxstep)
950  mv.set_name(mover_name)
951  movers.append(mv)
952 
953  return movers
954 
955 
956  def get_output(self):
957  self.m.update()
958  output = {}
959  score = self.weight * self.rs.unprotected_evaluate(None)
960  output["_TotalScore"] = str(score)
961  output["MembraneRestraint_" + self.label] = str(score)
962  output["MembraneRestraint_Z_" +
963  self.label] = str(self.z_center.get_nuisance())
964  return output
965 
966 
967 class FuzzyBoolean(object):
968 
969  '''
970  Fully Ambiguous Restraint that can be built using boolean logic
971  R. Pellarin. Pasteur Institute.
972  '''
973 
974  def __init__(self, p1, operator=None, p2=None):
975  '''
976  input a list of particles, the slope and theta of the sigmoid potential
977  theta is the cutoff distance for a protein-protein contact
978  '''
979  if isinstance(p1, FuzzyBoolean) and isinstance(p2, FuzzyBoolean):
980  self.operations = [p1, operator, p2]
981  self.value = None
982  else:
983  self.operations = []
984  self.value = p1
985 
986  def __or__(self, FuzzyBoolean2):
987  return FuzzyBoolean(self, self.or_, FuzzyBoolean2)
988 
989  def __and__(self, FuzzyBoolean2):
990  return FuzzyBoolean(self, self.and_, FuzzyBoolean2)
991 
992  def and_(self, a, b):
993  return a * b
994 
995  def or_(self, a, b):
996  return 1.0 - (1.0 - a) * (1.0 - b)
997 
998  def evaluate(self):
999 
1000  if len(self.operations) == 0:
1001  return self.value
1002  FuzzyBoolean1, op, FuzzyBoolean2 = self.operations
1003 
1004  return op(FuzzyBoolean1.evaluate(), FuzzyBoolean2.evaluate())
1005 
1006 
1008 
1009  '''
1010  Fully Ambiguous Restraint that can be built using boolean logic
1011  R. Pellarin. Pasteur Institute.
1012  '''
1013  plateau = 0.00000000000001
1014  theta = 5.0
1015  slope = 2.0
1016  innerslope = 0.01
1017 
1018  def __init__(self, m, p1, p2, operator=None):
1019  '''
1020  input a list of particles, the slope and theta of the sigmoid potential
1021  theta is the cutoff distance for a protein-protein contact
1022  '''
1023  IMP.Restraint.__init__(self, m, "FuzzyRestraint %1%")
1024  self.m = m
1025  self.min = sys.float_info.min
1026  if isinstance(p1, FuzzyRestraint) and isinstance(p2, FuzzyRestraint):
1027  self.operations = [p1, operator, p2]
1028  self.particle_pair = None
1029  elif isinstance(p1, FuzzyRestraint) and p2 is None:
1030  self.operations = [p1, operator, None]
1031  self.particle_pair = None
1032  else:
1033  self.operations = []
1034  self.particle_pair = (p1, p2)
1035 
1036  def __or__(self, FuzzyRestraint2):
1037  return FuzzyRestraint(self.m, self, FuzzyRestraint2, self.or_)
1038 
1039  def __and__(self, FuzzyRestraint2):
1040  return FuzzyRestraint(self.m, self, FuzzyRestraint2, self.and_)
1041 
1042  def __invert__(self):
1043  return FuzzyRestraint(self.m, self, None, self.invert_)
1044 
1045  def and_(self, a, b):
1046  c = a + b
1047  return c
1048 
1049  def or_(self, a, b):
1050  c = math.exp(-a) + math.exp(-b) - math.exp(-a - b)
1051  return -math.log(c)
1052 
1053  def invert_(self, a):
1054  c = 1.0 - math.exp(-a)
1055  return -math.log(c)
1056 
1057  def evaluate(self):
1058  if len(self.operations) == 0:
1059  return self.distance()
1060  FuzzyRestraint1, op, FuzzyRestraint2 = self.operations
1061 
1062  if FuzzyRestraint2 is not None:
1063  return op(FuzzyRestraint1.evaluate(), FuzzyRestraint2.evaluate())
1064  else:
1065  return op(FuzzyRestraint1.evaluate())
1066 
1067  def distance(self):
1068  d1 = IMP.core.XYZ(self.particle_pair[0])
1069  d2 = IMP.core.XYZ(self.particle_pair[1])
1070  d = IMP.core.get_distance(d1, d2)
1071  argvalue = (d-self.theta)/self.slope
1072  return -math.log(1.0 -(1.0-self.plateau)/(1.0+math.exp(-argvalue)))+self.innerslope*d
1073 
1074  def add_to_model(self):
1076 
1077  def unprotected_evaluate(self, da):
1078  return self.evaluate()
1079 
1080  def __str__(self):
1081  if len(self.operations) == 0:
1082  return str(self.particle_pair)
1083  FuzzyRestraint1, op, FuzzyRestraint2 = self.operations
1084  if FuzzyRestraint2 is not None:
1085  return str(FuzzyRestraint1) +str(op)+str(FuzzyRestraint2)
1086  else:
1087  return str(FuzzyRestraint1) +str(op)
1088 
1089  def do_get_inputs(self):
1090  if len(self.operations) == 0:
1091  return list(self.particle_pair)
1092  FuzzyRestraint1, op, FuzzyRestraint2 = self.operations
1093  if FuzzyRestraint2 is not None:
1094  return list(set(FuzzyRestraint1.do_get_inputs() +FuzzyRestraint2.do_get_inputs()))
1095  else:
1096  return list(set(FuzzyRestraint1.do_get_inputs()))
A base class for modifiers of ParticlesTemp.
A function that is harmonic over an interval.
Apply a SingletonFunction to a SingletonContainer to maintain an invariant.
this restraint allows ambiguous crosslinking between multiple copies excluding between symmetric copi...
void add_particles(RMF::FileHandle fh, const ParticlesTemp &hs)
Various classes to hold sets of particles.
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
A decorator for a particle which has bonds.
generate a connectivity restraint between domains setting up the restraint example: cr=restraints...
def get_full_graph
get the full graph of distances between every particle pair
Miscellaneous utilities.
Definition: tools.py:1
Add symmetric attribute to a particle.
Definition: Symmetric.h:23
ParticlesTemp get_particles(Model *m, const ParticleIndexes &ps)
A score on the distance between the surfaces of two spheres.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
generates and wraps a IMP.pmi.ConnectivityRestraint between domains example: cr=restraints.ConnectivityNetworkRestraint(simo,["CCC",(1,100,"TTT"),(100,150,"AAA")]) cr.add_to_model() cr.set_label("CR1")
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
Definition: RestraintSet.h:36
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:1579
a python restraint that computes the score for a composite of proteins Authors: G.
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def __init__
input a list of particles, the slope and theta of the sigmoid potential theta is the cutoff distance ...
Store a list of ParticleIndexes.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Definition: Mass.h:44
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:38
static Bonded setup_particle(Model *m, ParticleIndex pi)
A restraint for ambiguous cross-linking MS data and multiple state approach.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
handleparticles is a selection tuple compositeparticles is a list of selection tuples ...
Modify a set of continuous variables using a normal distribution.
Definition: NormalMover.h:20
Classes for writing output files and processing them.
Definition: output.py:1
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
Fully Ambiguous Restraint that can be built using boolean logic R.
def __init__
input a list of particles, the slope and theta of the sigmoid potential theta is the cutoff distance ...
VectorD< 3 > Vector3D
Definition: VectorD.h:395
Class to handle individual particles of a Model object.
Definition: Particle.h:41
Floats get_list_of_bipartite_minimum_sphere_distance(const ParticlesTemps &pss)
Definition: utilities.h:71
def sigmoid
a sigmoid function that scores the probability of a contact between two proteins
Restraint * create_connectivity_restraint(const Selections &s, double x0, double k, std::string name="Connectivity%1%")
Create a restraint connecting the selections.
Applies a PairScore to a Pair.
Definition: PairRestraint.h:29
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:700
Python classes to represent, score, sample and analyze models.
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
Fully Ambiguous Restraint that can be built using boolean logic R.
def __init__
input a list of particles, the slope and theta of the sigmoid potential theta is the cutoff distance ...
virtual ModelObjectsTemp do_get_inputs() const =0
A restraint is a term in an IMP ScoringFunction.
Definition: Restraint.h:52
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27