1 """@namespace IMP.EMageFit.domino_model
2 Classes to manage a model using DOMINO.
9 log = logging.getLogger(
"DominoModel")
40 Management of a model using DOMINO
43 def __init__(self, name="my model"):
45 self.model.set_name(name)
46 self.configuration_sampling_done =
False
47 self.assignments_sampling_done =
False
49 self.merge_tree =
None
51 self.measure_models =
False
54 self.create_coarse =
True
55 self.restraints = dict()
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
67 log.info(
"Adding restraint %s weight %s max_score %s",
68 name, weight, max_score)
70 raise ValueError(
"SQL database does not accept restraint names "
73 if max_score
is not None and max_score:
74 r.set_maximum_score(max_score)
76 self.restraints[name] = r
80 Aligns the set of structures considered for DOMINO sampling.
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
99 log.debug(
"Align the configurations read from the database before "
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()
108 "There are no particle states for %s" %
114 for rb, states
in zip(self.components_rbs, rb_states):
115 states.load_particle_state(0, rb)
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)
121 for j
in range(len(self.components_rbs)):
124 state_index = min(i, n_states[j] - 1)
125 rb_states[j].load_particle_state(
127 self.components_rbs[j])
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)
136 for i, rb
in enumerate(self.components_rbs):
140 distance, weight, stddev, max_score=
False):
142 Set a restraint on the maximum distance between 2 residues
143 @param id1 Name of the first component
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
150 @param residue2 Residue number for the aminoacid in the second
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".
168 log.info(
"Setting cross-linking restraint ")
169 log.info(
"%s", xlink.show())
170 self.xlinks_dict.add(xlink)
172 A = representation.get_component(self.assembly, xlink.first_id)
174 residue_index=xlink.first_residue)
175 p1 = s1.get_selected_particles()[0]
176 B = representation.get_component(self.assembly, xlink.second_id)
178 residue_index=xlink.second_residue)
179 p2 = s2.get_selected_particles()[0]
180 k = core.Harmonic.get_k_from_standard_deviation(stddev)
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)
192 max_sep_distance, max_penetration,
193 weight, max_score=1e10):
195 Set a restraint for geometric complementarity between 2 components
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)
205 log.info(
"Setting geometric complementarity restraint %s: %s - %s",
207 A = representation.get_component(self.assembly, name1)
208 B = representation.get_component(self.assembly, name2)
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)
219 Simplify the assembly with a coarse representation
220 @param n_residues Number of residues used for a coarse particle
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
229 def do_sampling(self, mode="assignments_heap_container", params=None):
232 @param mode The mode used to recover solutions from domino.
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"
247 if mode ==
"configuration":
248 self.configuration_set = self.sampler.create_sample()
250 log.info(
"found %s configurations. Time %s",
251 self.configuration_set.get_number_of_configurations(),
253 self.configuration_sampling_done =
True
254 elif mode ==
"assignments":
255 subset = self.rb_states_table.get_subset()
257 "Do sampling with assignments. Subset has %s elements: %s",
259 self.solution_assignments = \
260 self.sampler.get_sample_assignments(subset)
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)
270 root = self.merge_tree.get_vertices()[-1]
272 params.heap_solutions)
273 self.solution_assignments = container.get_assignments()
275 log.info(
"found %s assignments. Time %s",
276 len(self.solution_assignments), tf - t0)
277 self.assignments_sampling_done =
True
279 raise ValueError(
"Requested wrong mode for the DOMINO sampling")
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.
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])
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.
306 ordered_assignment = []
308 for rb
in self.components_rbs:
309 for i, particle
in enumerate(subset):
310 if rb.get_particle() == particle:
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]
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
330 info = text + [total_score] + scores
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
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
349 Load the positions of the components given by the assignment
350 @param assignment Assignment class
352 subset = self.rb_states_table.get_subset()
353 domino.load_particle_states(subset, assignment, self.rb_states_table)
357 Domino sampling that recovers the assignments for the root of the
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.
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")
370 subset = self.merge_tree.get_vertex_name(vertex)
371 log.info(
"computing assignments for vertex %s", subset)
374 self.restraint_cache,
"my_heap_assignments_container")
375 neighbors = self.merge_tree.get_out_neighbors(vertex)
376 if len(neighbors) == 0:
378 self.sampler.load_vertex_assignments(vertex, assignments_container)
380 if neighbors[0] > neighbors[1]:
383 neighbors = [neighbors[1], neighbors[0]]
387 self.sampler.load_vertex_assignments(vertex,
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
398 Recover the value of a restraint
399 @param name of the restraint
400 @param assignment Assignment class containing the assignment for
403 if not self.assignments_sampling_done:
404 raise ValueError(
"DominoModel: sampling not done")
405 log.debug(
"Computing restraint value for assignment %s", assignment)
407 for rb
in self.components_rbs:
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)
416 Sets the native model for benchmark, by reading the native
417 structure and set the rigid bodies.
419 self.measure_models =
True
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)
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)
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.
444 log.info(
"Set rigid body states for %s", rb.get_name())
445 RFs = [alg.ReferenceFrame3D(T)
for T
in transformations]
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())
456 Add a RestraintScoreSubsetFilterTable to DOMINO that allows to
457 reject assignments that have score worse than the thresholds for
460 log.info(
"Adding RestraintScoreSubsetFilterTable")
464 self.restraint_cache.add_restraints(list(self.restraints.values()))
466 rssft.set_name(
'myRestraintScoreSubsetFilterTable')
467 self.sampler.add_subset_filter_table(rssft)
471 Sets the components of an assembly, each one given as a separate
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.
476 if len(fn_pdbs) != len(names):
477 raise ValueError(
"Each PDB file needs a name")
479 self.assembly = representation.create_assembly(self.model, fn_pdbs,
481 self.components_rbs = representation.create_rigid_bodies(self.assembly)
482 self.not_optimized_rbs = []
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)
492 self.assembly = representation.create_assembly_from_pdb(self.model,
495 self.components_rbs = representation.create_rigid_bodies(self.assembly)
496 self.not_optimized_rbs = []
500 BranchAndBoundAssignmentsTable enumerate states based on provided
501 filtered using the provided the subset filter tables
503 log.info(
"Setting BranchAndBoundAssignmentsTable")
505 for i
in range(self.sampler.get_number_of_subset_filter_tables()):
506 ft = self.sampler.get_subset_filter_table(i)
509 self.sampler.set_assignments_table(at)
510 self.assignments_table = at
514 ExclusionSubsetFilterTable ensures that two particles are not given
515 the same state at the same time
517 log.info(
"Setting ExclusionSubsetFilterTable")
519 self.sampler.add_subset_filter_table(ft)
523 Creates a DOMINO sampler and adds the required filter tables
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")
530 self.sampler.set_restraints(list(self.restraints.values()))
531 self.sampler.set_log_level(IMP.TERSE)
532 self.sampler.set_merge_tree(self.merge_tree)
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
548 log.info(
"Reading merge tree from %s", fn)
549 ps = self.rb_states_table.get_particles()
550 log.debug(
"particles")
552 log.debug(
"%s", p.get_name())
553 rows = csv_related.read_csv(fn, delimiter=field_delim)
554 for i
in range(len(rows)):
556 rows[i] = [int(row[0]), int(row[1]), int(row[2]), row[3]]
562 names = row[3].split(unit_delim)
563 log.debug(
"row %s names %s", row, names)
567 l = [p
for p
in ps
if p.get_name() == name]
569 ValueError(
"More than one particle with same name" % names)
572 subset_names = [p.get_name()
for p
in particles]
573 log.debug(
"Merge tree Subset %s. Particles %s ", s, subset_names)
577 jt = domino.SubsetGraph()
578 jt.add_vertex(subsets[0])
579 mt = domino.get_merge_tree(jt)
581 for i
in range(1, len(subsets)):
582 mt.add_vertex(subsets[i])
589 mt.add_edge(vid, child_left)
590 if child_right != -1:
591 mt.add_edge(vid, child_right)
593 log.info(
"%s", self.merge_tree.show_graphviz())
597 Creates a merge tree from the restraints that are set already
599 rs = list(self.restraints.values())
600 ig = domino.get_interaction_graph(rs, self.rb_states_table)
604 jt = domino.get_junction_tree(ig)
606 self.merge_tree = domino.get_balanced_merge_tree(jt)
608 log.info(
"Balanced merge tree created")
609 log.info(
"%s", self.merge_tree.show_graphviz())
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
617 log.info(
"Writing merge to %s", fn)
618 if self.merge_tree
is None:
619 raise ValueError(
"No merge tree")
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]
630 Writes lines describing a subset (vertex of a merge tree)
631 @param w A csv.writer
632 @param v Vertex index
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]
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:
645 w.writerow([v, no_child_index, no_child_index, names])
647 if neighbors[0] > neighbors[1]:
648 neighbors = [neighbors[1], neighbors[0]]
649 w.writerow([v, neighbors[0], neighbors[1], names])
655 weight=1.0, max_score=
None, n_pairs=1,
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
671 c = representation.get_component(self.coarse_assembly, name)
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,
683 n_projections, weight, max_score=
False,
684 mode=
"fast", n_optimized=1):
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
702 log.info(
"Adding a em2D restraint: %s", fn_images_sel)
706 r = restraints.get_em2d_restraint(self.assembly,
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
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)
725 weight=1.0, ratio=0.05, stddev=2,
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
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
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()):
747 name =
"exc_vol_%s_%s" % (c1.get_name(), c2.get_name())
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
755 self.model,
"marker1 " + name)
757 self.model,
"marker2 " + name)
759 table_refiner.add_particle(marker1, ls1)
760 table_refiner.add_particle(marker2, ls2)
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)
778 restraint_name, distance=10,
779 weight=1.0, n_pairs=1, stddev=2, max_score=
None):
781 Set a pair_score restraint between the coarse representation
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
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
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))
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))
806 k = core.Harmonic.get_k_from_standard_deviation(stddev)
816 error_distance_allowed = 10
817 max_score = weight * n_pairs * \
818 temp_score.evaluate(distance + error_distance_allowed)
820 log.info(
"Setting pair score restraint for %s %s. k = %s, max_score "
821 "= %s, stddev %s", name1, name2, k, max_score, stddev)
827 Write the results of the DOMINO sampling to a SQLite database.
829 @param max_number Maximum number of results to write
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()
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]):
847 RFs = [rb.get_reference_frame()
for rb
in self.components_rbs]
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)
856 if self.measure_models:
857 assignment =
"native"
858 RFs = [rb.get_reference_frame()
for rb
in self.native_rbs]
861 db.add_native_record(assignment, RFs, total_score, scores)
866 """ Get the names of the restraints used for a subset
867 @param subset Subset of components of the assembly
869 return [r.get_name()
for r
870 in self.restraint_cache.get_restraints(subset, [])]
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
878 saved = [rb.get_reference_frame()
for rb
in self.components_rbs]
880 for rb, rn
in zip(self.components_rbs, self.native_rbs):
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]
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,
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))
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
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
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)
915 atom.write_pdb(self.assembly, fn_pdb)
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.
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)
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
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)
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
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)
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
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)
964 """Get a ScoringFunction that includes all restraints."""
966 list(self.restraints.values()))
970 Write the solution of a MonteCarlo run
971 @param fn_database Database file with the solution.
972 The database will contain only one record
977 for r
in self.restraints.values():
978 score = r.evaluate(
False)
979 rnames.append(r.get_name())
983 db.create(fn_database, overwrite=
True)
984 db.connect(fn_database)
985 db.add_results_table(rnames, self.measure_models)
986 RFs = [rb.get_reference_frame()
for rb
in self.components_rbs]
988 if(self.measure_models):
989 measures =
measure_model(self.assembly, self.native_assembly,
990 self.components_rbs, self.native_rbs)
992 db.add_record(solution_id,
"", RFs, total_score, scores, measures)
996 def print_restraints_values(self):
997 print(
"Restraints: Name, weight, value, maximum_value")
999 for r
in self.restraints.values():
1000 score = r.evaluate(
False)
1003 print(
"%20s %18f %18f" % (r.get_name(), r.get_weight(), score))
1004 total_score += score
1005 print(
"total_score:", total_score)
1010 "Anchor" a set of rigid bodies, by setting the position of one of them
1011 at the origin (0,0,0).
1012 @param components_rbs Rigid bodies of the components
1013 @param anchored - List of True/False values indicating if the components
1014 of the assembly are anchored. The function sets the FIRST anchored
1015 component in the (0,0,0) coordinate and moves the other
1016 components with the same translation. If all the values are False, the
1017 function does not do anything
1019 if True in anchored:
1021 for a, rb
in zip(anchored, components_rbs):
1025 center = anchored_rb.get_coordinates()
1026 log.info(
"Anchoring assembly at (0,0,0)")
1027 for rb
in components_rbs:
1028 T = rb.get_reference_frame().get_transformation_to()
1029 t = T.get_translation()
1030 R = T.get_rotation()
1031 log.debug(
"%s: Rerefence frame BEFORE %s",
1032 rb.get_name(), rb.get_reference_frame())
1033 T2 = alg.Transformation3D(R, t - center)
1034 rb.set_reference_frame(alg.ReferenceFrame3D(T2))
1035 log.debug(
"%s Rerefence frame AFTER %s",
1036 rb.get_name(), rb.get_reference_frame())
1041 Get the drms, cdrms and crmsd for a model
1042 @param assembly Hierarchy for an assembly
1043 @param native_assembly Hierarchy of the native structure
1044 @param rbs Rigid bodies for the components of the assembly
1045 @param native_rbs Rigid bodies for the components of the native
1047 @return cdrms, cdrms, crmsd
1048 - drms - typical drms measure
1049 - cdrms - drms for the centroids of the rigid bodies of the components
1051 - crmsd - rmsd for the centroids
1053 log.debug(
"Measure model")
1054 drms = comparisons.get_drms_for_backbone(assembly, native_assembly)
1055 RFs = [rb.get_reference_frame()
for rb
in rbs]
1056 vs = [r.get_transformation_to().get_translation()
for r
in RFs]
1057 nRFs = [rb.get_reference_frame()
for rb
in native_rbs]
1058 nvs = [r.get_transformation_to().get_translation()
for r
in nRFs]
1059 cdrms = atom.get_drms(vs, nvs)
1060 crmsd, RFs = alignments.align_centroids_using_pca(RFs, nRFs)
1061 log.debug(
"drms %s cdrms %s crmsd %s", drms, cdrms, crmsd)
1062 return [drms, cdrms, crmsd]
1067 Return a list of the coordinates of all the members of the rigid bodies
1069 if len(rigid_bodies) == 0:
1071 "get_coordinates: There are not rigid bodies to get coordinates from")
1073 for rb
in rigid_bodies:
1075 rb_members = [m
for m
in rb.get_members()
1076 if not core.RigidBody.get_is_setup(m.get_particle())]
1077 coords.extend([m.get_coordinates()
for m
in rb_members])
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.
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.
Sample best solutions using Domino.
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.
def get_coordinates
Return a list of the coordinates of all the members of the rigid bodies.
Utility functions to handle IO.
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.
Class defining a cross-link.
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.
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.
A lookup based particle refiner.
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.
Description of crosslinking restraints as a python dictionary.
def align_rigid_bodies_states
Aligns the set of structures considered for DOMINO sampling.
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.
Class to handle individual particles of a Model object.
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.
Applies a PairScore to a Pair.
Output IMP model data in various file formats.
Functionality for loading, creating, manipulating and scoring atomic structures.
Utility functions to handle cross links.
Select hierarchy particles identified by the biological name.
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.