1 """@namespace IMP.EMageFit.domino_model
2 Classes to manage a model using DOMINO.
27 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
68 log.info(
"Adding restraint %s weight %s max_score %s",
69 name, weight, max_score)
71 raise ValueError(
"SQL database does not accept restraint names "
74 if max_score
is not None and max_score:
75 r.set_maximum_score(max_score)
77 self.restraints[name] = r
81 Aligns the set of structures considered for DOMINO sampling.
83 1) reads the possible reference frames that the
84 rb_states_table stores for each rigid body. This table
85 must be filled before using this function. Usually it is
86 filled with the results from a previous Monte Carlo
87 sampling. If the solutions from Monte Carlo seem to have
88 the same structure but they are not aligned to each other,
89 this function can help setting better starting positions
91 2) Gets the first state for each of the rigid bodies and sets
92 a reference structure using such states.
93 3) Aligns all the rest of the structures respect to the
95 4) Replaces the values of the reference frames stored in the
96 rb_states_table with the new values obtained from the
97 alignments. It does it for all states of a rigid body.
98 @note: If this function is applied, the parameters "anchor"
99 and "fixed" are ignored, as they are superseded by the use
100 of the alignments calculated here.
102 log.debug(
"Align the configurations read from the database before "
106 for rb
in self.components_rbs:
107 ps = self.rb_states_table.get_particle_states(rb)
108 n = ps.get_number_of_particle_states()
111 "There are no particle states for %s" %
117 for rb, states
in zip(self.components_rbs, rb_states):
118 states.load_particle_state(0, rb)
120 aligned_transformations = [[]
for rb
in self.components_rbs]
121 for i
in range(1, max(n_states)):
122 log.debug(
"Aligning rigid bodies configuration %s" % i)
124 for j
in range(len(self.components_rbs)):
127 state_index = min(i, n_states[j] - 1)
128 rb_states[j].load_particle_state(
130 self.components_rbs[j])
132 T = alg.get_transformation_aligning_first_to_second(
133 coords, reference_coords)
134 for j, rb
in enumerate(self.components_rbs):
135 t = rb.get_reference_frame().get_transformation_to()
136 new_t = alg.compose(T, t)
137 aligned_transformations[j].append(new_t)
139 for i, rb
in enumerate(self.components_rbs):
143 distance, weight, stddev, max_score=
False):
145 Set a restraint on the maximum distance between 2 residues
146 @param id1 Name of the first component
148 @param residue1 Residue number for the aminoacid in the first
149 component.The number is the number in the PDB file, not the
150 number relative to the beginning of the chain
151 @param id2 Name of the second component
153 @param residue2 Residue number for the aminoacid in the second
155 @param distance Maximum distance tolerated
156 @param weight Weight of the restraint
157 @param stddev Standard deviation used to approximate the
158 HarmonicUpperBound function to a Gaussian
159 @param max_score See help for add_restraint(). If none is given,
160 the maximum score is set to allow a maximum distance of 10 Angstrom
161 greater than the parameter "distance".
171 log.info(
"Setting cross-linking restraint ")
172 log.info(
"%s", xlink.show())
173 self.xlinks_dict.add(xlink)
175 A = representation.get_component(self.assembly, xlink.first_id)
177 residue_index=xlink.first_residue)
178 p1 = s1.get_selected_particles()[0]
179 B = representation.get_component(self.assembly, xlink.second_id)
181 residue_index=xlink.second_residue)
182 p2 = s2.get_selected_particles()[0]
183 k = core.Harmonic.get_k_from_standard_deviation(stddev)
188 error_distance_allowed = 100
189 max_score = weight * \
190 score.evaluate(distance + error_distance_allowed)
191 log.info(
"%s, max_score %s", xlink.show(), max_score)
195 max_sep_distance, max_penetration,
196 weight, max_score=1e10):
198 Set a restraint for geometric complementarity between 2 components
200 @param name2 - The restraint is applied to this components
201 @param rname - The name given to the restraint
202 @param max_sep_distance - maximum distance between molecules
203 tolerated by the restraint
204 @param max_penetration - Maximum penetration allowed (angstrom)
208 log.info(
"Setting geometric complementarity restraint %s: %s - %s",
210 A = representation.get_component(self.assembly, name1)
211 B = representation.get_component(self.assembly, name2)
213 atom.get_leaves(A), atom.get_leaves(B), rname)
214 restraint.set_maximum_separation(max_sep_distance)
215 restraint.set_maximum_penetration(max_penetration)
216 log.debug(
"Maximum separation: %s Maximum penetration score: %s",
217 max_sep_distance, max_penetration)
222 Simplify the assembly with a coarse representation
223 @param n_residues Number of residues used for a coarse particle
225 if self.create_coarse:
226 log.info(
"Creating simplified assembly")
227 self.coarse_assembly = representation.create_simplified_assembly(
228 self.assembly, self.components_rbs, n_residues)
229 self.create_coarse =
False
231 def do_sampling(self, mode="assignments_heap_container", params=None):
234 @param mode The mode used to recover solutions from domino.
236 the value "configuration", "assignments", or
237 "assignments_heap_container". The mode
238 "configuration" recovers all possible configurations
239 (takes a while to generate them).
240 The mode "assignments" only provides the assignments
241 (state numbers) of the solutions. It is faster but requires
242 generating the solutions afterwards
243 The mode "assignments_heap_container" selects the best solutions
244 after each merging in DOMINO, discarding the rest.
245 In practice I used the mode "assignments_heap_container"
249 if mode ==
"configuration":
250 self.configuration_set = self.sampler.create_sample()
252 log.info(
"found %s configurations. Time %s",
253 self.configuration_set.get_number_of_configurations(),
255 self.configuration_sampling_done =
True
256 elif mode ==
"assignments":
257 subset = self.rb_states_table.get_subset()
259 "Do sampling with assignments. Subset has %s elements: %s",
261 self.solution_assignments = \
262 self.sampler.get_sample_assignments(subset)
264 log.info(
"found %s assignments. Time %s",
265 len(self.solution_assignments), tf - t0)
266 self.assignments_sampling_done =
True
267 elif mode ==
"assignments_heap_container":
268 subset = self.rb_states_table.get_subset()
269 log.info(
"DOMINO sampling an assignments_heap_container. "
270 "Subset has %s elements: %s", len(subset), subset)
272 root = self.merge_tree.get_vertices()[-1]
274 params.heap_solutions)
275 self.solution_assignments = container.get_assignments()
277 log.info(
"found %s assignments. Time %s",
278 len(self.solution_assignments), tf - t0)
279 self.assignments_sampling_done =
True
281 raise ValueError(
"Requested wrong mode for the DOMINO sampling")
285 Get the formatted text for an assignment
286 @param subset Subset of components of the assembly
287 @param assignment Assignment class with the states for the subset
288 @note: The order in assignment and subset is not always the order
289 of the rigid bodies in self.components_rbs. The function reorders
290 the terms in the assignment so there is correspondence.
292 ordered_assignment = []
293 for rb
in self.components_rbs:
294 for i, particle
in enumerate(subset):
295 if rb.get_particle() == particle:
296 ordered_assignment.append(assignment[i])
297 text =
"|".join([str(i)
for i
in ordered_assignment])
302 Return a line with texts for an assignment and the
303 the reference frames of the RigidBodies in the subset.
304 @param subset Subset of components of the assembly
305 @param assignment Assignment class with the states for the subset
306 @note: see the get_assignment_assignment_text() note.
308 ordered_assignment = []
310 for rb
in self.components_rbs:
311 for i, particle
in enumerate(subset):
312 if rb.get_particle() == particle:
314 ordered_assignment.append(j)
315 pstates = self.rb_states_table.get_particle_states(rb)
316 pstates.load_particle_state(j, rb.get_particle())
317 RFs.append(rb.get_reference_frame())
318 RFs_text = unit_delim.join(
319 [io.ReferenceFrameToText(ref).get_text()
for ref
in RFs])
320 assignment_numbers =
"|".join([str(i)
for i
in ordered_assignment])
321 return [assignment_numbers, RFs_text]
325 Info related to an assignment. Returns a list with text for the
326 assignment and all the scores for the restraints
327 @param subset Subset of components of the assembly
328 @param assignment Assignment class with the states for the subset
332 info = text + [total_score] + scores
337 Returns a list with the values for the restraints on a subset of
338 the components of the assembly
339 @param subset Subset of components of the assembly
340 @param assignment Assignment class with the states for the subset
342 restraints_to_evaluate = \
343 self.restraint_cache.get_restraints(subset, [])
344 scores = [self.restraint_cache.get_score(r, subset, assignment)
345 for r
in restraints_to_evaluate]
346 total_score = sum(scores)
347 return scores, total_score
351 Load the positions of the components given by the assignment
352 @param assignment Assignment class
354 subset = self.rb_states_table.get_subset()
355 domino.load_particle_states(subset, assignment, self.rb_states_table)
359 Domino sampling that recovers the assignments for the root of the
361 conserving only the best k scores for each vertex of the tree.
362 @param[in] vertex Vertex with the root of the current
363 merge tree. This function is recursive.
366 if self.sampler.get_number_of_subset_filter_tables() == 0:
367 raise ValueError(
"No subset filter tables")
368 if self.merge_tree
is None:
369 raise ValueError(
"No merge tree")
372 subset = self.merge_tree.get_vertex_name(vertex)
373 log.info(
"computing assignments for vertex %s", subset)
376 subset, k, self.restraint_cache,
"my_heap_assignments_container")
377 neighbors = self.merge_tree.get_out_neighbors(vertex)
378 if len(neighbors) == 0:
380 self.sampler.load_vertex_assignments(vertex, assignments_container)
382 if neighbors[0] > neighbors[1]:
385 neighbors = [neighbors[1], neighbors[0]]
389 self.sampler.load_vertex_assignments(vertex,
392 assignments_container)
393 tf = time.time() - t0
394 log.info(
"Merge tree vertex: %s assignments: %s Time %s sec.", subset,
395 assignments_container.get_number_of_assignments(), tf)
396 return assignments_container
400 Recover the value of a restraint
401 @param name of the restraint
402 @param assignment Assignment class containing the assignment for
405 if not self.assignments_sampling_done:
406 raise ValueError(
"DominoModel: sampling not done")
407 log.debug(
"Computing restraint value for assignment %s", assignment)
409 for rb
in self.components_rbs:
411 rb.get_reference_frame().get_transformation_to())
412 val = self.restraints[name].evaluate(
False)
413 log.debug(
"restraint %s = %s", name, val)
418 Sets the native model for benchmark, by reading the native
419 structure and set the rigid bodies.
421 self.measure_models =
True
423 if hasattr(params.benchmark,
"fn_pdb_native"):
424 self.native_assembly = representation.create_assembly_from_pdb(
425 self.native_model, params.benchmark.fn_pdb_native,
427 elif hasattr(params.benchmark,
"fn_pdbs_native"):
428 self.native_assembly = representation.create_assembly(
429 self.native_model, params.benchmark.fn_pdbs_native,
433 "set_native_assembly_for_benchmark: Requested "
434 "to set a native assembly for benchmark but did not provide "
435 "its pdb, or the pdbs of the components")
436 self.native_rbs = representation.create_rigid_bodies(
437 self.native_assembly)
442 Add a set of reference frames as possible states for a rigid body
443 @param rb The rigid body
444 @param transformations Transformations are used to build the
445 reference frames for the rigid body.
447 log.info(
"Set rigid body states for %s", rb.get_name())
448 RFs = [alg.ReferenceFrame3D(T)
for T
in transformations]
452 self.rb_states_table.set_particle_states(rb, rb_states)
453 st = self.rb_states_table.get_particle_states(rb)
454 log.info(
"RigidBodyStates created %s",
455 st.get_number_of_particle_states())
459 Add a RestraintScoreSubsetFilterTable to DOMINO that allows to
460 reject assignments that have score worse than the thresholds for
463 log.info(
"Adding RestraintScoreSubsetFilterTable")
467 self.restraint_cache.add_restraints(list(self.restraints.values()))
469 rssft.set_name(
'myRestraintScoreSubsetFilterTable')
470 self.sampler.add_subset_filter_table(rssft)
474 Sets the components of an assembly, each one given as a separate
476 @param fn_pdbs List with the names of the pdb files
477 @param names List with the names of the components of the assembly.
479 if len(fn_pdbs) != len(names):
480 raise ValueError(
"Each PDB file needs a name")
482 self.assembly = representation.create_assembly(self.model, fn_pdbs,
484 self.components_rbs = representation.create_rigid_bodies(self.assembly)
485 self.not_optimized_rbs = []
489 Sets an assembly from a single PDB file. The function assumes that
490 the components of the assembly are the chains of the PDB file.
491 @param fn_pdb PDB with the file for the assembly
492 @param names Names for each of the components (chains)
495 self.assembly = representation.create_assembly_from_pdb(self.model,
498 self.components_rbs = representation.create_rigid_bodies(self.assembly)
499 self.not_optimized_rbs = []
503 BranchAndBoundAssignmentsTable enumerate states based on provided
504 filtered using the provided the subset filter tables
506 log.info(
"Setting BranchAndBoundAssignmentsTable")
508 for i
in range(self.sampler.get_number_of_subset_filter_tables()):
509 ft = self.sampler.get_subset_filter_table(i)
512 self.sampler.set_assignments_table(at)
513 self.assignments_table = at
517 ExclusionSubsetFilterTable ensures that two particles are not given
518 the same state at the same time
520 log.info(
"Setting ExclusionSubsetFilterTable")
522 self.sampler.add_subset_filter_table(ft)
526 Creates a DOMINO sampler and adds the required filter tables
528 if self.merge_tree
is None:
529 raise ValueError(
"Merge tree for DOMINO not set. It is required "
530 "to setup the sampler")
531 log.info(
"Domino sampler")
533 self.sampler.set_restraints(list(self.restraints.values()))
534 self.sampler.set_log_level(IMP.TERSE)
535 self.sampler.set_merge_tree(self.merge_tree)
542 Read and create a merge tree from a file.
543 The format of the file is the format written by write merge_tree()
544 It does not matter the order of the indices in the file, as
545 they are sorted by this function.
546 The only condition is that the index for the vertex that is the
547 root of the tree MUST be the largest. This is guaranteed when
548 creating a merge tree with create_merge_tree()
549 @param fn File containing the merge tree
551 log.info(
"Reading merge tree from %s", fn)
552 ps = self.rb_states_table.get_particles()
553 log.debug(
"particles")
555 log.debug(
"%s", p.get_name())
556 rows = csv_related.read_csv(fn, delimiter=field_delim)
557 for i
in range(len(rows)):
559 rows[i] = [int(row[0]), int(row[1]), int(row[2]), row[3]]
565 names = row[3].split(unit_delim)
566 log.debug(
"row %s names %s", row, names)
570 lp = [p
for p
in ps
if p.get_name() == name]
572 ValueError(
"More than one particle with same name" % names)
575 subset_names = [p.get_name()
for p
in particles]
576 log.debug(
"Merge tree Subset %s. Particles %s ", s, subset_names)
580 jt = domino.SubsetGraph()
581 jt.add_vertex(subsets[0])
582 mt = domino.get_merge_tree(jt)
584 for i
in range(1, len(subsets)):
585 mt.add_vertex(subsets[i])
592 mt.add_edge(vid, child_left)
593 if child_right != -1:
594 mt.add_edge(vid, child_right)
596 log.info(
"%s", self.merge_tree.show_graphviz())
600 Creates a merge tree from the restraints that are set already
602 rs = list(self.restraints.values())
603 ig = domino.get_interaction_graph(rs, self.rb_states_table)
607 jt = domino.get_junction_tree(ig)
609 self.merge_tree = domino.get_balanced_merge_tree(jt)
611 log.info(
"Balanced merge tree created")
612 log.info(
"%s", self.merge_tree.show_graphviz())
616 Writes the current merge tree to file. The format is:
617 parent | child_left | child_right | labels
618 @param fn File for storing the merge tree
620 log.info(
"Writing merge to %s", fn)
621 if self.merge_tree
is None:
622 raise ValueError(
"No merge tree")
625 "# Vertex index | Child0 | Child1 | label "
626 "(names of the particles)\n")
627 w = csv.writer(f, delimiter=field_delim)
628 root = self.merge_tree.get_vertices()[-1]
634 Writes lines describing a subset (vertex of a merge tree)
635 @param w A csv.writer
636 @param v Vertex index
639 subset = self.merge_tree.get_vertex_name(v)
640 names = [p.get_name()
for p
in subset]
641 reps = [x
for x
in names
if names.count(x) > 1]
644 "The names of the particles in the subset %s are not unique"
646 names = unit_delim.join(names)
647 neighbors = self.merge_tree.get_out_neighbors(v)
648 if len(neighbors) == 0:
650 w.writerow([v, no_child_index, no_child_index, names])
652 if neighbors[0] > neighbors[1]:
653 neighbors = [neighbors[1], neighbors[0]]
654 w.writerow([v, neighbors[0], neighbors[1], names])
660 weight=1.0, max_score=
None, n_pairs=1,
663 Set a connectivity restraint between the elements of the assembly
664 @param names List with all the elements to be connected
665 @param distance Maximum distance tolerated by the restraint
666 @param rname Name for the restraint
667 @param weight Weight for the restraint
668 @param max_score Maximum score for the restraint
669 @param n_pairs Number of pairs of particles used in the
670 KClosePairScore of the connectivity restraint
671 @param stddev Standard deviation used to approximate the
672 HarmonicUpperBound function to a Gaussian
676 c = representation.get_component(self.coarse_assembly, name)
678 log.info(
"Connectivity restraint for %s", components)
679 spring_constant = core.Harmonic.get_k_from_standard_deviation(stddev)
680 if max_score
is None:
681 max_score = weight * spring_constant * n_pairs * (stddev) ** 2
682 r = restraints.get_connectivity_restraint(
683 components, distance, n_pairs,
688 n_projections, weight, max_score=
False,
689 mode=
"fast", n_optimized=1):
691 Set a Restraint computing the em2d score.
692 @param name Name for the restraint
693 @param fn_images_sel Selection file with the names of the images
694 used for the restraint
695 @param pixel_size - pixel size in the images
696 @param resolution - resolution used to downsample the projections
697 of the model when projecting
698 @param weight Weight of the restraint
699 @param max_score - Maximum value tolerated for the restraint
700 @param n_projections - Number of projections to generate
701 when projecting the model to do the coarse alignments
702 @param mode - Mode used for computing the restraints.
703 @param n_optimized - number of results from the coarse
704 alignment that are refined with the Simplex algorithm
705 to get a more accurate value for the restraint
707 log.info(
"Adding a em2D restraint: %s", fn_images_sel)
711 r = restraints.get_em2d_restraint(self.assembly,
719 Set a part of the model as not optimized (it does not move during
720 the model optimization)
721 @param name of the component to optimized
723 if name
not in self.names:
724 raise ValueError(
"DominoModel: There is not component "
725 "in the assembly with this name")
726 rb_name = representation.get_rb_name(name)
727 self.not_optimized_rbs.append(rb_name)
730 self, distance=10, weight=1.0, ratio=0.05, stddev=2,
733 Creates an excluded volume restraint between all the components
734 of the coarse assembly.See help for
735 @param distance Maximum distance tolerated between particles
736 @param weight Weight for the restraint
738 @param max_score Maximum value for the restraint. If the parameter
739 is None, an automatic value is computed (using the ratio).
740 @param ratio Fraction of the number of possible pairs of
741 particles that are tolerated to overlap. The maximum
742 score is modified according to this ratio.
743 I have found that ratios of 0.05-0.1 work well allowing for
744 some interpenetration
745 @param stddev Standard deviation used to approximate the
746 HarmonicUpperBound function to a Gaussian
748 k = core.Harmonic.get_k_from_standard_deviation(stddev)
749 for i, c1
in enumerate(self.coarse_assembly.get_children()):
750 for j, c2
in enumerate(self.coarse_assembly.get_children()):
752 name =
"exc_vol_%s_%s" % (c1.get_name(), c2.get_name())
754 ls1 = atom.get_leaves(c1)
755 ls2 = atom.get_leaves(c2)
756 possible_pairs = len(ls1) * len(ls2)
757 n_pairs = possible_pairs * ratio
760 self.model,
"marker1 " + name)
762 self.model,
"marker2 " + name)
764 table_refiner.add_particle(marker1, ls1)
765 table_refiner.add_particle(marker2, ls2)
775 minimum_distance_allowed = 0
776 max_score = weight * n_pairs * \
777 score.evaluate(minimum_distance_allowed)
779 "Setting close_pairs_excluded_volume_restraint(), "
780 "max_score %s", max_score)
784 restraint_name, distance=10,
785 weight=1.0, n_pairs=1, stddev=2,
788 Set a pair_score restraint between the coarse representation
790 @param name1 Name of the first component
791 @param name2 Name of the second component
792 @param restraint_name Name for the restraint
793 @param distance Maximum distance tolerated between particles
794 @param weight Weight of the restraint
796 @param max_score Maximum value tolerated for the restraint
797 @param stddev Standard deviation used to approximate the
798 HarmonicUpperBound function to a Gaussian
805 A = representation.get_component(self.coarse_assembly, name1)
806 marker1 =
IMP.Particle(self.model,
"marker1 " + restraint_name)
807 table_refiner.add_particle(marker1, atom.get_leaves(A))
809 B = representation.get_component(self.coarse_assembly, name2)
810 marker2 =
IMP.Particle(self.model,
"marker2 " + restraint_name)
811 table_refiner.add_particle(marker2, atom.get_leaves(B))
813 k = core.Harmonic.get_k_from_standard_deviation(stddev)
823 error_distance_allowed = 10
824 max_score = weight * n_pairs * \
825 temp_score.evaluate(distance + error_distance_allowed)
827 log.info(
"Setting pair score restraint for %s %s. k = %s, max_score "
828 "= %s, stddev %s", name1, name2, k, max_score, stddev)
834 Write the results of the DOMINO sampling to a SQLite database.
836 @param max_number Maximum number of results to write
838 log.info(
"Creating the database of solutions")
839 if self.measure_models
and not hasattr(self,
"native_model"):
840 raise ValueError(
"The native model has not been set")
842 db.create(fn_database, overwrite=
True)
843 db.connect(fn_database)
844 subset = self.rb_states_table.get_subset()
846 db.add_results_table(restraints_names, self.measure_models)
847 n = len(self.solution_assignments)
848 if max_number
is not None:
849 n = min(n, max_number)
850 for sol_id, assignment
in enumerate(self.solution_assignments[0:n]):
854 RFs = [rb.get_reference_frame()
for rb
in self.components_rbs]
858 if self.measure_models:
859 measures =
measure_model(self.assembly, self.native_assembly,
860 self.components_rbs, self.native_rbs)
861 db.add_record(sol_id, txt, RFs, total_score, scores, measures)
863 if self.measure_models:
864 assignment =
"native"
865 RFs = [rb.get_reference_frame()
for rb
in self.native_rbs]
868 db.add_native_record(assignment, RFs, total_score, scores)
873 """ Get the names of the restraints used for a subset
874 @param subset Subset of components of the assembly
876 return [r.get_name()
for r
877 in self.restraint_cache.get_restraints(subset, [])]
881 Get values of the restraints for the native structure
882 @param restraints_names Names of the restraints requested
883 @return a list with the values of all scores, and the total score
885 saved = [rb.get_reference_frame()
for rb
in self.components_rbs]
887 for rb, rn
in zip(self.components_rbs, self.native_rbs):
889 rb_members = [m
for m
in rb.get_members()
890 if not core.RigidBody.get_is_setup(m.get_particle())]
891 rn_members = [m
for m
in rn.get_members()
892 if not core.RigidBody.get_is_setup(m.get_particle())]
893 rb_coords = [m.get_coordinates()
for m
in rb_members]
894 rn_coords = [m.get_coordinates()
for m
in rn_members]
897 if len(rn_coords) != len(rb_coords):
898 raise ValueError(
"Mismatch in the number of members. "
899 "Reference %d Aligned %d "
900 % (len(rn_coords), len(rb_coords)))
901 T = alg.get_transformation_aligning_first_to_second(rb_coords,
903 t = rb.get_reference_frame().get_transformation_to()
904 new_t = alg.compose(T, t)
905 rb.set_reference_frame(alg.ReferenceFrame3D(new_t))
907 for name
in restraints_names:
908 scores.append(self.restraints[name].evaluate(
False))
909 total_score = sum(scores)
910 representation.set_reference_frames(self.components_rbs, saved)
911 return scores, total_score
915 Write the solution in the assignment
916 @param assignment Assignment class with the states for the subset
917 @param fn_pdb File to write the model
919 if not self.assignments_sampling_done:
920 raise ValueError(
"DominoModel: sampling not done")
921 log.info(
"Writing PDB %s for assignment %s", fn_pdb, assignment)
923 atom.write_pdb(self.assembly, fn_pdb)
927 Write a file with a configuration for the model (configuration
928 here means a configuration in DOMINO)
929 @param n Index of the configuration desired.
932 if not self.configuration_sampling_done:
933 raise ValueError(
"DominoModel: sampling not done")
934 log.debug(
"Writing PDB %s for configuration %s", fn_pdb, n)
936 atom.write_pdb(self.assembly, fn_pdb)
940 Write a PDB file with a solution given by a set of reference frames
941 @param RFs Reference frames for the elements of the complex
942 @param fn_pdb File to write the solution
944 log.debug(
"Writing PDB %s for reference frames %s", fn_pdb, RFs)
945 representation.set_reference_frames(self.components_rbs, RFs)
946 atom.write_pdb(self.assembly, fn_pdb)
950 Write a separate PDB for each of the elements
951 @param RFs Reference frames for the elements of the complex
952 @param fn_base base string to build the names of the PDBs files
954 log.debug(
"Writing PDBs with basename %s", fn_base)
955 representation.set_reference_frames(self.components_rbs, RFs)
956 for i, ch
in enumerate(self.assembly.get_children()):
957 atom.write_pdb(ch, fn_base +
"component-%02d.pdb" % i)
961 Write one component of the assembly
962 @param component_index Position of the component in the assembly
963 @param ref Reference frame of the component
964 @param fn_pdb Output PDB
966 (
"Writing PDB %s for component %s", fn_pdb, component_index)
967 self.components_rbs[component_index].set_reference_frame(ref)
968 hierarchy_component = self.assembly.get_child(component_index)
969 atom.write_pdb(hierarchy_component, fn_pdb)
972 """Get a ScoringFunction that includes all restraints."""
974 list(self.restraints.values()))
978 Write the solution of a MonteCarlo run
979 @param fn_database Database file with the solution.
980 The database will contain only one record
985 for r
in self.restraints.values():
986 score = r.evaluate(
False)
987 rnames.append(r.get_name())
991 db.create(fn_database, overwrite=
True)
992 db.connect(fn_database)
993 db.add_results_table(rnames, self.measure_models)
994 RFs = [rb.get_reference_frame()
for rb
in self.components_rbs]
996 if self.measure_models:
997 measures =
measure_model(self.assembly, self.native_assembly,
998 self.components_rbs, self.native_rbs)
1000 db.add_record(solution_id,
"", RFs, total_score, scores, measures)
1004 def print_restraints_values(self):
1005 print(
"Restraints: Name, weight, value, maximum_value")
1007 for r
in self.restraints.values():
1008 score = r.evaluate(
False)
1011 print(
"%20s %18f %18f" % (r.get_name(), r.get_weight(), score))
1012 total_score += score
1013 print(
"total_score:", total_score)
1018 "Anchor" a set of rigid bodies, by setting the position of one of them
1019 at the origin (0,0,0).
1020 @param components_rbs Rigid bodies of the components
1021 @param anchored - List of True/False values indicating if the
1022 components of the assembly are anchored. The function sets the
1023 FIRST anchored component in the (0,0,0) coordinate and moves
1024 the other components with the same translation. If all the
1025 values are False, the function does not do anything
1027 if True in anchored:
1029 for a, rb
in zip(anchored, components_rbs):
1033 center = anchored_rb.get_coordinates()
1034 log.info(
"Anchoring assembly at (0,0,0)")
1035 for rb
in components_rbs:
1036 T = rb.get_reference_frame().get_transformation_to()
1037 t = T.get_translation()
1038 R = T.get_rotation()
1039 log.debug(
"%s: Rerefence frame BEFORE %s",
1040 rb.get_name(), rb.get_reference_frame())
1041 T2 = alg.Transformation3D(R, t - center)
1042 rb.set_reference_frame(alg.ReferenceFrame3D(T2))
1043 log.debug(
"%s Rerefence frame AFTER %s",
1044 rb.get_name(), rb.get_reference_frame())
1049 Get the drms, cdrms and crmsd for a model
1050 @param assembly Hierarchy for an assembly
1051 @param native_assembly Hierarchy of the native structure
1052 @param rbs Rigid bodies for the components of the assembly
1053 @param native_rbs Rigid bodies for the components of the native
1055 @return cdrms, cdrms, crmsd
1056 - drms - typical drms measure
1057 - cdrms - drms for the centroids of the rigid bodies of the components
1059 - crmsd - rmsd for the centroids
1061 log.debug(
"Measure model")
1062 drms = comparisons.get_drms_for_backbone(assembly, native_assembly)
1063 RFs = [rb.get_reference_frame()
for rb
in rbs]
1064 vs = [r.get_transformation_to().get_translation()
for r
in RFs]
1065 nRFs = [rb.get_reference_frame()
for rb
in native_rbs]
1066 nvs = [r.get_transformation_to().get_translation()
for r
in nRFs]
1067 cdrms = atom.get_drms(vs, nvs)
1068 crmsd, RFs = alignments.align_centroids_using_pca(RFs, nRFs)
1069 log.debug(
"drms %s cdrms %s crmsd %s", drms, cdrms, crmsd)
1070 return [drms, cdrms, crmsd]
1075 Return a list of the coordinates of all the members
1078 if len(rigid_bodies) == 0:
1080 "get_coordinates: There are no rigid bodies to get "
1083 for rb
in rigid_bodies:
1085 rb_members = [m
for m
in rb.get_members()
1086 if not core.RigidBody.get_is_setup(m.get_particle())]
1087 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.
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.
Create a scoring function on a list of restraints.
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.
Score a pair of particles based on the distance between their centers.
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.
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.