IMP logo
IMP Reference Guide  2.10.0
The Integrative Modeling Platform
domino_model.py
1 """@namespace IMP.EMageFit.domino_model
2  Classes to manage a model using DOMINO.
3 """
4 
5 import sys
6 import time
7 import csv
8 import logging
9 log = logging.getLogger("DominoModel")
10 
11 import IMP
12 import IMP.domino as domino
13 import IMP.core as core
14 import IMP.container as container
15 import IMP.display as display
16 import IMP.atom as atom
17 import IMP.algebra as alg
18 import IMP.em2d as em2d
19 import IMP.multifit as multifit
20 
21 import IMP.EMageFit.imp_general.comparisons as comparisons
22 import IMP.EMageFit.imp_general.alignments as alignments
23 import IMP.EMageFit.imp_general.io as io
24 import IMP.EMageFit.imp_general.representation as representation
25 
26 import IMP.EMageFit.restraints as restraints
27 import IMP.EMageFit.buildxlinks as buildxlinks
28 import IMP.EMageFit.solutions_io as solutions_io
29 import IMP.EMageFit.csv_related as csv_related
30 
31 
32 field_delim = "," # separate fields in the output text file
33 unit_delim = "/" # separate units within a field (eg, reference frames).
34  # It is used in output text files and in databse files
35 
36 
38 
39  """
40  Management of a model using DOMINO
41  """
42 
43  def __init__(self, name="my model"):
44  self.model = IMP.Model()
45  self.model.set_name(name)
46  self.configuration_sampling_done = False
47  self.assignments_sampling_done = False
48  self.rb_states_table = domino.ParticleStatesTable()
49  self.merge_tree = None
50  # measure_models is True only when doing benchmark
51  self.measure_models = False
52  # The first time that create_coarse is found and its value is True
53  # a coarse version of the assembly is built
54  self.create_coarse = True
55  self.restraints = dict()
56  self.xlinks_dict = buildxlinks.XlinksDict()
57 
58  def add_restraint(self, r, name, weight, max_score=False):
59  """
60  Adds a restraint to the model
61  @param r An IMP.Restraint object
62  @param name Name for the restraint
63  @param weight Weight for the restraint
64  @param max_score Maximum score allowed for the restraint. If
65  max_score is False, there is no limit for the value of the restraint
66  """
67  log.info("Adding restraint %s weight %s max_score %s",
68  name, weight, max_score)
69  if "-" in name:
70  raise ValueError("SQL database does not accept restraint names "
71  "containing -")
72  r.set_name(name)
73  if max_score is not None and max_score:
74  r.set_maximum_score(max_score)
75  r.set_weight(weight)
76  self.restraints[name] = r
77 
79  """
80  Aligns the set of structures considered for DOMINO sampling.
81  The function:
82  1) reads the possible reference frames that the
83  rb_states_table stores for each rigid body. This table
84  must be filled before using this function. Usually it is
85  filled with the results from a previous Monte Carlo sampling.
86  If the solutions from Monte Carlo seem to have the same structure
87  but they are not aligned to each other, this function can
88  help setting better starting positions to use with DOMINO.
89  2) Gets the first state for each of the rigid bodies and sets
90  a reference structure using such states.
91  3) Aligns all the rest of the structures respect to the reference
92  4) Replaces the values of the reference frames stored in the
93  rb_states_table with the new values obtained from the alignments.
94  It does it for all states of a rigid body.
95  @note: If this function is applied, the parameters "anchor" and "fixed"
96  are ignored, as they are superseded by the use of the alignments
97  calculated here.
98  """
99  log.debug("Align the configurations read from the database before "
100  "running DOMINO")
101  rb_states = []
102  n_states = []
103  for rb in self.components_rbs:
104  ps = self.rb_states_table.get_particle_states(rb)
105  n = ps.get_number_of_particle_states()
106  if n == 0:
107  raise ValueError(
108  "There are no particle states for %s" %
109  rb.get_name())
110  n_states.append(n)
111  rb_states.append(ps)
112  # coordinates of the first configuration (there is at least one for all
113  # the rbs)
114  for rb, states in zip(self.components_rbs, rb_states):
115  states.load_particle_state(0, rb)
116  reference_coords = get_coordinates(self.components_rbs)
117  aligned_transformations = [[] for rb in self.components_rbs]
118  for i in range(1, max(n_states)):
119  log.debug("Aligning rigid bodies configuration %s" % i)
120  # load the configuration
121  for j in range(len(self.components_rbs)):
122  # if there are no more particle states for this rigid body, use
123  # the last
124  state_index = min(i, n_states[j] - 1)
125  rb_states[j].load_particle_state(
126  state_index,
127  self.components_rbs[j])
128  coords = get_coordinates(self.components_rbs)
129  T = alg.get_transformation_aligning_first_to_second(
130  coords, reference_coords)
131  for j, rb in enumerate(self.components_rbs):
132  t = rb.get_reference_frame().get_transformation_to()
133  new_t = alg.compose(T, t)
134  aligned_transformations[j].append(new_t)
135  # set the new possible transformations
136  for i, rb in enumerate(self.components_rbs):
137  self.set_rb_states(rb, aligned_transformations[i])
138 
139  def set_xlink_restraint(self, id1, chain1, residue1, id2, chain2, residue2,
140  distance, weight, stddev, max_score=False):
141  """
142  Set a restraint on the maximum distance between 2 residues
143  @param id1 Name of the first component
144  @param chain1
145  @param residue1 Residue number for the aminoacid in the first
146  component.The number is the number in the PDB file, not the
147  number relative to the beginning of the chain
148  @param id2 Name of the second component
149  @param chain2
150  @param residue2 Residue number for the aminoacid in the second
151  component.
152  @param distance Maximum distance tolerated
153  @param weight Weight of the restraint
154  @param stddev Standard deviation used to approximate the
155  HarmonicUpperBound function to a Gaussian
156  @param max_score See help for add_restraint(). If none is given,
157  the maximum score is set to allow a maximum distance of 10 Angstrom
158  greater than the parameter "distance".
159  """
160  xlink = buildxlinks.Xlink(
161  id1,
162  chain1,
163  residue1,
164  id2,
165  chain2,
166  residue2,
167  distance)
168  log.info("Setting cross-linking restraint ")
169  log.info("%s", xlink.show())
170  self.xlinks_dict.add(xlink)
171  # setup restraint
172  A = representation.get_component(self.assembly, xlink.first_id)
173  s1 = IMP.atom.Selection(A, chain=xlink.first_chain,
174  residue_index=xlink.first_residue)
175  p1 = s1.get_selected_particles()[0]
176  B = representation.get_component(self.assembly, xlink.second_id)
177  s2 = IMP.atom.Selection(B, chain=xlink.second_chain,
178  residue_index=xlink.second_residue)
179  p2 = s2.get_selected_particles()[0]
180  k = core.Harmonic.get_k_from_standard_deviation(stddev)
181  score = core.HarmonicUpperBound(xlink.distance, k)
182  pair_score = IMP.core.DistancePairScore(score)
183  r = IMP.core.PairRestraint(self.model, pair_score, (p1, p2))
184  if not max_score:
185  error_distance_allowed = 100
186  max_score = weight * \
187  score.evaluate(distance + error_distance_allowed)
188  log.info("%s, max_score %s", xlink.show(), max_score)
189  self.add_restraint(r, xlink.get_name(), weight, max_score)
190 
191  def set_complementarity_restraint(self, name1, name2, rname,
192  max_sep_distance, max_penetration,
193  weight, max_score=1e10):
194  """
195  Set a restraint for geometric complementarity between 2 components
196  @param name1 name of
197  @param name2 - The restraint is applied to this components
198  @param rname - The name given to the restraint
199  @param max_sep_distance - maximum distance between molecules
200  tolerated by the restraint
201  @param max_penetration - Maximum penetration allowed (angstrom)
202  @param weight
203  @param max_score
204  """
205  log.info("Setting geometric complementarity restraint %s: %s - %s",
206  rname, name1, name2)
207  A = representation.get_component(self.assembly, name1)
208  B = representation.get_component(self.assembly, name2)
209  restraint = multifit.ComplementarityRestraint(atom.get_leaves(A),
210  atom.get_leaves(B), rname)
211  restraint.set_maximum_separation(max_sep_distance)
212  restraint.set_maximum_penetration(max_penetration)
213  log.debug("Maximum separation: %s Maximum penetration score: %s",
214  max_sep_distance, max_penetration)
215  self.add_restraint(restraint, rname, weight, max_score)
216 
217  def create_coarse_assembly(self, n_residues):
218  """
219  Simplify the assembly with a coarse representation
220  @param n_residues Number of residues used for a coarse particle
221  """
222  if self.create_coarse:
223  log.info("Creating simplified assembly")
224  self.coarse_assembly = \
225  representation.create_simplified_assembly(self.assembly,
226  self.components_rbs, n_residues)
227  self.create_coarse = False
228 
229  def do_sampling(self, mode="assignments_heap_container", params=None):
230  """
231  Do Domino sampling
232  @param mode The mode used to recover solutions from domino.
233  It can have
234  the value "configuration", "assignments", or
235  "assignments_heap_container". The mode
236  "configuration" recovers all possible configurations
237  (takes a while to generate them).
238  The mode "assignments" only provides the assignments
239  (state numbers) of the solutions. It is faster but requires
240  generating the solutions afterwards
241  The mode "assignments_heap_container" selects the best solutions
242  after each merging in DOMINO, discarding the rest.
243  In practice I used the mode "assignments_heap_container"
244  @param params
245  """
246  t0 = time.time()
247  if mode == "configuration":
248  self.configuration_set = self.sampler.create_sample()
249  tf = time.time()
250  log.info("found %s configurations. Time %s",
251  self.configuration_set.get_number_of_configurations(),
252  tf - t0)
253  self.configuration_sampling_done = True
254  elif mode == "assignments":
255  subset = self.rb_states_table.get_subset()
256  log.info(
257  "Do sampling with assignments. Subset has %s elements: %s",
258  len(subset), subset)
259  self.solution_assignments = \
260  self.sampler.get_sample_assignments(subset)
261  tf = time.time()
262  log.info("found %s assignments. Time %s",
263  len(self.solution_assignments), tf - t0)
264  self.assignments_sampling_done = True
265  elif mode == "assignments_heap_container":
266  subset = self.rb_states_table.get_subset()
267  log.info("DOMINO sampling an assignments_heap_container. "
268  "Subset has %s elements: %s", len(subset), subset)
269  # last vertex is the root of the merge tree
270  root = self.merge_tree.get_vertices()[-1]
271  container = self.get_assignments_with_heap(root,
272  params.heap_solutions)
273  self.solution_assignments = container.get_assignments()
274  tf = time.time()
275  log.info("found %s assignments. Time %s",
276  len(self.solution_assignments), tf - t0)
277  self.assignments_sampling_done = True
278  else:
279  raise ValueError("Requested wrong mode for the DOMINO sampling")
280 
281  def get_assignment_text(self, assignment, subset):
282  """
283  Get the formatted text for an assignment
284  @param subset Subset of components of the assembly
285  @param assignment Assignment class with the states for the subset
286  @note: The order in assignment and subset is not always the order
287  of the rigid bodies in self.components_rbs. The function reorders
288  the terms in the assignment so there is correspondence.
289  """
290  ordered_assignment = []
291  for rb in self.components_rbs:
292  for i, particle in enumerate(subset):
293  if rb.get_particle() == particle:
294  ordered_assignment.append(assignment[i])
295  text = "|".join([str(i) for i in ordered_assignment])
296  return text
297 
298  def get_assignment_and_RFs(self, assignment, subset):
299  """
300  Return a line with texts for an assignment and the
301  the reference frames of the RigidBodies in the subset.
302  @param subset Subset of components of the assembly
303  @param assignment Assignment class with the states for the subset
304  @note: see the get_assignment_assignment_text() note.
305  """
306  ordered_assignment = []
307  RFs = []
308  for rb in self.components_rbs:
309  for i, particle in enumerate(subset):
310  if rb.get_particle() == particle:
311  j = assignment[i]
312  ordered_assignment.append(j)
313  pstates = self.rb_states_table.get_particle_states(rb)
314  pstates.load_particle_state(j, rb.get_particle())
315  RFs.append(rb.get_reference_frame())
316  RFs_text = unit_delim.join(
317  [io.ReferenceFrameToText(ref).get_text() for ref in RFs])
318  assignment_numbers = "|".join([str(i) for i in ordered_assignment])
319  return [assignment_numbers, RFs_text]
320 
321  def get_assignment_info(self, assignment, subset):
322  """
323  Info related to an assignment. Returns a list with text for the
324  assignment and all the scores for the restraints
325  @param subset Subset of components of the assembly
326  @param assignment Assignment class with the states for the subset
327  """
328  text = self.get_assignment_and_RFs(assignment, subset)
329  scores, total_score = self.get_assignment_scores(assignment, subset)
330  info = text + [total_score] + scores
331  return info
332 
333  def get_assignment_scores(self, assignment, subset):
334  """
335  Returns a list with the values for the restraints on a subset of
336  the components of the assembly
337  @param subset Subset of components of the assembly
338  @param assignment Assignment class with the states for the subset
339  """
340  restraints_to_evaluate = \
341  self.restraint_cache.get_restraints(subset, [])
342  scores = [self.restraint_cache.get_score(r, subset, assignment)
343  for r in restraints_to_evaluate]
344  total_score = sum(scores)
345  return scores, total_score
346 
347  def load_state(self, assignment):
348  """
349  Load the positions of the components given by the assignment
350  @param assignment Assignment class
351  """
352  subset = self.rb_states_table.get_subset()
353  domino.load_particle_states(subset, assignment, self.rb_states_table)
354 
355  def get_assignments_with_heap(self, vertex, k=0):
356  """
357  Domino sampling that recovers the assignments for the root of the
358  merge tree, but
359  conserving only the best k scores for each vertex of the tree.
360  @param[in] vertex Vertex with the root of the current merge tree. This
361  function is recursive.
362  @param[in] k
363  """
364  if(self.sampler.get_number_of_subset_filter_tables() == 0):
365  raise ValueError("No subset filter tables")
366  if(self.merge_tree is None):
367  raise ValueError("No merge tree")
368  # In the merge tree, the names of the vertices are the subsets.
369  # The type of the vertex name is a domino.Subset
370  subset = self.merge_tree.get_vertex_name(vertex)
371  log.info("computing assignments for vertex %s", subset)
372  t0 = time.time()
373  assignments_container = domino.HeapAssignmentContainer(subset, k,
374  self.restraint_cache, "my_heap_assignments_container")
375  neighbors = self.merge_tree.get_out_neighbors(vertex)
376  if len(neighbors) == 0: # A leaf
377  # Fill the container with the assignments for the leaf
378  self.sampler.load_vertex_assignments(vertex, assignments_container)
379  else:
380  if neighbors[0] > neighbors[1]:
381  # get_vertex_assignments() methods in domino
382  # expects the children in sorted order
383  neighbors = [neighbors[1], neighbors[0]]
384  # recurse on the two children
385  container_child0 = self.get_assignments_with_heap(neighbors[0], k)
386  container_child1 = self.get_assignments_with_heap(neighbors[1], k)
387  self.sampler.load_vertex_assignments(vertex,
388  container_child0,
389  container_child1,
390  assignments_container)
391  tf = time.time() - t0
392  log.info("Merge tree vertex: %s assignments: %s Time %s sec.", subset,
393  assignments_container.get_number_of_assignments(), tf)
394  return assignments_container
395 
396  def get_restraint_value_for_assignment(self, assignment, name):
397  """
398  Recover the value of a restraint
399  @param name of the restraint
400  @param assignment Assignment class containing the assignment for
401  the solution desired
402  """
403  if not self.assignments_sampling_done:
404  raise ValueError("DominoModel: sampling not done")
405  log.debug("Computing restraint value for assignment %s", assignment)
406  self.load_state(assignment)
407  for rb in self.components_rbs:
408  log.debug("rb - %s",
409  rb.get_reference_frame().get_transformation_to())
410  val = self.restraints[name].evaluate(False)
411  log.debug("restraint %s = %s", restraint.get_name(), val)
412  return val
413 
415  """
416  Sets the native model for benchmark, by reading the native
417  structure and set the rigid bodies.
418  """
419  self.measure_models = True
420  self.native_model = IMP.Model()
421  if hasattr(params.benchmark, "fn_pdb_native"):
422  self.native_assembly = \
423  representation.create_assembly_from_pdb(self.native_model,
424  params.benchmark.fn_pdb_native, params.names)
425  elif hasattr(params.benchmark, "fn_pdbs_native"):
426  self.native_assembly = \
427  representation.create_assembly(self.native_model,
428  params.benchmark.fn_pdbs_native, params.names)
429  else:
430  raise ValueError("set_native_assembly_for_benchmark: Requested "
431  "to set a native assembly for benchmark but did not provide "
432  "its pdb, or the pdbs of the components")
433  self.native_rbs = representation.create_rigid_bodies(
434  self.native_assembly)
435  anchor_assembly(self.native_rbs, params.anchor)
436 
437  def set_rb_states(self, rb, transformations):
438  """
439  Add a set of reference frames as possible states for a rigid body
440  @param rb The rigid body
441  @param transformations Transformations are used to build the
442  reference frames for the rigid body.
443  """
444  log.info("Set rigid body states for %s", rb.get_name())
445  RFs = [alg.ReferenceFrame3D(T) for T in transformations]
446  for st in RFs:
447  log.debug("%s", st)
448  rb_states = domino.RigidBodyStates(RFs)
449  self.rb_states_table.set_particle_states(rb, rb_states)
450  st = self.rb_states_table.get_particle_states(rb)
451  log.info("RigidBodyStates created %s",
452  st.get_number_of_particle_states())
453 
455  """
456  Add a RestraintScoreSubsetFilterTable to DOMINO that allows to
457  reject assignments that have score worse than the thresholds for
458  the restraints
459  """
460  log.info("Adding RestraintScoreSubsetFilterTable")
461  # Restraint Cache is a centralized container of the restraints and
462  # their scores
463  self.restraint_cache = domino.RestraintCache(self.rb_states_table)
464  self.restraint_cache.add_restraints(self.restraints.values())
465  rssft = domino.RestraintScoreSubsetFilterTable(self.restraint_cache)
466  rssft.set_name('myRestraintScoreSubsetFilterTable')
467  self.sampler.add_subset_filter_table(rssft)
468 
469  def set_assembly_components(self, fn_pdbs, names):
470  """
471  Sets the components of an assembly, each one given as a separate
472  PDB file.
473  @param fn_pdbs List with the names of the pdb files
474  @param names List with the names of the components of the assembly.
475  """
476  if len(fn_pdbs) != len(names):
477  raise ValueError("Each PDB file needs a name")
478  self.names = names
479  self.assembly = representation.create_assembly(self.model, fn_pdbs,
480  names)
481  self.components_rbs = representation.create_rigid_bodies(self.assembly)
482  self.not_optimized_rbs = []
483 
484  def set_assembly(self, fn_pdb, names):
485  """
486  Sets an assembly from a single PDB file. The function assumes that
487  the components of the assembly are the chains of the PDB file.
488  @param fn_pdb PDB with the file for the assembly
489  @param names Names for each of the components (chains)
490  of the assembly
491  """
492  self.assembly = representation.create_assembly_from_pdb(self.model,
493  fn_pdb,
494  names)
495  self.components_rbs = representation.create_rigid_bodies(self.assembly)
496  self.not_optimized_rbs = []
497 
499  """
500  BranchAndBoundAssignmentsTable enumerate states based on provided
501  filtered using the provided the subset filter tables
502  """
503  log.info("Setting BranchAndBoundAssignmentsTable")
504  fts = []
505  for i in range(self.sampler.get_number_of_subset_filter_tables()):
506  ft = self.sampler.get_subset_filter_table(i)
507  fts.append(ft)
508  at = domino.BranchAndBoundAssignmentsTable(self.rb_states_table, fts)
509  self.sampler.set_assignments_table(at)
510  self.assignments_table = at
511 
513  """
514  ExclusionSubsetFilterTable ensures that two particles are not given
515  the same state at the same time
516  """
517  log.info("Setting ExclusionSubsetFilterTable")
518  ft = domino.ExclusionSubsetFilterTable(self.rb_states_table)
519  self.sampler.add_subset_filter_table(ft)
520 
522  """
523  Creates a DOMINO sampler and adds the required filter tables
524  """
525  if self.merge_tree is None:
526  raise ValueError("Merge tree for DOMINO not set. It is required "
527  "to setup the sampler")
528  log.info("Domino sampler")
529  self.sampler = domino.DominoSampler(self.model, self.rb_states_table)
530  self.sampler.set_restraints(self.restraints.values())
531  self.sampler.set_log_level(IMP.TERSE)
532  self.sampler.set_merge_tree(self.merge_tree)
536 
537  def read_merge_tree(self, fn):
538  """
539  Read and create a merge tree from a file.
540  The format of the file is the format written by write merge_tree()
541  It does not matter the order of the indices in the file, as
542  they are sorted by this function.
543  The only condition is that the index for the vertex that is the
544  root of the tree MUST be the largest. This is guaranteed when
545  creating a merge tree with create_merge_tree()
546  @param fn File containing the merge tree
547  """
548  log.info("Reading merge tree from %s", fn)
549  ps = self.rb_states_table.get_particles()
550  log.debug("particles")
551  for p in ps:
552  log.debug("%s", p.get_name())
553  rows = csv_related.read_csv(fn, delimiter=field_delim)
554  for i in range(len(rows)):
555  row = rows[i]
556  rows[i] = [int(row[0]), int(row[1]), int(row[2]), row[3]]
557  # Sort rows by vertice index
558  rows.sort()
559  subsets = []
560  # build subsets from the rows text
561  for row in rows:
562  names = row[3].split(unit_delim)
563  log.debug("row %s names %s", row, names)
564  # get particles in the subset
565  particles = []
566  for name in names:
567  l = [p for p in ps if p.get_name() == name]
568  if (len(l) > 1):
569  ValueError("More than one particle with same name" % names)
570  particles.extend(l)
571  s = domino.Subset(particles)
572  subset_names = [p.get_name() for p in particles]
573  log.debug("Merge tree Subset %s. Particles %s ", s, subset_names)
574  subsets.append(s)
575  # The vertex with the largest index is the root.
576  # trick: get the merge tree object from a tree with only one node ...
577  jt = domino.SubsetGraph()
578  jt.add_vertex(subsets[0])
579  mt = domino.get_merge_tree(jt)
580  # ... and add the rest of the vertices
581  for i in range(1, len(subsets)):
582  mt.add_vertex(subsets[i]) # the name of the vertex is a subset
583  # set edges
584  for row in rows:
585  vid = row[0]
586  child_left = row[1]
587  child_right = row[2]
588  if child_left != -1:
589  mt.add_edge(vid, child_left)
590  if child_right != -1:
591  mt.add_edge(vid, child_right)
592  self.merge_tree = mt
593  log.info("%s", self.merge_tree.show_graphviz())
594 
595  def create_merge_tree(self):
596  """
597  Creates a merge tree from the restraints that are set already
598  """
599  rs = self.restraints.values()
600  ig = domino.get_interaction_graph(rs, self.rb_states_table)
601 # pruned_dep = IMP.get_pruned_dependency_graph(self.model)
602 # IMP.show_graphviz(pruned_dep)
603 # IMP.show_graphviz(ig)
604  jt = domino.get_junction_tree(ig)
605 # IMP.show_graphviz(jt)
606  self.merge_tree = domino.get_balanced_merge_tree(jt)
607 # IMP.show_graphviz(self.merge_tree)
608  log.info("Balanced merge tree created")
609  log.info("%s", self.merge_tree.show_graphviz())
610 
611  def write_merge_tree(self, fn):
612  """
613  Writes the current merge tree to file. The format is:
614  parent | child_left | child_right | labels
615  @param fn File for storing the merge tree
616  """
617  log.info("Writing merge to %s", fn)
618  if self.merge_tree is None:
619  raise ValueError("No merge tree")
620  f = open(fn, "w")
621  f.write(
622  "# Vertex index | Child0 | Child1 | label (names of the particles)\n")
623  w = csv.writer(f, delimiter=field_delim)
624  root = self.merge_tree.get_vertices()[-1]
625  self.write_subset(root, w)
626  f.close()
627 
628  def write_subset(self, v, w):
629  """
630  Writes lines describing a subset (vertex of a merge tree)
631  @param w A csv.writer
632  @param v Vertex index
633  """
634  no_child_index = -1
635  subset = self.merge_tree.get_vertex_name(v)
636  names = [p.get_name() for p in subset]
637  reps = [x for x in names if names.count(x) > 1] # repeated names
638  if(len(reps)):
639  raise ValueError("The names of the "
640  "particles in the subset %s are not unique" % subset)
641  names = unit_delim.join(names)
642  neighbors = self.merge_tree.get_out_neighbors(v)
643  if len(neighbors) == 0: # A leaf
644  # Fill the container with the assignments for the leaf
645  w.writerow([v, no_child_index, no_child_index, names])
646  else:
647  if neighbors[0] > neighbors[1]:
648  neighbors = [neighbors[1], neighbors[0]]
649  w.writerow([v, neighbors[0], neighbors[1], names])
650  self.write_subset(neighbors[0], w)
651  self.write_subset(neighbors[1], w)
652 
653  def set_connectivity_restraint(self, names, rname,
654  distance=10,
655  weight=1.0, max_score=None, n_pairs=1,
656  stddev=2):
657  """
658  Set a connectivity restraint between the elements of the assembly
659  @param names List with all the elements to be connected
660  @param distance Maximum distance tolerated by the restraint
661  @param rname Name for the restraint
662  @param weight Weight for the restraint
663  @param max_score Maximum score for the restraint
664  @param n_pairs Number of pairs of particles used in the
665  KClosePairScore of the connectivity restraint
666  @param stddev Standard deviation used to approximate the
667  HarmonicUpperBound function to a Gaussian
668  """
669  components = []
670  for name in names:
671  c = representation.get_component(self.coarse_assembly, name)
672  components.append(c)
673  log.info("Connectivity restraint for %s", components)
674  spring_constant = core.Harmonic.get_k_from_standard_deviation(stddev)
675  if max_score is None:
676  max_score = weight * spring_constant * n_pairs * (stddev) ** 2
677  r = restraints.get_connectivity_restraint(
678  components, distance, n_pairs,
679  spring_constant)
680  self.add_restraint(r, rname, weight, max_score)
681 
682  def set_em2d_restraint(self, name, fn_images_sel, pixel_size, resolution,
683  n_projections, weight, max_score=False,
684  mode="fast", n_optimized=1):
685  """
686  Set a Restraint computing the em2d score.
687  @param name Name for the restraint
688  @param fn_images_sel Selection file with the names of the images
689  used for the restraint
690  @param pixel_size - pixel size in the images
691  @param resolution - resolution used to downsample the projections
692  of the model when projecting
693  @param weight Weight of the restraint
694  @param max_score - Maximum value tolerated for the restraint
695  @param n_projections - Number of projections to generate
696  when projecting the model to do the coarse alignments
697  @param mode - Mode used for computing the restraints.
698  @param n_optimized - number of results from the coarse
699  alignment that are refined with the Simplex algorithm
700  to get a more accurate value for the restraint
701  """
702  log.info("Adding a em2D restraint: %s", fn_images_sel)
703  restraint_params = em2d.Em2DRestraintParameters(pixel_size,
704  resolution,
705  n_projections)
706  r = restraints.get_em2d_restraint(self.assembly,
707  fn_images_sel,
708  restraint_params,
709  mode, n_optimized)
710  self.add_restraint(r, name, weight, max_score)
711 
712  def set_not_optimized(self, name):
713  """
714  Set a part of the model as not optimized (it does not move during
715  the model optimization)
716  @param name of the component to optimized
717  """
718  if name not in self.names:
719  raise ValueError("DominoModel: There is not component "
720  "in the assembly with this name")
721  rb_name = representation.get_rb_name(name)
722  self.not_optimized_rbs.append(rb_name)
723 
724  def set_close_pairs_excluded_volume_restraint(self, distance=10,
725  weight=1.0, ratio=0.05, stddev=2,
726  max_score=False):
727  """
728  Creates an excluded volume restraint between all the components
729  of the coarse assembly.See help for
730  @param distance Maximum distance tolerated between particles
731  @param weight Weight for the restraint
732  @param stddev
733  @param max_score Maximum value for the restraint. If the parameter
734  is None, an automatic value is computed (using the ratio).
735  @param ratio Fraction of the number of possible pairs of
736  particles that are tolerated to overlap. The maximum
737  score is modified according to this ratio.
738  I have found that ratios of 0.05-0.1 work well allowing for
739  some interpenetration
740  @param stddev Standard deviation used to approximate the
741  HarmonicUpperBound function to a Gaussian
742  """
743  k = core.Harmonic.get_k_from_standard_deviation(stddev)
744  for i, c1 in enumerate(self.coarse_assembly.get_children()):
745  for j, c2 in enumerate(self.coarse_assembly.get_children()):
746  if j > i:
747  name = "exc_vol_%s_%s" % (c1.get_name(), c2.get_name())
748 
749  ls1 = atom.get_leaves(c1)
750  ls2 = atom.get_leaves(c2)
751  possible_pairs = len(ls1) * len(ls2)
752  n_pairs = possible_pairs * ratio
753 
754  marker1 = IMP.Particle(
755  self.model, "marker1 " + name)
756  marker2 = IMP.Particle(
757  self.model, "marker2 " + name)
758  table_refiner = core.TableRefiner()
759  table_refiner.add_particle(marker1, ls1)
760  table_refiner.add_particle(marker2, ls2)
761  score = core.HarmonicLowerBound(distance, k)
762  pair_score = core.SphereDistancePairScore(score)
763  close_pair_score = core.ClosePairsPairScore(pair_score,
764  table_refiner,
765  distance)
766  r = core.PairRestraint(self.model, close_pair_score,
767  (marker1, marker2))
768 
769  if not max_score:
770  minimum_distance_allowed = 0
771  max_score = weight * n_pairs * \
772  score.evaluate(minimum_distance_allowed)
773  log.debug("Setting close_pairs_excluded_volume_restraint(), "
774  "max_score %s", max_score)
775  self.add_restraint(r, name, weight, max_score)
776 
777  def set_pair_score_restraint(self, name1, name2,
778  restraint_name, distance=10,
779  weight=1.0, n_pairs=1, stddev=2, max_score=None):
780  """
781  Set a pair_score restraint between the coarse representation
782  of two components
783  @param name1 Name of the first component
784  @param name2 Name of the second component
785  @param restraint_name Name for the restraint
786  @param distance Maximum distance tolerated between particles
787  @param weight Weight of the restraint
788  @param n_pairs
789  @param max_score Maximum value tolerated for the restraint
790  @param stddev Standard deviation used to approximate the
791  HarmonicUpperBound function to a Gaussian
792  """
793  table_refiner = core.TableRefiner()
794 
795  # The particles in A are attached to a marker particle via a refiner.
796  # When the refiner gets a request for marker1, it returns the attached
797  # particles
798  A = representation.get_component(self.coarse_assembly, name1)
799  marker1 = IMP.Particle(self.model, "marker1 " + restraint_name)
800  table_refiner.add_particle(marker1, atom.get_leaves(A))
801  # same for B
802  B = representation.get_component(self.coarse_assembly, name2)
803  marker2 = IMP.Particle(self.model, "marker2 " + restraint_name)
804  table_refiner.add_particle(marker2, atom.get_leaves(B))
805 
806  k = core.Harmonic.get_k_from_standard_deviation(stddev)
808  # The score is set for the n_pairs closest pairs of particles
809  pair_score = core.KClosePairsPairScore(score, table_refiner, n_pairs)
810  # When KClosePairsPairScore is called, the refiner will provide the
811  # particles for A and B
812  if not max_score:
813  # Build a maximum score based on the function type that is used,
814  # an HarmonicUpperBound
815  temp_score = core.HarmonicUpperBound(distance, k)
816  error_distance_allowed = 10
817  max_score = weight * n_pairs * \
818  temp_score.evaluate(distance + error_distance_allowed)
819 
820  log.info("Setting pair score restraint for %s %s. k = %s, max_score "
821  "= %s, stddev %s", name1, name2, k, max_score, stddev)
822  r = core.PairRestraint(self.model, pair_score, (marker1, marker2))
823  self.add_restraint(r, restraint_name, weight, max_score)
824 
825  def write_solutions_database(self, fn_database, max_number=None):
826  """
827  Write the results of the DOMINO sampling to a SQLite database.
828  @param fn_database
829  @param max_number Maximum number of results to write
830  """
831  log.info("Creating the database of solutions")
832  if self.measure_models and not hasattr(self, "native_model"):
833  raise ValueError("The native model has not been set")
835  db.create(fn_database, overwrite=True)
836  db.connect(fn_database)
837  subset = self.rb_states_table.get_subset()
838  restraints_names = self.get_restraints_names_used(subset)
839  db.add_results_table(restraints_names, self.measure_models)
840  n = len(self.solution_assignments)
841  if max_number is not None:
842  n = min(n, max_number)
843  for sol_id, assignment in enumerate(self.solution_assignments[0:n]):
844  txt = self.get_assignment_text(assignment, subset)
845  # load the positions of the components in the state
846  self.load_state(assignment)
847  RFs = [rb.get_reference_frame() for rb in self.components_rbs]
848  scores, total_score = self.get_assignment_scores(
849  assignment, subset)
850  measures = None
851  if(self.measure_models):
852  measures = measure_model(self.assembly, self.native_assembly,
853  self.components_rbs, self.native_rbs)
854  db.add_record(sol_id, txt, RFs, total_score, scores, measures)
855  # Add native state if doing benchmarking
856  if self.measure_models:
857  assignment = "native"
858  RFs = [rb.get_reference_frame() for rb in self.native_rbs]
859  scores, total_score = self.get_restraints_for_native(
860  restraints_names)
861  db.add_native_record(assignment, RFs, total_score, scores)
862  db.save_records()
863  db.close()
864 
865  def get_restraints_names_used(self, subset):
866  """ Get the names of the restraints used for a subset
867  @param subset Subset of components of the assembly
868  """
869  return [r.get_name() for r
870  in self.restraint_cache.get_restraints(subset, [])]
871 
872  def get_restraints_for_native(self, restraints_names):
873  """
874  Get values of the restraints for the native structure
875  @param restraints_names Names of the restraints requested
876  @return a list with the values of all scores, and the total score
877  """
878  saved = [rb.get_reference_frame() for rb in self.components_rbs]
879  # Set the state of the model to native
880  for rb, rn in zip(self.components_rbs, self.native_rbs):
881  # remove sub-rigid bodies
882  rb_members = [m for m in rb.get_members()
883  if not core.RigidBody.get_is_setup(m.get_particle())]
884  rn_members = [m for m in rn.get_members()
885  if not core.RigidBody.get_is_setup(m.get_particle())]
886  rb_coords = [m.get_coordinates() for m in rb_members]
887  rn_coords = [m.get_coordinates() for m in rn_members]
888 
889  # align and put rb in the position of rn
890  if len(rn_coords) != len(rb_coords):
891  raise ValueError("Mismatch in the number of members. "
892  "Reference %d Aligned %d " % (len(rn_coords), len(rb_coords)))
893  T = alg.get_transformation_aligning_first_to_second(rb_coords,
894  rn_coords)
895  t = rb.get_reference_frame().get_transformation_to()
896  new_t = alg.compose(T, t)
897  rb.set_reference_frame(alg.ReferenceFrame3D(new_t))
898  scores = []
899  for name in restraints_names:
900  scores.append(self.restraints[name].evaluate(False))
901  total_score = sum(scores)
902  representation.set_reference_frames(self.components_rbs, saved)
903  return scores, total_score
904 
905  def write_pdb_for_assignment(self, assignment, fn_pdb):
906  """
907  Write the solution in the assignment
908  @param assignment Assignment class with the states for the subset
909  @param fn_pdb File to write the model
910  """
911  if not self.assignments_sampling_done:
912  raise ValueError("DominoModel: sampling not done")
913  log.info("Writting PDB %s for assignment %s", fn_pdb, assignment)
914  self.load_state(assignment)
915  atom.write_pdb(self.assembly, fn_pdb)
916 
917  def write_pdb_for_configuration(self, n, fn_pdb):
918  """
919  Write a file with a configuration for the model (configuration
920  here means a configuration in DOMINO)
921  @param n Index of the configuration desired.
922  @param fn_pdb
923  """
924  if not self.configuration_sampling_done:
925  raise ValueError("DominoModel: sampling not done")
926  log.debug("Writting PDB %s for configuration %s", fn_pdb, n)
927  configuration_set.load_configuration(n)
928  atom.write_pdb(self.assembly, fn_pdb)
929 
930  def write_pdb_for_reference_frames(self, RFs, fn_pdb):
931  """
932  Write a PDB file with a solution given by a set of reference frames
933  @param RFs Reference frames for the elements of the complex
934  @param fn_pdb File to write the solution
935  """
936  log.debug("Writting PDB %s for reference frames %s", fn_pdb, RFs)
937  representation.set_reference_frames(self.components_rbs, RFs)
938  atom.write_pdb(self.assembly, fn_pdb)
939 
940  def write_pdbs_for_reference_frames(self, RFs, fn_base):
941  """
942  Write a separate PDB for each of the elements
943  @param RFs Reference frames for the elements of the complex
944  @param fn_base base string to build the names of the PDBs files
945  """
946  log.debug("Writting PDBs with basename %s", fn_base)
947  representation.set_reference_frames(self.components_rbs, RFs)
948  for i, ch in enumerate(self.assembly.get_children()):
949  atom.write_pdb(ch, fn_base + "component-%02d.pdb" % i)
950 
951  def write_pdb_for_component(self, component_index, ref, fn_pdb):
952  """
953  Write one component of the assembly
954  @param component_index Position of the component in the assembly
955  @param ref Reference frame of the component
956  @param fn_pdb Output PDB
957  """
958  ("Writting PDB %s for component %s", fn_pdb, component_index)
959  self.components_rbs[component_index].set_reference_frame(ref)
960  hierarchy_component = self.assembly.get_child(component_index)
961  atom.write_pdb(hierarchy_component, fn_pdb)
962 
964  """Get a ScoringFunction that includes all restraints."""
965  return IMP.core.RestraintsScoringFunction(self.restraints.values())
966 
967  def write_monte_carlo_solution(self, fn_database):
968  """
969  Write the solution of a MonteCarlo run
970  @param fn_database Database file with the solution.
971  The database will contain only one record
972  """
973  total_score = 0
974  rnames = []
975  scores = []
976  for r in self.restraints.values():
977  score = r.evaluate(False)
978  rnames.append(r.get_name())
979  scores.append(score)
980  total_score += score
982  db.create(fn_database, overwrite=True)
983  db.connect(fn_database)
984  db.add_results_table(rnames, self.measure_models)
985  RFs = [rb.get_reference_frame() for rb in self.components_rbs]
986  measures = None
987  if(self.measure_models):
988  measures = measure_model(self.assembly, self.native_assembly,
989  self.components_rbs, self.native_rbs)
990  solution_id = 0
991  db.add_record(solution_id, "", RFs, total_score, scores, measures)
992  db.save_records()
993  db.close()
994 
995  def print_restraints_values(self):
996  print("Restraints: Name, weight, value, maximum_value")
997  total_score = 0
998  for r in self.restraints.values():
999  score = r.evaluate(False)
1000  # print "%20s %18f %18f %18f" % (r.get_name(), r.get_weight(),
1001  # score, r.get_maximum_score())
1002  print("%20s %18f %18f" % (r.get_name(), r.get_weight(), score))
1003  total_score += score
1004  print("total_score:", total_score)
1005 
1006 
1007 def anchor_assembly(components_rbs, anchored):
1008  """
1009  "Anchor" a set of rigid bodies, by setting the position of one of them
1010  at the origin (0,0,0).
1011  @param components_rbs Rigid bodies of the components
1012  @param anchored - List of True/False values indicating if the components
1013  of the assembly are anchored. The function sets the FIRST anchored
1014  component in the (0,0,0) coordinate and moves the other
1015  components with the same translation. If all the values are False, the
1016  function does not do anything
1017  """
1018  if True in anchored:
1019  anchored_rb = None
1020  for a, rb in zip(anchored, components_rbs):
1021  if a:
1022  anchored_rb = rb
1023  break
1024  center = anchored_rb.get_coordinates()
1025  log.info("Anchoring assembly at (0,0,0)")
1026  for rb in components_rbs:
1027  T = rb.get_reference_frame().get_transformation_to()
1028  t = T.get_translation()
1029  R = T.get_rotation()
1030  log.debug("%s: Rerefence frame BEFORE %s",
1031  rb.get_name(), rb.get_reference_frame())
1032  T2 = alg.Transformation3D(R, t - center)
1033  rb.set_reference_frame(alg.ReferenceFrame3D(T2))
1034  log.debug("%s Rerefence frame AFTER %s",
1035  rb.get_name(), rb.get_reference_frame())
1036 
1037 
1038 def measure_model(assembly, native_assembly, rbs, native_rbs):
1039  """
1040  Get the drms, cdrms and crmsd for a model
1041  @param assembly Hierarchy for an assembly
1042  @param native_assembly Hierarchy of the native structure
1043  @param rbs Rigid bodies for the components of the assembly
1044  @param native_rbs Rigid bodies for the components of the native
1045  structure
1046  @return cdrms, cdrms, crmsd
1047  - drms - typical drms measure
1048  - cdrms - drms for the centroids of the rigid bodies of the components
1049  of the assembly
1050  - crmsd - rmsd for the centroids
1051  """
1052  log.debug("Measure model")
1053  drms = comparisons.get_drms_for_backbone(assembly, native_assembly)
1054  RFs = [rb.get_reference_frame() for rb in rbs]
1055  vs = [r.get_transformation_to().get_translation() for r in RFs]
1056  nRFs = [rb.get_reference_frame() for rb in native_rbs]
1057  nvs = [r.get_transformation_to().get_translation() for r in nRFs]
1058  cdrms = atom.get_drms(vs, nvs)
1059  crmsd, RFs = alignments.align_centroids_using_pca(RFs, nRFs)
1060  log.debug("drms %s cdrms %s crmsd %s", drms, cdrms, crmsd)
1061  return [drms, cdrms, crmsd]
1062 
1063 
1064 def get_coordinates(rigid_bodies):
1065  """
1066  Return a list of the coordinates of all the members of the rigid bodies
1067  """
1068  if len(rigid_bodies) == 0:
1069  raise ValueError(
1070  "get_coordinates: There are not rigid bodies to get coordinates from")
1071  coords = []
1072  for rb in rigid_bodies:
1073  # remove sub-rigid bodies
1074  rb_members = [m for m in rb.get_members()
1075  if not core.RigidBody.get_is_setup(m.get_particle())]
1076  coords.extend([m.get_coordinates() for m in rb_members])
1077  return coords
def write_pdb_for_component
Write one component of the assembly.
def get_assignment_and_RFs
Return a line with texts for an assignment and the the reference frames of the RigidBodies in the sub...
Utility functions to handle restraints.
Definition: restraints.py:1
A harmonic upper bound on the distance between two spheres.
Lower bound harmonic function (non-zero when feature < mean)
Compute the complementarity between two molecules.
def set_rb_states
Add a set of reference frames as possible states for a rigid body.
Apply a score to a fixed number of close pairs from the two sets.
def add_exclusion_filter_table
ExclusionSubsetFilterTable ensures that two particles are not given the same state at the same time...
Restraints using electron microscopy 2D images (class averages).
def get_restraints_names_used
Get the names of the restraints used for a subset.
def write_solutions_database
Write the results of the DOMINO sampling to a SQLite database.
Various classes to hold sets of particles.
Upper bound harmonic function (non-zero when feature > mean)
Apply the score to all pairs whose spheres are within a distance threshold.
def get_assignment_scores
Returns a list with the values for the restraints on a subset of the components of the assembly...
Utility functions to handle representation.
def anchor_assembly
"Anchor" a set of rigid bodies, by setting the position of one of them at the origin (0...
Utility functions to handle alignments.
Definition: alignments.py:1
Sample best solutions using Domino.
Definition: DominoSampler.h:32
def set_close_pairs_excluded_volume_restraint
Creates an excluded volume restraint between all the components of the coarse assembly.See help for.
def set_xlink_restraint
Set a restraint on the maximum distance between 2 residues.
def get_assignment_text
Get the formatted text for an assignment.
def load_state
Load the positions of the components given by the assignment.
Filter a configuration of the subset using the Model thresholds.
Represent a subset of the particles being optimized.
Definition: Subset.h:33
def get_coordinates
Return a list of the coordinates of all the members of the rigid bodies.
Utility functions to handle IO.
Definition: io.py:1
A score on the distance between the surfaces of two spheres.
def create_merge_tree
Creates a merge tree from the restraints that are set already.
Management of a model using DOMINO.
Definition: domino_model.py:37
def get_restraint_value_for_assignment
Recover the value of a restraint.
def write_merge_tree
Writes the current merge tree to file.
def set_connectivity_restraint
Set a connectivity restraint between the elements of the assembly.
def write_pdbs_for_reference_frames
Write a separate PDB for each of the elements.
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
def write_pdb_for_reference_frames
Write a PDB file with a solution given by a set of reference frames.
Do not allow two particles to be in the same state.
def create_coarse_assembly
Simplify the assembly with a coarse representation.
def write_subset
Writes lines describing a subset (vertex of a merge tree)
def get_restraints_for_native
Get values of the restraints for the native structure.
Utility functions to store and retrieve solution information.
Definition: solutions_io.py:1
A lookup based particle refiner.
Definition: TableRefiner.h:21
def set_native_assembly_for_benchmark
Sets the native model for benchmark, by reading the native structure and set the rigid bodies...
def get_scoring_function
Get a ScoringFunction that includes all restraints.
def do_sampling
Do Domino sampling.
def write_pdb_for_configuration
Write a file with a configuration for the model (configuration here means a configuration in DOMINO) ...
def set_complementarity_restraint
Set a restraint for geometric complementarity between 2 components.
def measure_model
Get the drms, cdrms and crmsd for a model.
Fitting atomic structures into a cryo-electron microscopy density map.
Class for managing the results of the experiments.
def set_assembly
Sets an assembly from a single PDB file.
def get_assignment_info
Info related to an assignment.
def read_merge_tree
Read and create a merge tree from a file.
def setup_domino_sampler
Creates a DOMINO sampler and adds the required filter tables.
def align_rigid_bodies_states
Aligns the set of structures considered for DOMINO sampling.
Definition: domino_model.py:78
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...
Parameters used by Em2DRestraint and ProjectionFinder.
def get_assignments_with_heap
Domino sampling that recovers the assignments for the root of the merge tree, but conserving only the...
def set_pair_score_restraint
Set a pair_score restraint between the coarse representation of two components.
Utility functions for comparisons.
Definition: comparisons.py:1
Class to handle individual particles of a Model object.
Definition: Particle.h:41
def set_em2d_restraint
Set a Restraint computing the em2d score.
def add_branch_and_bound_assignments_table
BranchAndBoundAssignmentsTable enumerate states based on provided filtered using the provided the sub...
def set_not_optimized
Set a part of the model as not optimized (it does not move during the model optimization) ...
Store a set of k top scoring assignments.
def add_restraint
Adds a restraint to the model.
Definition: domino_model.py:58
Applies a PairScore to a Pair.
Definition: PairRestraint.h:29
Output IMP model data in various file formats.
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
def write_pdb_for_assignment
Write the solution in the assignment.
Divide-and-conquer inferential optimization in discrete space.
def set_assembly_components
Sets the components of an assembly, each one given as a separate PDB file.
def add_restraint_score_filter_table
Add a RestraintScoreSubsetFilterTable to DOMINO that allows to reject assignments that have score wor...
def write_monte_carlo_solution
Write the solution of a MonteCarlo run.