1 """@namespace IMP.pmi.macros
2 Protocols for sampling structures and analyzing them.
5 from __future__
import print_function, division
19 from operator
import itemgetter
20 from collections
import defaultdict
28 class _RMFRestraints(object):
29 """All restraints that are written out to the RMF file"""
30 def __init__(self, model, user_restraints):
32 self._user_restraints = user_restraints
if user_restraints
else []
35 return (len(self._user_restraints)
36 + self._rmf_rs.get_number_of_restraints())
40 __nonzero__ = __bool__
42 def __getitem__(self, i):
43 class FakePMIWrapper(object):
44 def __init__(self, r):
45 self.r = IMP.RestraintSet.get_from(r)
46 get_restraint =
lambda self: self.r
47 lenuser = len(self._user_restraints)
49 return self._user_restraints[i]
50 elif 0 <= i - lenuser < self._rmf_rs.get_number_of_restraints():
51 r = self._rmf_rs.get_restraint(i - lenuser)
52 return FakePMIWrapper(r)
54 raise IndexError(
"Out of range")
58 """A macro to help setup and run replica exchange.
59 Supports Monte Carlo and molecular dynamics.
60 Produces trajectory RMF files, best PDB structures,
61 and output stat files.
67 monte_carlo_sample_objects=
None,
68 molecular_dynamics_sample_objects=
None,
70 rmf_output_objects=
None,
71 crosslink_restraints=
None,
72 monte_carlo_temperature=1.0,
73 simulated_annealing=
False,
74 simulated_annealing_minimum_temperature=1.0,
75 simulated_annealing_maximum_temperature=2.5,
76 simulated_annealing_minimum_temperature_nframes=100,
77 simulated_annealing_maximum_temperature_nframes=100,
78 replica_exchange_minimum_temperature=1.0,
79 replica_exchange_maximum_temperature=2.5,
80 replica_exchange_swap=
True,
82 number_of_best_scoring_models=500,
85 molecular_dynamics_steps=10,
86 molecular_dynamics_max_time_step=1.0,
87 number_of_frames=1000,
88 save_coordinates_mode=
"lowest_temperature",
89 nframes_write_coordinates=1,
90 write_initial_rmf=
True,
91 initial_rmf_name_suffix=
"initial",
92 stat_file_name_suffix=
"stat",
93 best_pdb_name_suffix=
"model",
95 do_create_directories=
True,
96 global_output_directory=
"./",
99 replica_stat_file_suffix=
"stat_replica",
100 em_object_for_rmf=
None,
102 replica_exchange_object=
None,
105 @param model The IMP model
106 @param representation PMI.representation.Representation object
107 (or list of them, for multi-state modeling)
108 @param root_hier Instead of passing Representation, pass the System hierarchy
109 @param monte_carlo_sample_objects Objects for MC sampling; a list of
110 structural components (generally the representation) that will
111 be moved and restraints with parameters that need to
113 For PMI2: just pass flat list of movers
114 @param molecular_dynamics_sample_objects Objects for MD sampling
115 For PMI2: just pass flat list of particles
116 @param output_objects A list of structural objects and restraints
117 that will be included in output (ie, statistics "stat" files). Any object
118 that provides a get_output() method can be used here. If None is passed
119 the macro will not write stat files.
120 @param rmf_output_objects A list of structural objects and restraints
121 that will be included in rmf. Any object
122 that provides a get_output() method can be used here.
123 @param crosslink_restraints List of cross-link restraints that will
124 be included in output RMF files (for visualization).
125 @param monte_carlo_temperature MC temp (may need to be optimized
126 based on post-sampling analysis)
127 @param simulated_annealing If True, perform simulated annealing
128 @param simulated_annealing_minimum_temperature Should generally be
129 the same as monte_carlo_temperature.
130 @param simulated_annealing_minimum_temperature_nframes Number of
131 frames to compute at minimum temperature.
132 @param simulated_annealing_maximum_temperature_nframes Number of
134 temps > simulated_annealing_maximum_temperature.
135 @param replica_exchange_minimum_temperature Low temp for REX; should
136 generally be the same as monte_carlo_temperature.
137 @param replica_exchange_maximum_temperature High temp for REX
138 @param replica_exchange_swap Boolean, enable disable temperature
140 @param num_sample_rounds Number of rounds of MC/MD per cycle
141 @param number_of_best_scoring_models Number of top-scoring PDB models
142 to keep around for analysis
143 @param monte_carlo_steps Number of MC steps per round
144 @param self_adaptive self adaptive scheme for monte carlo movers
145 @param molecular_dynamics_steps Number of MD steps per round
146 @param molecular_dynamics_max_time_step Max time step for MD
147 @param number_of_frames Number of REX frames to run
148 @param save_coordinates_mode string: how to save coordinates.
149 "lowest_temperature" (default) only the lowest temperatures is saved
150 "25th_score" all replicas whose score is below the 25th percentile
151 "50th_score" all replicas whose score is below the 50th percentile
152 "75th_score" all replicas whose score is below the 75th percentile
153 @param nframes_write_coordinates How often to write the coordinates
155 @param write_initial_rmf Write the initial configuration
156 @param global_output_directory Folder that will be created to house
158 @param test_mode Set to True to avoid writing any files, just test one frame.
165 self.output_objects = output_objects
166 self.rmf_output_objects=rmf_output_objects
167 self.representation = representation
169 if isinstance(representation, list):
170 self.is_multi_state =
True
171 self.root_hiers = [r.prot
for r
in representation]
172 self.vars[
"number_of_states"] = len(representation)
174 self.is_multi_state =
False
175 self.root_hier = representation.prot
176 self.vars[
"number_of_states"] = 1
178 and root_hier.get_name() ==
'System':
180 if self.output_objects
is not None:
182 if self.rmf_output_objects
is not None:
184 self.root_hier = root_hier
185 states = IMP.atom.get_by_type(root_hier,IMP.atom.STATE_TYPE)
186 self.vars[
"number_of_states"] = len(states)
188 self.root_hiers = states
189 self.is_multi_state =
True
191 self.root_hier = root_hier
192 self.is_multi_state =
False
194 raise Exception(
"Must provide representation or System hierarchy (root_hier)")
196 self._rmf_restraints = _RMFRestraints(model, crosslink_restraints)
197 self.em_object_for_rmf = em_object_for_rmf
198 self.monte_carlo_sample_objects = monte_carlo_sample_objects
199 self.vars[
"self_adaptive"]=self_adaptive
200 if sample_objects
is not None:
201 self.monte_carlo_sample_objects+=sample_objects
202 self.molecular_dynamics_sample_objects=molecular_dynamics_sample_objects
203 self.replica_exchange_object = replica_exchange_object
204 self.molecular_dynamics_max_time_step = molecular_dynamics_max_time_step
205 self.vars[
"monte_carlo_temperature"] = monte_carlo_temperature
207 "replica_exchange_minimum_temperature"] = replica_exchange_minimum_temperature
209 "replica_exchange_maximum_temperature"] = replica_exchange_maximum_temperature
210 self.vars[
"replica_exchange_swap"] = replica_exchange_swap
211 self.vars[
"simulated_annealing"]=\
213 self.vars[
"simulated_annealing_minimum_temperature"]=\
214 simulated_annealing_minimum_temperature
215 self.vars[
"simulated_annealing_maximum_temperature"]=\
216 simulated_annealing_maximum_temperature
217 self.vars[
"simulated_annealing_minimum_temperature_nframes"]=\
218 simulated_annealing_minimum_temperature_nframes
219 self.vars[
"simulated_annealing_maximum_temperature_nframes"]=\
220 simulated_annealing_maximum_temperature_nframes
222 self.vars[
"num_sample_rounds"] = num_sample_rounds
224 "number_of_best_scoring_models"] = number_of_best_scoring_models
225 self.vars[
"monte_carlo_steps"] = monte_carlo_steps
226 self.vars[
"molecular_dynamics_steps"]=molecular_dynamics_steps
227 self.vars[
"number_of_frames"] = number_of_frames
228 if not save_coordinates_mode
in [
"lowest_temperature",
"25th_score",
"50th_score",
"75th_score"]:
229 raise Exception(
"save_coordinates_mode has unrecognized value")
231 self.vars[
"save_coordinates_mode"] = save_coordinates_mode
232 self.vars[
"nframes_write_coordinates"] = nframes_write_coordinates
233 self.vars[
"write_initial_rmf"] = write_initial_rmf
234 self.vars[
"initial_rmf_name_suffix"] = initial_rmf_name_suffix
235 self.vars[
"best_pdb_name_suffix"] = best_pdb_name_suffix
236 self.vars[
"stat_file_name_suffix"] = stat_file_name_suffix
237 self.vars[
"do_clean_first"] = do_clean_first
238 self.vars[
"do_create_directories"] = do_create_directories
239 self.vars[
"global_output_directory"] = global_output_directory
240 self.vars[
"rmf_dir"] = rmf_dir
241 self.vars[
"best_pdb_dir"] = best_pdb_dir
242 self.vars[
"atomistic"] = atomistic
243 self.vars[
"replica_stat_file_suffix"] = replica_stat_file_suffix
244 self.vars[
"geometries"] =
None
245 self.test_mode = test_mode
248 if self.vars[
"geometries"]
is None:
249 self.vars[
"geometries"] = list(geometries)
251 self.vars[
"geometries"].extend(geometries)
254 print(
"ReplicaExchange0: it generates initial.*.rmf3, stat.*.out, rmfs/*.rmf3 for each replica ")
255 print(
"--- it stores the best scoring pdb models in pdbs/")
256 print(
"--- the stat.*.out and rmfs/*.rmf3 are saved only at the lowest temperature")
257 print(
"--- variables:")
258 keys = list(self.vars.keys())
261 print(
"------", v.ljust(30), self.vars[v])
263 def get_replica_exchange_object(self):
264 return self.replica_exchange_object
266 def _add_provenance(self, sampler_md, sampler_mc):
267 """Record details about the sampling in the IMP Hierarchies"""
268 if not self.is_multi_state
or self.pmi2:
269 output_hierarchies = [self.root_hier]
271 output_hierarchies = self.root_hiers
275 method =
"Molecular Dynamics"
276 iterations += self.vars[
"molecular_dynamics_steps"]
278 method =
"Hybrid MD/MC" if sampler_md
else "Monte Carlo"
279 iterations += self.vars[
"monte_carlo_steps"]
281 if iterations == 0
or self.vars[
"number_of_frames"] == 0:
283 iterations *= self.vars[
"num_sample_rounds"]
285 for h
in output_hierarchies:
286 pi = self.model.add_particle(
"sampling")
288 self.model, pi, method, self.vars[
"number_of_frames"],
290 p.set_number_of_replicas(
291 self.replica_exchange_object.get_number_of_replicas())
292 IMP.pmi.tools._add_pmi_provenance(h)
295 def execute_macro(self):
296 temp_index_factor = 100000.0
300 if self.monte_carlo_sample_objects
is not None:
301 print(
"Setting up MonteCarlo")
303 self.monte_carlo_sample_objects,
304 self.vars[
"monte_carlo_temperature"])
305 if self.vars[
"simulated_annealing"]:
306 tmin=self.vars[
"simulated_annealing_minimum_temperature"]
307 tmax=self.vars[
"simulated_annealing_maximum_temperature"]
308 nfmin=self.vars[
"simulated_annealing_minimum_temperature_nframes"]
309 nfmax=self.vars[
"simulated_annealing_maximum_temperature_nframes"]
310 sampler_mc.set_simulated_annealing(tmin,tmax,nfmin,nfmax)
311 if self.vars[
"self_adaptive"]:
312 sampler_mc.set_self_adaptive(isselfadaptive=self.vars[
"self_adaptive"])
313 if self.output_objects
is not None:
314 self.output_objects.append(sampler_mc)
315 if self.rmf_output_objects
is not None:
316 self.rmf_output_objects.append(sampler_mc)
317 samplers.append(sampler_mc)
320 if self.molecular_dynamics_sample_objects
is not None:
321 print(
"Setting up MolecularDynamics")
323 self.molecular_dynamics_sample_objects,
324 self.vars[
"monte_carlo_temperature"],
325 maximum_time_step=self.molecular_dynamics_max_time_step)
326 if self.vars[
"simulated_annealing"]:
327 tmin=self.vars[
"simulated_annealing_minimum_temperature"]
328 tmax=self.vars[
"simulated_annealing_maximum_temperature"]
329 nfmin=self.vars[
"simulated_annealing_minimum_temperature_nframes"]
330 nfmax=self.vars[
"simulated_annealing_maximum_temperature_nframes"]
331 sampler_md.set_simulated_annealing(tmin,tmax,nfmin,nfmax)
332 if self.output_objects
is not None:
333 self.output_objects.append(sampler_md)
334 if self.rmf_output_objects
is not None:
335 self.rmf_output_objects.append(sampler_md)
336 samplers.append(sampler_md)
339 print(
"Setting up ReplicaExchange")
342 "replica_exchange_minimum_temperature"],
344 "replica_exchange_maximum_temperature"],
346 replica_exchange_object=self.replica_exchange_object)
347 self.replica_exchange_object = rex.rem
349 myindex = rex.get_my_index()
350 if self.output_objects
is not None:
351 self.output_objects.append(rex)
352 if self.rmf_output_objects
is not None:
353 self.rmf_output_objects.append(rex)
357 min_temp_index = int(min(rex.get_temperatures()) * temp_index_factor)
361 globaldir = self.vars[
"global_output_directory"] +
"/"
362 rmf_dir = globaldir + self.vars[
"rmf_dir"]
363 pdb_dir = globaldir + self.vars[
"best_pdb_dir"]
365 if not self.test_mode:
366 if self.vars[
"do_clean_first"]:
369 if self.vars[
"do_create_directories"]:
372 os.makedirs(globaldir)
380 if not self.is_multi_state:
386 for n
in range(self.vars[
"number_of_states"]):
388 os.makedirs(pdb_dir +
"/" + str(n))
395 if self.output_objects
is not None:
396 self.output_objects.append(sw)
397 if self.rmf_output_objects
is not None:
398 self.rmf_output_objects.append(sw)
400 print(
"Setting up stat file")
402 low_temp_stat_file = globaldir + \
403 self.vars[
"stat_file_name_suffix"] +
"." + str(myindex) +
".out"
406 if not self.test_mode:
409 if not self.test_mode:
410 if self.output_objects
is not None:
411 output.init_stat2(low_temp_stat_file,
413 extralabels=[
"rmf_file",
"rmf_frame_index"])
415 print(
"Stat file writing is disabled")
417 if self.rmf_output_objects
is not None:
418 print(
"Stat info being written in the rmf file")
420 print(
"Setting up replica stat file")
421 replica_stat_file = globaldir + \
422 self.vars[
"replica_stat_file_suffix"] +
"." + str(myindex) +
".out"
423 if not self.test_mode:
424 output.init_stat2(replica_stat_file, [rex], extralabels=[
"score"])
426 if not self.test_mode:
427 print(
"Setting up best pdb files")
428 if not self.is_multi_state:
429 if self.vars[
"number_of_best_scoring_models"] > 0:
430 output.init_pdb_best_scoring(pdb_dir +
"/" +
431 self.vars[
"best_pdb_name_suffix"],
434 "number_of_best_scoring_models"],
435 replica_exchange=
True)
436 output.write_psf(pdb_dir +
"/" +
"model.psf",pdb_dir +
"/" +
437 self.vars[
"best_pdb_name_suffix"]+
".0.pdb")
439 if self.vars[
"number_of_best_scoring_models"] > 0:
440 for n
in range(self.vars[
"number_of_states"]):
441 output.init_pdb_best_scoring(pdb_dir +
"/" + str(n) +
"/" +
442 self.vars[
"best_pdb_name_suffix"],
445 "number_of_best_scoring_models"],
446 replica_exchange=
True)
447 output.write_psf(pdb_dir +
"/" + str(n) +
"/" +
"model.psf",pdb_dir +
"/" + str(n) +
"/" +
448 self.vars[
"best_pdb_name_suffix"]+
".0.pdb")
451 if self.em_object_for_rmf
is not None:
452 if not self.is_multi_state
or self.pmi2:
453 output_hierarchies = [
455 self.em_object_for_rmf.get_density_as_hierarchy(
458 output_hierarchies = self.root_hiers
459 output_hierarchies.append(
460 self.em_object_for_rmf.get_density_as_hierarchy())
462 if not self.is_multi_state
or self.pmi2:
463 output_hierarchies = [self.root_hier]
465 output_hierarchies = self.root_hiers
468 if not self.test_mode:
469 print(
"Setting up and writing initial rmf coordinate file")
470 init_suffix = globaldir + self.vars[
"initial_rmf_name_suffix"]
471 output.init_rmf(init_suffix +
"." + str(myindex) +
".rmf3",
472 output_hierarchies,listofobjects=self.rmf_output_objects)
473 if self._rmf_restraints:
474 output.add_restraints_to_rmf(
475 init_suffix +
"." + str(myindex) +
".rmf3",
476 self._rmf_restraints)
477 output.write_rmf(init_suffix +
"." + str(myindex) +
".rmf3")
478 output.close_rmf(init_suffix +
"." + str(myindex) +
".rmf3")
482 if not self.test_mode:
483 mpivs=IMP.pmi.samplers.MPI_values(self.replica_exchange_object)
487 self._add_provenance(sampler_md, sampler_mc)
489 if not self.test_mode:
490 print(
"Setting up production rmf files")
491 rmfname = rmf_dir +
"/" + str(myindex) +
".rmf3"
492 output.init_rmf(rmfname, output_hierarchies, geometries=self.vars[
"geometries"],
493 listofobjects = self.rmf_output_objects)
495 if self._rmf_restraints:
496 output.add_restraints_to_rmf(rmfname, self._rmf_restraints)
498 ntimes_at_low_temp = 0
502 self.replica_exchange_object.set_was_used(
True)
503 nframes = self.vars[
"number_of_frames"]
506 for i
in range(nframes):
510 for nr
in range(self.vars[
"num_sample_rounds"]):
511 if sampler_md
is not None:
513 self.vars[
"molecular_dynamics_steps"])
514 if sampler_mc
is not None:
515 sampler_mc.optimize(self.vars[
"monte_carlo_steps"])
517 self.model).evaluate(
False)
518 mpivs.set_value(
"score",score)
519 output.set_output_entry(
"score", score)
523 my_temp_index = int(rex.get_my_temp() * temp_index_factor)
525 if self.vars[
"save_coordinates_mode"] ==
"lowest_temperature":
526 save_frame=(min_temp_index == my_temp_index)
527 elif self.vars[
"save_coordinates_mode"] ==
"25th_score":
528 score_perc=mpivs.get_percentile(
"score")
529 save_frame=(score_perc*100.0<=25.0)
530 elif self.vars[
"save_coordinates_mode"] ==
"50th_score":
531 score_perc=mpivs.get_percentile(
"score")
532 save_frame=(score_perc*100.0<=50.0)
533 elif self.vars[
"save_coordinates_mode"] ==
"75th_score":
534 score_perc=mpivs.get_percentile(
"score")
535 save_frame=(score_perc*100.0<=75.0)
538 if save_frame
or not self.test_mode:
542 print(
"--- frame %s score %s " % (str(i), str(score)))
544 if not self.test_mode:
545 if i % self.vars[
"nframes_write_coordinates"]==0:
546 print(
'--- writing coordinates')
547 if self.vars[
"number_of_best_scoring_models"] > 0:
548 output.write_pdb_best_scoring(score)
549 output.write_rmf(rmfname)
550 output.set_output_entry(
"rmf_file", rmfname)
551 output.set_output_entry(
"rmf_frame_index", ntimes_at_low_temp)
553 output.set_output_entry(
"rmf_file", rmfname)
554 output.set_output_entry(
"rmf_frame_index",
'-1')
555 if self.output_objects
is not None:
556 output.write_stat2(low_temp_stat_file)
557 ntimes_at_low_temp += 1
559 if not self.test_mode:
560 output.write_stat2(replica_stat_file)
561 if self.vars[
"replica_exchange_swap"]:
562 rex.swap_temp(i, score)
563 for p, state
in IMP.pmi.tools._all_protocol_outputs(
564 [self.representation],
565 self.root_hier
if self.pmi2
else None):
566 p.add_replica_exchange(state, self)
568 if not self.test_mode:
569 print(
"closing production rmf files")
570 output.close_rmf(rmfname)
576 """A macro to build a IMP::pmi::topology::System based on a TopologyReader object.
577 Easily create multi-state systems by calling this macro
578 repeatedly with different TopologyReader objects!
579 A useful function is get_molecules() which returns the PMI Molecules grouped by state
580 as a dictionary with key = (molecule name), value = IMP.pmi.topology.Molecule
581 Quick multi-state system:
584 reader1 = IMP.pmi.topology.TopologyReader(tfile1)
585 reader2 = IMP.pmi.topology.TopologyReader(tfile2)
586 bs = IMP.pmi.macros.BuildSystem(model)
587 bs.add_state(reader1)
588 bs.add_state(reader2)
589 bs.execute_macro() # build everything including degrees of freedom
590 IMP.atom.show_molecular_hierarchy(bs.get_hierarchy())
591 ### now you have a two state system, you add restraints etc
593 @note The "domain name" entry of the topology reader is not used.
594 All molecules are set up by the component name, but split into rigid bodies
599 sequence_connectivity_scale=4.0,
600 force_create_gmm_files=
False,
603 @param model An IMP Model
604 @param sequence_connectivity_scale For scaling the connectivity restraint
605 @param force_create_gmm_files If True, will sample and create GMMs
606 no matter what. If False, will only sample if the
607 files don't exist. If number of Gaussians is zero, won't
609 @param resolutions The resolutions to build for structured regions
614 self._domain_res = []
616 self.force_create_gmm_files = force_create_gmm_files
617 self.resolutions = resolutions
621 keep_chain_id=
False, fasta_name_map=
None):
622 """Add a state using the topology info in a IMP::pmi::topology::TopologyReader object.
623 When you are done adding states, call execute_macro()
624 @param reader The TopologyReader object
625 @param fasta_name_map dictionary for converting protein names found in the fasta file
627 state = self.system.create_state()
628 self._readers.append(reader)
629 these_domain_res = {}
631 chain_ids = string.ascii_uppercase+string.ascii_lowercase+
'0123456789'
636 for molname
in reader.get_molecules():
637 copies = reader.get_molecules()[molname].domains
638 for nc,copyname
in enumerate(copies):
639 print(
"BuildSystem.add_state: setting up molecule %s copy number %s" % (molname,str(nc)))
640 copy = copies[copyname]
642 if keep_chain_id ==
True:
643 chain_id = copy[0].chain
645 chain_id = chain_ids[numchain]
648 fasta_flag=copy[0].fasta_flag
649 if fasta_flag ==
"RNA" or fasta_flag ==
"DNA": is_nucleic=
True
651 print(
"BuildSystem.add_state: molecule %s sequence has %s residues" % (molname,len(seq)))
652 orig_mol = state.create_molecule(molname,
658 print(
"BuildSystem.add_state: creating a copy for molecule %s" % molname)
659 mol = orig_mol.create_copy(chain_id)
662 for domainnumber,domain
in enumerate(copy):
663 print(
"BuildSystem.add_state: ---- setting up domain %s of molecule %s" % (domainnumber,molname))
666 these_domains[domain.get_unique_name()] = domain
667 if domain.residue_range==[]
or domain.residue_range
is None:
668 domain_res = mol.get_residues()
670 start = domain.residue_range[0]+domain.pdb_offset
671 if domain.residue_range[1]==
'END':
672 end = len(mol.sequence)
674 end = domain.residue_range[1]+domain.pdb_offset
675 domain_res = mol.residue_range(start-1,end-1)
676 print(
"BuildSystem.add_state: -------- domain %s of molecule %s extends from residue %s to residue %s " % (domainnumber,molname,start,end))
677 if domain.pdb_file==
"BEADS":
678 print(
"BuildSystem.add_state: -------- domain %s of molecule %s represented by BEADS " % (domainnumber,molname))
679 mol.add_representation(
681 resolutions=[domain.bead_size],
682 setup_particles_as_densities=(
683 domain.em_residues_per_gaussian!=0),
684 color = domain.color)
685 these_domain_res[domain.get_unique_name()] = (set(),domain_res)
686 elif domain.pdb_file==
"IDEAL_HELIX":
687 print(
"BuildSystem.add_state: -------- domain %s of molecule %s represented by IDEAL_HELIX " % (domainnumber,molname))
688 mol.add_representation(
690 resolutions=self.resolutions,
692 density_residues_per_component=domain.em_residues_per_gaussian,
693 density_prefix=domain.density_prefix,
694 density_force_compute=self.force_create_gmm_files,
695 color = domain.color)
696 these_domain_res[domain.get_unique_name()] = (domain_res,set())
698 print(
"BuildSystem.add_state: -------- domain %s of molecule %s represented by pdb file %s " % (domainnumber,molname,domain.pdb_file))
699 domain_atomic = mol.add_structure(domain.pdb_file,
701 domain.residue_range,
704 domain_non_atomic = domain_res - domain_atomic
705 if not domain.em_residues_per_gaussian:
706 mol.add_representation(domain_atomic,
707 resolutions=self.resolutions,
708 color = domain.color)
709 if len(domain_non_atomic)>0:
710 mol.add_representation(domain_non_atomic,
711 resolutions=[domain.bead_size],
712 color = domain.color)
714 print(
"BuildSystem.add_state: -------- domain %s of molecule %s represented by gaussians " % (domainnumber,molname))
715 mol.add_representation(
717 resolutions=self.resolutions,
718 density_residues_per_component=domain.em_residues_per_gaussian,
719 density_prefix=domain.density_prefix,
720 density_force_compute=self.force_create_gmm_files,
721 color = domain.color)
722 if len(domain_non_atomic)>0:
723 mol.add_representation(domain_non_atomic,
724 resolutions=[domain.bead_size],
725 setup_particles_as_densities=
True,
726 color = domain.color)
727 these_domain_res[domain.get_unique_name()] = (
728 domain_atomic,domain_non_atomic)
729 self._domain_res.append(these_domain_res)
730 self._domains.append(these_domains)
731 print(
'BuildSystem.add_state: State',len(self.system.states),
'added')
734 """Return list of all molecules grouped by state.
735 For each state, it's a dictionary of Molecules where key is the molecule name
737 return [s.get_molecules()
for s
in self.system.get_states()]
739 def get_molecule(self,molname,copy_index=0,state_index=0):
740 return self.system.get_states()[state_index].
get_molecules()[molname][copy_index]
742 def execute_macro(self, max_rb_trans=4.0, max_rb_rot=0.04, max_bead_trans=4.0, max_srb_trans=4.0,max_srb_rot=0.04):
743 """Builds representations and sets up degrees of freedom"""
744 print(
"BuildSystem.execute_macro: building representations")
745 self.root_hier = self.system.build()
747 print(
"BuildSystem.execute_macro: setting up degrees of freedom")
749 for nstate,reader
in enumerate(self._readers):
750 rbs = reader.get_rigid_bodies()
751 srbs = reader.get_super_rigid_bodies()
752 csrbs = reader.get_chains_of_super_rigid_bodies()
755 domains_in_rbs = set()
757 print(
"BuildSystem.execute_macro: -------- building rigid body %s" % (str(rblist)))
758 all_res = IMP.pmi.tools.OrderedSet()
759 bead_res = IMP.pmi.tools.OrderedSet()
761 domain = self._domains[nstate][dname]
762 print(
"BuildSystem.execute_macro: -------- adding %s" % (str(dname )))
763 all_res|=self._domain_res[nstate][dname][0]
764 bead_res|=self._domain_res[nstate][dname][1]
765 domains_in_rbs.add(dname)
767 print(
"BuildSystem.execute_macro: -------- \
768 creating rigid body with max_trans %s \
769 max_rot %s non_rigid_max_trans %s" \
770 % (str(max_rb_trans),str(max_rb_rot),str(max_bead_trans)))
771 self.dof.create_rigid_body(all_res,
772 nonrigid_parts=bead_res,
773 max_trans=max_rb_trans,
775 nonrigid_max_trans=max_bead_trans)
778 for dname, domain
in self._domains[nstate].items():
779 if dname
not in domains_in_rbs:
780 if domain.pdb_file !=
"BEADS":
782 "Making %s flexible. This may distort the "
783 "structure; consider making it rigid" % dname,
785 self.dof.create_flexible_beads(
786 self._domain_res[nstate][dname][1],
787 max_trans=max_bead_trans)
791 print(
"BuildSystem.execute_macro: -------- building super rigid body %s" % (str(srblist)))
792 all_res = IMP.pmi.tools.OrderedSet()
793 for dname
in srblist:
794 print(
"BuildSystem.execute_macro: -------- adding %s" % (str(dname )))
795 all_res|=self._domain_res[nstate][dname][0]
796 all_res|=self._domain_res[nstate][dname][1]
798 print(
"BuildSystem.execute_macro: -------- \
799 creating super rigid body with max_trans %s max_rot %s " \
800 % (str(max_srb_trans),str(max_srb_rot)))
801 self.dof.create_super_rigid_body(all_res,max_trans=max_srb_trans,max_rot=max_srb_rot)
804 for csrblist
in csrbs:
805 all_res = IMP.pmi.tools.OrderedSet()
806 for dname
in csrblist:
807 all_res|=self._domain_res[nstate][dname][0]
808 all_res|=self._domain_res[nstate][dname][1]
809 all_res = list(all_res)
810 all_res.sort(key=
lambda r:r.get_index())
812 self.dof.create_main_chain_mover(all_res)
813 return self.root_hier,self.dof
818 """A macro for running all the basic operations of analysis.
819 Includes clustering, precision analysis, and making ensemble density maps.
820 A number of plots are also supported.
823 merge_directories=[
"./"],
824 stat_file_name_suffix=
"stat",
825 best_pdb_name_suffix=
"model",
827 do_create_directories=
True,
828 global_output_directory=
"output/",
829 replica_stat_file_suffix=
"stat_replica",
830 global_analysis_result_directory=
"./analysis/",
833 @param model The IMP model
834 @param stat_file_name_suffix
835 @param merge_directories The directories containing output files
836 @param best_pdb_name_suffix
837 @param do_clean_first
838 @param do_create_directories
839 @param global_output_directory Where everything is
840 @param replica_stat_file_suffix
841 @param global_analysis_result_directory
842 @param test_mode If True, nothing is changed on disk
846 from mpi4py
import MPI
847 self.comm = MPI.COMM_WORLD
848 self.rank = self.comm.Get_rank()
849 self.number_of_processes = self.comm.size
852 self.number_of_processes = 1
854 self.test_mode = test_mode
855 self._protocol_output = []
856 self.cluster_obj =
None
858 stat_dir = global_output_directory
861 for rd
in merge_directories:
862 stat_files = glob.glob(os.path.join(rd,stat_dir,
"stat.*.out"))
863 if len(stat_files)==0:
864 warnings.warn(
"no stat files found in %s"
865 % os.path.join(rd, stat_dir),
867 self.stat_files += stat_files
870 """Capture details of the modeling protocol.
871 @param p an instance of IMP.pmi.output.ProtocolOutput or a subclass.
874 self._protocol_output.append((p, p._last_state))
877 score_key=
"SimplifiedModel_Total_Score_None",
878 rmf_file_key=
"rmf_file",
879 rmf_file_frame_key=
"rmf_frame_index",
882 nframes_trajectory=10000):
883 """ Get a trajectory of the modeling run, for generating demonstrative movies
884 @param score_key The score for ranking models
885 @param rmf_file_key Key pointing to RMF filename
886 @param rmf_file_frame_key Key pointing to RMF frame number
887 @param outputdir The local output directory used in the run
888 @param get_every Extract every nth frame
889 @param nframes_trajectory Total number of frames of the trajectory
891 from operator
import itemgetter
899 rmf_file_list=trajectory_models[0]
900 rmf_file_frame_list=trajectory_models[1]
901 score_list=list(map(float, trajectory_models[2]))
903 max_score=max(score_list)
904 min_score=min(score_list)
907 bins=[(max_score-min_score)*math.exp(-float(i))+min_score
for i
in range(nframes_trajectory)]
908 binned_scores=[
None]*nframes_trajectory
909 binned_model_indexes=[-1]*nframes_trajectory
911 for model_index,s
in enumerate(score_list):
912 bins_score_diffs=[abs(s-b)
for b
in bins]
913 bin_index=min(enumerate(bins_score_diffs), key=itemgetter(1))[0]
914 if binned_scores[bin_index]==
None:
915 binned_scores[bin_index]=s
916 binned_model_indexes[bin_index]=model_index
918 old_diff=abs(binned_scores[bin_index]-bins[bin_index])
919 new_diff=abs(s-bins[bin_index])
920 if new_diff < old_diff:
921 binned_scores[bin_index]=s
922 binned_model_indexes[bin_index]=model_index
925 print(binned_model_indexes)
928 def _expand_ambiguity(self,prot,d):
929 """If using PMI2, expand the dictionary to include copies as ambiguous options
930 This also keeps the states separate.
935 if '..' in key
or (isinstance(val, tuple)
and len(val) >= 3):
938 states = IMP.atom.get_by_type(prot,IMP.atom.STATE_TYPE)
939 if isinstance(val, tuple):
947 for nst
in range(len(states)):
949 copies = sel.get_selected_particles(with_representation=
False)
951 for nc
in range(len(copies)):
953 newdict[
'%s.%i..%i'%(name,nst,nc)] = (start,stop,name,nc,nst)
955 newdict[
'%s..%i'%(name,nc)] = (start,stop,name,nc,nst)
962 score_key=
"SimplifiedModel_Total_Score_None",
963 rmf_file_key=
"rmf_file",
964 rmf_file_frame_key=
"rmf_frame_index",
969 alignment_components=
None,
970 number_of_best_scoring_models=10,
971 rmsd_calculation_components=
None,
972 distance_matrix_file=
'distances.mat',
973 load_distance_matrix_file=
False,
974 skip_clustering=
False,
975 number_of_clusters=1,
977 exit_after_display=
True,
979 first_and_last_frames=
None,
980 density_custom_ranges=
None,
981 write_pdb_with_centered_coordinates=
False,
983 """ Get the best scoring models, compute a distance matrix, cluster them, and create density maps.
984 Tuple format: "molname" just the molecule, or (start,stop,molname,copy_num(optional),state_num(optional)
985 Can pass None for copy or state to ignore that field.
986 If you don't pass a specific copy number
987 @param score_key The score for ranking models. Use "Total_Score"
988 instead of default for PMI2.
989 @param rmf_file_key Key pointing to RMF filename
990 @param rmf_file_frame_key Key pointing to RMF frame number
991 @param state_number State number to analyze
992 @param prefiltervalue Only include frames where the
993 score key is below this value
994 @param feature_keys Keywords for which you want to
995 calculate average, medians, etc.
996 If you pass "Keyname" it'll include everything that matches "*Keyname*"
997 @param outputdir The local output directory used in the run
998 @param alignment_components Dictionary with keys=groupname, values are tuples
999 for aligning the structures e.g. {"Rpb1": (20,100,"Rpb1"),"Rpb2":"Rpb2"}
1000 @param number_of_best_scoring_models Num models to keep per run
1001 @param rmsd_calculation_components For calculating RMSD
1002 (same format as alignment_components)
1003 @param distance_matrix_file Where to store/read the distance matrix
1004 @param load_distance_matrix_file Try to load the distance matrix file
1005 @param skip_clustering Just extract the best scoring models
1007 @param number_of_clusters Number of k-means clusters
1008 @param display_plot Display the distance matrix
1009 @param exit_after_display Exit after displaying distance matrix
1010 @param get_every Extract every nth frame
1011 @param first_and_last_frames A tuple with the first and last frames to be
1012 analyzed. Values are percentages!
1013 Default: get all frames
1014 @param density_custom_ranges For density calculation
1015 (same format as alignment_components)
1016 @param write_pdb_with_centered_coordinates
1017 @param voxel_size Used for the density output
1021 self._outputdir = os.path.abspath(outputdir)
1022 self._number_of_clusters = number_of_clusters
1023 for p, state
in self._protocol_output:
1024 p.add_replica_exchange_analysis(state, self, density_custom_ranges)
1035 if not load_distance_matrix_file:
1036 if len(self.stat_files)==0: print(
"ERROR: no stat file found in the given path");
return
1037 my_stat_files = IMP.pmi.tools.chunk_list_into_segments(
1038 self.stat_files,self.number_of_processes)[self.rank]
1042 orig_score_key = score_key
1043 if score_key
not in po.get_keys():
1044 if 'Total_Score' in po.get_keys():
1045 score_key =
'Total_Score'
1047 "Using 'Total_Score' instead of "
1048 "'SimplifiedModel_Total_Score_None' for the score key",
1050 for k
in [orig_score_key,score_key,rmf_file_key,rmf_file_frame_key]:
1051 if k
in feature_keys:
1053 "no need to pass " + k +
" to feature_keys.",
1055 feature_keys.remove(k)
1063 get_every, provenance=prov)
1064 rmf_file_list=best_models[0]
1065 rmf_file_frame_list=best_models[1]
1066 score_list=best_models[2]
1067 feature_keyword_list_dict=best_models[3]
1073 if self.number_of_processes > 1:
1077 rmf_file_frame_list)
1078 for k
in feature_keyword_list_dict:
1080 feature_keyword_list_dict[k])
1083 score_rmf_tuples = list(zip(score_list,
1085 rmf_file_frame_list,
1086 list(range(len(score_list)))))
1088 if density_custom_ranges:
1089 for k
in density_custom_ranges:
1090 if not isinstance(density_custom_ranges[k], list):
1091 raise Exception(
"Density custom ranges: values must be lists of tuples")
1094 if first_and_last_frames
is not None:
1095 nframes = len(score_rmf_tuples)
1096 first_frame = int(first_and_last_frames[0] * nframes)
1097 last_frame = int(first_and_last_frames[1] * nframes)
1098 if last_frame > len(score_rmf_tuples):
1100 score_rmf_tuples = score_rmf_tuples[first_frame:last_frame]
1103 best_score_rmf_tuples = sorted(score_rmf_tuples,
1104 key=
lambda x: float(x[0]))[:number_of_best_scoring_models]
1105 best_score_rmf_tuples=[t+(n,)
for n,t
in enumerate(best_score_rmf_tuples)]
1107 prov.append(IMP.pmi.io.FilterProvenance(
"Best scoring",
1108 0, number_of_best_scoring_models))
1110 best_score_feature_keyword_list_dict = defaultdict(list)
1111 for tpl
in best_score_rmf_tuples:
1113 for f
in feature_keyword_list_dict:
1114 best_score_feature_keyword_list_dict[f].append(
1115 feature_keyword_list_dict[f][index])
1116 my_best_score_rmf_tuples = IMP.pmi.tools.chunk_list_into_segments(
1117 best_score_rmf_tuples,
1118 self.number_of_processes)[self.rank]
1121 prot_ahead = IMP.pmi.analysis.get_hiers_from_rmf(self.model,
1123 my_best_score_rmf_tuples[0][1])[0]
1125 if rmsd_calculation_components
is not None:
1126 tmp = self._expand_ambiguity(prot_ahead,rmsd_calculation_components)
1127 if tmp!=rmsd_calculation_components:
1128 print(
'Detected ambiguity, expand rmsd components to',tmp)
1129 rmsd_calculation_components = tmp
1130 if alignment_components
is not None:
1131 tmp = self._expand_ambiguity(prot_ahead,alignment_components)
1132 if tmp!=alignment_components:
1133 print(
'Detected ambiguity, expand alignment components to',tmp)
1134 alignment_components = tmp
1140 my_best_score_rmf_tuples[0],
1141 rmsd_calculation_components,
1142 state_number=state_number)
1144 my_best_score_rmf_tuples,
1145 alignment_components,
1146 rmsd_calculation_components,
1147 state_number=state_number)
1152 all_coordinates=got_coords[0]
1153 alignment_coordinates=got_coords[1]
1154 rmsd_coordinates=got_coords[2]
1155 rmf_file_name_index_dict=got_coords[3]
1156 all_rmf_file_names=got_coords[4]
1162 if density_custom_ranges:
1164 density_custom_ranges,
1167 dircluster=os.path.join(outputdir,
"all_models."+str(n))
1173 os.mkdir(dircluster)
1176 clusstat=open(os.path.join(dircluster,
"stat."+str(self.rank)+
".out"),
"w")
1177 for cnt,tpl
in enumerate(my_best_score_rmf_tuples):
1179 rmf_frame_number=tpl[2]
1182 for key
in best_score_feature_keyword_list_dict:
1183 tmp_dict[key]=best_score_feature_keyword_list_dict[key][index]
1186 prots,rs = IMP.pmi.analysis.get_hiers_and_restraints_from_rmf(
1191 linking_successful=IMP.pmi.analysis.link_hiers_and_restraints_to_rmf(
1197 if not linking_successful:
1204 states = IMP.atom.get_by_type(prots[0],IMP.atom.STATE_TYPE)
1205 prot = states[state_number]
1207 prot = prots[state_number]
1212 coords_f1=alignment_coordinates[cnt]
1214 coords_f2=alignment_coordinates[cnt]
1217 transformation = Ali.align()[1]
1239 out_pdb_fn=os.path.join(dircluster,str(cnt)+
"."+str(self.rank)+
".pdb")
1240 out_rmf_fn=os.path.join(dircluster,str(cnt)+
"."+str(self.rank)+
".rmf3")
1241 o.init_pdb(out_pdb_fn,prot)
1242 o.write_pdb(out_pdb_fn,
1243 translate_to_geometric_center=write_pdb_with_centered_coordinates)
1245 tmp_dict[
"local_pdb_file_name"]=os.path.basename(out_pdb_fn)
1246 tmp_dict[
"rmf_file_full_path"]=rmf_name
1247 tmp_dict[
"local_rmf_file_name"]=os.path.basename(out_rmf_fn)
1248 tmp_dict[
"local_rmf_frame_number"]=0
1250 clusstat.write(str(tmp_dict)+
"\n")
1255 h.set_name(
"System")
1257 o.init_rmf(out_rmf_fn, [h], rs)
1259 o.init_rmf(out_rmf_fn, [prot],rs)
1261 o.write_rmf(out_rmf_fn)
1262 o.close_rmf(out_rmf_fn)
1264 if density_custom_ranges:
1265 DensModule.add_subunits_density(prot)
1267 if density_custom_ranges:
1268 DensModule.write_mrc(path=dircluster)
1275 if self.number_of_processes > 1:
1281 rmf_file_name_index_dict)
1283 alignment_coordinates)
1290 [best_score_feature_keyword_list_dict,
1291 rmf_file_name_index_dict],
1297 print(
"setup clustering class")
1300 for n, model_coordinate_dict
in enumerate(all_coordinates):
1301 template_coordinate_dict = {}
1303 if alignment_components
is not None and len(self.cluster_obj.all_coords) == 0:
1305 self.cluster_obj.set_template(alignment_coordinates[n])
1306 self.cluster_obj.fill(all_rmf_file_names[n], rmsd_coordinates[n])
1307 print(
"Global calculating the distance matrix")
1310 self.cluster_obj.dist_matrix()
1314 self.cluster_obj.do_cluster(number_of_clusters)
1317 self.cluster_obj.plot_matrix(figurename=os.path.join(outputdir,
'dist_matrix.pdf'))
1318 if exit_after_display:
1320 self.cluster_obj.save_distance_matrix_file(file_name=distance_matrix_file)
1327 print(
"setup clustering class")
1329 self.cluster_obj.load_distance_matrix_file(file_name=distance_matrix_file)
1330 print(
"clustering with %s clusters" % str(number_of_clusters))
1331 self.cluster_obj.do_cluster(number_of_clusters)
1332 [best_score_feature_keyword_list_dict,
1333 rmf_file_name_index_dict] = self.load_objects(
".macro.pkl")
1336 self.cluster_obj.plot_matrix(figurename=os.path.join(outputdir,
'dist_matrix.pdf'))
1337 if exit_after_display:
1339 if self.number_of_processes > 1:
1347 print(self.cluster_obj.get_cluster_labels())
1348 for n, cl
in enumerate(self.cluster_obj.get_cluster_labels()):
1349 print(
"rank %s " % str(self.rank))
1350 print(
"cluster %s " % str(n))
1351 print(
"cluster label %s " % str(cl))
1352 print(self.cluster_obj.get_cluster_label_names(cl))
1353 cluster_size = len(self.cluster_obj.get_cluster_label_names(cl))
1354 cluster_prov = prov + \
1355 [IMP.pmi.io.ClusterProvenance(cluster_size)]
1358 if density_custom_ranges:
1360 density_custom_ranges,
1363 dircluster = outputdir +
"/cluster." + str(n) +
"/"
1365 os.mkdir(dircluster)
1369 rmsd_dict = {
"AVERAGE_RMSD":
1370 str(self.cluster_obj.get_cluster_label_average_rmsd(cl))}
1371 clusstat = open(dircluster +
"stat.out",
"w")
1372 for k, structure_name
in enumerate(self.cluster_obj.get_cluster_label_names(cl)):
1375 tmp_dict.update(rmsd_dict)
1376 index = rmf_file_name_index_dict[structure_name]
1377 for key
in best_score_feature_keyword_list_dict:
1379 key] = best_score_feature_keyword_list_dict[
1385 rmf_name = structure_name.split(
"|")[0]
1386 rmf_frame_number = int(structure_name.split(
"|")[1])
1387 clusstat.write(str(tmp_dict) +
"\n")
1391 prots,rs = IMP.pmi.analysis.get_hiers_and_restraints_from_rmf(
1396 linking_successful = IMP.pmi.analysis.link_hiers_and_restraints_to_rmf(
1402 if not linking_successful:
1408 states = IMP.atom.get_by_type(prots[0],IMP.atom.STATE_TYPE)
1409 prot = states[state_number]
1411 prot = prots[state_number]
1417 model_index = self.cluster_obj.get_model_index_from_name(
1419 transformation = self.cluster_obj.get_transformation_to_first_member(
1439 if density_custom_ranges:
1440 DensModule.add_subunits_density(prot)
1445 o.init_pdb(dircluster + str(k) +
".pdb", prot)
1446 o.write_pdb(dircluster + str(k) +
".pdb")
1451 h.set_name(
"System")
1453 o.init_rmf(dircluster + str(k) +
".rmf3", [h], rs)
1455 o.init_rmf(dircluster + str(k) +
".rmf3", [prot],rs)
1456 o.write_rmf(dircluster + str(k) +
".rmf3")
1457 o.close_rmf(dircluster + str(k) +
".rmf3")
1462 if density_custom_ranges:
1463 DensModule.write_mrc(path=dircluster)
1466 if self.number_of_processes>1:
1469 def get_cluster_rmsd(self,cluster_num):
1470 if self.cluster_obj
is None:
1472 return self.cluster_obj.get_cluster_label_average_rmsd(cluster_num)
1474 def save_objects(self, objects, file_name):
1476 outf = open(file_name,
'wb')
1477 pickle.dump(objects, outf)
1480 def load_objects(self, file_name):
1482 inputf = open(file_name,
'rb')
1483 objects = pickle.load(inputf)
1490 This class contains analysis utilities to investigate ReplicaExchange results.
1503 Construction of the Class.
1504 @param model IMP.Model()
1505 @param stat_files list of string. Can be ascii stat files, rmf files names
1506 @param best_models Integer. Number of best scoring models, if None: all models will be read
1507 @param alignment boolean (Default=True). Align before computing the rmsd.
1511 self.best_models=best_models
1524 self.clusters.append(c)
1525 for n0
in range(len(self.stath0)):
1527 self.pairwise_rmsd={}
1528 self.pairwise_molecular_assignment={}
1529 self.alignment=alignment
1530 self.symmetric_molecules={}
1531 self.issymmetricsel={}
1533 self.molcopydict0=IMP.pmi.tools.get_molecules_dictionary_by_copy(
IMP.atom.get_leaves(self.stath0))
1534 self.molcopydict1=IMP.pmi.tools.get_molecules_dictionary_by_copy(
IMP.atom.get_leaves(self.stath1))
1538 Setup the selection onto which the rmsd is computed
1539 @param kwargs use IMP.atom.Selection keywords
1547 Store names of symmetric molecules
1549 self.symmetric_molecules[molecule_name]=0
1554 Setup the selection onto which the alignment is computed
1555 @param kwargs use IMP.atom.Selection keywords
1563 def clean_clusters(self):
1564 for c
in self.clusters: del c
1568 def cluster(self, rmsd_cutoff=10, metric=IMP.atom.get_rmsd):
1570 Cluster the models based on RMSD.
1571 @param rmsd_cutoff Float the distance cutoff in Angstrom
1572 @param metric (Default=IMP.atom.get_rmsd) the metric that will be used to compute rmsds
1574 self.clean_clusters()
1575 not_clustered=set(range(len(self.stath1)))
1576 while len(not_clustered)>0:
1577 self.
aggregate(not_clustered, rmsd_cutoff, metric)
1583 Refine the clusters by merging the ones whose centers are close
1584 @param rmsd_cutoff cutoff distance in Angstorms
1587 clusters_copy=self.clusters
1588 for c0,c1
in itertools.combinations(self.clusters,2):
1589 if c0.center_index
is None:
1591 if c1.center_index
is None:
1593 d0=self.stath0[c0.center_index]
1594 d1=self.stath1[c1.center_index]
1595 rmsd, molecular_assignment = self.
rmsd()
1596 if rmsd <= rmsd_cutoff:
1597 if c1
in self.clusters:
1598 clusters_copy.remove(c1)
1600 self.clusters=clusters_copy
1607 def set_cluster_assignments(self, cluster_ids):
1608 if len(cluster_ids)!=len(self.stath0):
1609 raise ValueError(
'cluster ids has to be same length as number of frames')
1612 for i
in sorted(list(set(cluster_ids))):
1614 for i, (idx, d)
in enumerate(zip(cluster_ids, self.stath0)):
1615 self.clusters[idx].add_member(i,d)
1619 Return the model data from a cluster
1620 @param cluster IMP.pmi.output.Cluster object
1629 Save the data for the whole models into a pickle file
1630 @param filename string
1632 self.stath0.save_data(filename)
1636 Set the data from an external IMP.pmi.output.Data
1637 @param data IMP.pmi.output.Data
1639 self.stath0.data=data
1640 self.stath1.data=data
1644 Load the data from an external pickled file
1645 @param filename string
1647 self.stath0.load_data(filename)
1648 self.stath1.load_data(filename)
1649 self.best_models=len(self.stath0)
1651 def add_cluster(self,rmf_name_list):
1653 print(
"creating cluster index "+str(len(self.clusters)))
1654 self.clusters.append(c)
1655 current_len=len(self.stath0)
1657 for rmf
in rmf_name_list:
1658 print(
"adding rmf "+rmf)
1659 self.stath0.add_stat_file(rmf)
1660 self.stath1.add_stat_file(rmf)
1662 for n0
in range(current_len,len(self.stath0)):
1669 Save the clusters into a pickle file
1670 @param filename string
1673 import cPickle
as pickle
1676 fl=open(filename,
'wb')
1677 pickle.dump(self.clusters,fl)
1681 Load the clusters from a pickle file
1682 @param filename string
1683 @param append bool (Default=False), if True. append the clusters to the ones currently present
1686 import cPickle
as pickle
1689 fl=open(filename,
'rb')
1690 self.clean_clusters()
1692 self.clusters+=pickle.load(fl)
1694 self.clusters=pickle.load(fl)
1703 Compute the cluster center for a given cluster
1705 member_distance=defaultdict(float)
1707 for n0,n1
in itertools.combinations(cluster.members,2):
1710 rmsd, _ = self.
rmsd()
1711 member_distance[n0]+=rmsd
1713 if len(member_distance)>0:
1714 cluster.center_index=min(member_distance, key=member_distance.get)
1716 cluster.center_index=cluster.members[0]
1720 Save the coordinates of the current cluster a single rmf file
1722 print(
"saving coordinates",cluster)
1725 if rmf_name
is None:
1726 rmf_name=prefix+
'/'+str(cluster.cluster_id)+
".rmf3"
1728 d1=self.stath1[cluster.members[0]]
1730 o.init_rmf(rmf_name, [self.stath1])
1731 for n1
in cluster.members:
1735 if self.alignment: self.align()
1736 o.write_rmf(rmf_name)
1738 o.close_rmf(rmf_name)
1742 remove structures that are similar
1743 append it to a new cluster
1745 print(
"pruning models")
1748 remaining=range(1,len(self.stath1),10)
1750 while len(remaining)>0:
1751 d0=self.stath0[selected]
1753 for n1
in remaining:
1755 if self.alignment: self.align()
1759 print(
"pruning model %s, similar to model %s, rmsd %s"%(str(n1),str(selected),str(d)))
1760 remaining=[x
for x
in remaining
if x
not in rm]
1761 if len(remaining)==0:
break
1762 selected=remaining[0]
1763 filtered.append(selected)
1766 self.clusters.append(c)
1776 Compute the precision of a cluster
1782 if not cluster.center_index
is None:
1783 members1=[cluster.center_index]
1785 members1=cluster.members
1789 for n1
in cluster.members:
1794 tmp_rmsd, _ = self.
rmsd()
1799 precision=rmsd/npairs
1800 cluster.precision=precision
1805 Compute the bipartite precision (ie the cross-precision)
1806 between two clusters
1810 for cn0,n0
in enumerate(cluster1.members):
1812 for cn1,n1
in enumerate(cluster2.members):
1814 tmp_rmsd, _ =self.
rmsd()
1815 if verbose: print(
"--- rmsd between structure %s and structure %s is %s"%(str(cn0),str(cn1),str(tmp_rmsd)))
1818 precision=rmsd/npairs
1821 def rmsf(self,cluster,molecule,copy_index=0,state_index=0,cluster_ref=None,step=1):
1823 Compute the Root mean square fluctuations
1824 of a molecule in a cluster
1825 Returns an IMP.pmi.tools.OrderedDict() where the keys are the residue indexes and the value is the rmsf
1827 rmsf=IMP.pmi.tools.OrderedDict()
1830 if cluster_ref
is not None:
1831 if not cluster_ref.center_index
is None:
1832 members0 = [cluster_ref.center_index]
1834 members0 = cluster_ref.members
1836 if not cluster.center_index
is None:
1837 members0 = [cluster.center_index]
1839 members0 = cluster.members
1842 copy_index=copy_index,state_index=state_index)
1843 ps0=s0.get_selected_particles()
1857 for n1
in cluster.members[::step]:
1859 print(
"--- rmsf %s %s"%(str(n0),str(n1)))
1862 s1=
IMP.atom.Selection(self.stath1,molecule=molecule,residue_indexes=residue_indexes,resolution=1,
1863 copy_index=copy_index,state_index=state_index)
1864 ps1 = s1.get_selected_particles()
1867 if self.alignment: self.align()
1868 for n,(p0,p1)
in enumerate(zip(ps0,ps1)):
1869 r=residue_indexes[n]
1881 for stath
in [self.stath0,self.stath1]:
1882 if molecule
not in self.symmetric_molecules:
1884 copy_index=copy_index,state_index=state_index)
1887 state_index=state_index)
1889 ps = s.get_selected_particles()
1898 def save_densities(self,cluster,density_custom_ranges,voxel_size=5,reference="Absolute", prefix="./",step=1):
1903 for n1
in cluster.members[::step]:
1904 print(
"density "+str(n1))
1907 if self.alignment: self.align()
1908 dens.add_subunits_density(self.stath1)
1910 dens.write_mrc(path=prefix+
'/',suffix=str(cluster.cluster_id))
1913 def contact_map(self,cluster,contact_threshold=15,log_scale=False,consolidate=False,molecules=None,prefix='./',reference="Absolute"):
1916 import matplotlib.pyplot
as plt
1917 import matplotlib.cm
as cm
1918 from scipy.spatial.distance
import cdist
1920 if molecules
is None:
1924 unique_copies=[mol
for mol
in mols
if mol.get_copy_index() == 0]
1925 mol_names_unique=dict((mol.get_name(),mol)
for mol
in unique_copies)
1926 total_len_unique=sum(max(mol.get_residue_indexes())
for mol
in unique_copies)
1935 seqlen=max(mol.get_residue_indexes())
1936 index_dict[mol] = range(prev_stop, prev_stop + seqlen)
1940 for mol
in unique_copies:
1941 seqlen=max(mol.get_residue_indexes())
1942 index_dict[mol] = range(prev_stop, prev_stop + seqlen)
1946 for ncl,n1
in enumerate(cluster.members):
1950 coord_dict=IMP.pmi.tools.OrderedDict()
1952 mol_name=mol.get_name()
1953 copy_index=mol.get_copy_index()
1954 rindexes = mol.get_residue_indexes()
1955 coords = np.ones((max(rindexes),3))
1956 for rnum
in rindexes:
1958 selpart = sel.get_selected_particles()
1959 if len(selpart) == 0:
continue
1960 selpart = selpart[0]
1961 coords[rnum - 1, :] =
IMP.core.XYZ(selpart).get_coordinates()
1962 coord_dict[mol]=coords
1965 coords=np.concatenate(list(coord_dict.values()))
1966 dists = cdist(coords, coords)
1967 binary_dists = np.where((dists <= contact_threshold) & (dists >= 1.0), 1.0, 0.0)
1969 binary_dists_dict={}
1971 len1 = max(mol1.get_residue_indexes())
1973 name1=mol1.get_name()
1974 name2=mol2.get_name()
1975 dists = cdist(coord_dict[mol1], coord_dict[mol2])
1976 if (name1, name2)
not in binary_dists_dict:
1977 binary_dists_dict[(name1, name2)] = np.zeros((len1,len1))
1978 binary_dists_dict[(name1, name2)] += np.where((dists <= contact_threshold) & (dists >= 1.0), 1.0, 0.0)
1979 binary_dists=np.zeros((total_len_unique,total_len_unique))
1981 for name1,name2
in binary_dists_dict:
1982 r1=index_dict[mol_names_unique[name1]]
1983 r2=index_dict[mol_names_unique[name2]]
1984 binary_dists[min(r1):max(r1)+1,min(r2):max(r2)+1] = np.where((binary_dists_dict[(name1, name2)]>=1.0),1.0,0.0)
1989 contact_freqs = binary_dists
1991 dist_maps.append(dists)
1992 av_dist_map += dists
1993 contact_freqs += binary_dists
1998 contact_freqs =-np.log(1.0-1.0/(len(cluster)+1)*contact_freqs)
2000 contact_freqs =1.0/len(cluster)*contact_freqs
2001 av_dist_map=1.0/len(cluster)*contact_freqs
2003 fig = plt.figure(figsize=(100, 100))
2004 ax = fig.add_subplot(111)
2007 gap_between_components=50
2019 prot_list=list(zip(*sorted_tuple))[1]
2022 prot_list=list(zip(*sorted_tuple))[1]
2024 prot_listx = prot_list
2025 nresx = gap_between_components + \
2026 sum([max(mol.get_residue_indexes())
2027 + gap_between_components
for mol
in prot_listx])
2030 prot_listy = prot_list
2031 nresy = gap_between_components + \
2032 sum([max(mol.get_residue_indexes())
2033 + gap_between_components
for mol
in prot_listy])
2038 res = gap_between_components
2039 for mol
in prot_listx:
2040 resoffsetx[mol] = res
2041 res += max(mol.get_residue_indexes())
2043 res += gap_between_components
2047 res = gap_between_components
2048 for mol
in prot_listy:
2049 resoffsety[mol] = res
2050 res += max(mol.get_residue_indexes())
2052 res += gap_between_components
2054 resoffsetdiagonal = {}
2055 res = gap_between_components
2056 for mol
in IMP.pmi.tools.OrderedSet(prot_listx + prot_listy):
2057 resoffsetdiagonal[mol] = res
2058 res += max(mol.get_residue_indexes())
2059 res += gap_between_components
2064 for n, prot
in enumerate(prot_listx):
2065 res = resoffsetx[prot]
2067 for proty
in prot_listy:
2068 resy = resoffsety[proty]
2069 endy = resendy[proty]
2070 ax.plot([res, res], [resy, endy], linestyle=
'-',color=
'gray', lw=0.4)
2071 ax.plot([end, end], [resy, endy], linestyle=
'-',color=
'gray', lw=0.4)
2072 xticks.append((float(res) + float(end)) / 2)
2077 for n, prot
in enumerate(prot_listy):
2078 res = resoffsety[prot]
2080 for protx
in prot_listx:
2081 resx = resoffsetx[protx]
2082 endx = resendx[protx]
2083 ax.plot([resx, endx], [res, res], linestyle=
'-',color=
'gray', lw=0.4)
2084 ax.plot([resx, endx], [end, end], linestyle=
'-',color=
'gray', lw=0.4)
2085 yticks.append((float(res) + float(end)) / 2)
2090 tmp_array = np.zeros((nresx, nresy))
2092 for px
in prot_listx:
2093 for py
in prot_listy:
2094 resx = resoffsetx[px]
2095 lengx = resendx[px] - 1
2096 resy = resoffsety[py]
2097 lengy = resendy[py] - 1
2098 indexes_x = index_dict[px]
2099 minx = min(indexes_x)
2100 maxx = max(indexes_x)
2101 indexes_y = index_dict[py]
2102 miny = min(indexes_y)
2103 maxy = max(indexes_y)
2104 tmp_array[resx:lengx,resy:lengy] = contact_freqs[minx:maxx,miny:maxy]
2105 ret[(px,py)]=np.argwhere(contact_freqs[minx:maxx,miny:maxy] == 1.0)+1
2107 cax = ax.imshow(tmp_array,
2112 interpolation=
'nearest')
2114 ax.set_xticks(xticks)
2115 ax.set_xticklabels(xlabels, rotation=90)
2116 ax.set_yticks(yticks)
2117 ax.set_yticklabels(ylabels)
2118 plt.setp(ax.get_xticklabels(), fontsize=6)
2119 plt.setp(ax.get_yticklabels(), fontsize=6)
2122 fig.set_size_inches(0.005 * nresx, 0.005 * nresy)
2123 [i.set_linewidth(2.0)
for i
in ax.spines.values()]
2128 plt.savefig(prefix+
"/contact_map."+str(cluster.cluster_id)+
".pdf", dpi=300,transparent=
"False")
2132 def plot_rmsd_matrix(self,filename):
2134 self.compute_all_pairwise_rmsd()
2135 distance_matrix = np.zeros(
2136 (len(self.stath0), len(self.stath1)))
2137 for (n0,n1)
in self.pairwise_rmsd:
2138 distance_matrix[n0, n1] = self.pairwise_rmsd[(n0,n1)]
2140 import matplotlib
as mpl
2142 import matplotlib.pylab
as pl
2143 from scipy.cluster
import hierarchy
as hrc
2145 fig = pl.figure(figsize=(10,8))
2146 ax = fig.add_subplot(212)
2147 dendrogram = hrc.dendrogram(
2148 hrc.linkage(distance_matrix),
2151 leaves_order = dendrogram[
'leaves']
2152 ax.set_xlabel(
'Model')
2153 ax.set_ylabel(
'RMSD [Angstroms]')
2155 ax2 = fig.add_subplot(221)
2157 distance_matrix[leaves_order,
2160 interpolation=
'nearest')
2161 cb = fig.colorbar(cax)
2162 cb.set_label(
'RMSD [Angstroms]')
2163 ax2.set_xlabel(
'Model')
2164 ax2.set_ylabel(
'Model')
2166 pl.savefig(filename, dpi=300)
2175 Update the cluster id numbers
2177 for n,c
in enumerate(self.clusters):
2180 def get_molecule(self, hier, name, copy):
2188 self.seldict0=IMP.pmi.tools.get_selections_dictionary(self.sel0_rmsd.get_selected_particles())
2189 self.seldict1=IMP.pmi.tools.get_selections_dictionary(self.sel1_rmsd.get_selected_particles())
2190 for mol
in self.seldict0:
2191 for sel
in self.seldict0[mol]:
2192 self.issymmetricsel[sel]=
False
2193 for mol
in self.symmetric_molecules:
2194 self.symmetric_molecules[mol]=len(self.seldict0[mol])
2195 for sel
in self.seldict0[mol]:
2196 self.issymmetricsel[sel]=
True
2203 for rb
in self.rbs1:
2206 for bead
in self.beads1:
2211 def aggregate(self, idxs, rmsd_cutoff=10, metric=IMP.atom.get_rmsd):
2213 initial filling of the clusters.
2216 print(
"clustering model "+str(n0))
2217 d0 = self.stath0[n0]
2219 print(
"creating cluster index "+str(len(self.clusters)))
2220 self.clusters.append(c)
2222 clustered = set([n0])
2224 print(
"--- trying to add model "+str(n1)+
" to cluster "+str(len(self.clusters)))
2225 d1 = self.stath1[n1]
2226 if self.alignment: self.align()
2227 rmsd, _ = self.
rmsd(metric=metric)
2228 if rmsd<rmsd_cutoff:
2229 print(
"--- model "+str(n1)+
" added, rmsd="+str(rmsd))
2233 print(
"--- model "+str(n1)+
" NOT added, rmsd="+str(rmsd))
2238 merge the clusters that have close members
2239 @param rmsd_cutoff cutoff distance in Angstorms
2245 for c0, c1
in filter(
lambda x: len(x[0].members)>1, itertools.combinations(self.clusters, 2)):
2246 n0, n1 = [c.members[0]
for c
in (c0,c1)]
2247 d0 = self.stath0[n0]
2248 d1 = self.stath1[n1]
2249 rmsd, _ = self.
rmsd()
2251 to_merge.append((c0,c1))
2253 for c0, c
in reversed(to_merge):
2257 self.clusters = [c
for c
in filter(
lambda x: len(x.members)>0, self.clusters)]
2262 returns true if c0 and c1 have members that are closer than rmsd_cutoff
2264 print(
"check close members for clusters "+str(c0.cluster_id)+
" and "+str(c1.cluster_id))
2265 for n0, n1
in itertools.product(c0.members[1:], c1.members):
2266 d0 = self.stath0[n0]
2267 d1 = self.stath1[n1]
2268 rmsd, _ = self.
rmsd(metric=metric)
2269 if rmsd<rmsd_cutoff:
2285 a function that returns the permutation best_sel of sels0 that minimizes metric
2287 best_rmsd2 = float(
'inf')
2289 if self.issymmetricsel[sels0[0]]:
2292 for offset
in range(N):
2293 order=[(offset+i)%N
for i
in range(N)]
2294 sels = [sels0[(offset+i)%N]
for i
in range(N)]
2297 r=metric(sel0, sel1)
2300 if rmsd2 < best_rmsd2:
2305 for sels
in itertools.permutations(sels0):
2307 for sel0, sel1
in itertools.takewhile(
lambda x: rmsd2<best_rmsd2, zip(sels, sels1)):
2308 r=metric(sel0, sel1)
2310 if rmsd2 < best_rmsd2:
2324 return best_sel, best_rmsd2
2326 def compute_all_pairwise_rmsd(self):
2327 for d0
in self.stath0:
2328 for d1
in self.stath1:
2329 rmsd, _ = self.
rmsd()
2331 def rmsd(self,metric=IMP.atom.get_rmsd):
2333 Computes the RMSD. Resolves ambiguous pairs assignments
2336 n0=self.stath0.current_index
2337 n1=self.stath1.current_index
2338 if ((n0,n1)
in self.pairwise_rmsd)
and ((n0,n1)
in self.pairwise_molecular_assignment):
2339 return self.pairwise_rmsd[(n0,n1)], self.pairwise_molecular_assignment[(n0,n1)]
2347 seldict_best_order={}
2348 molecular_assignment={}
2349 for molname, sels0
in self.seldict0.items():
2350 sels_best_order, best_rmsd2 = self.
rmsd_helper(sels0, self.seldict1[molname], metric)
2352 Ncoords = len(sels_best_order[0].get_selected_particles())
2353 Ncopies = len(self.seldict1[molname])
2354 total_rmsd += Ncoords*best_rmsd2
2355 total_N += Ncoords*Ncopies
2357 for sel0, sel1
in zip(sels_best_order, self.seldict1[molname]):
2358 p0 = sel0.get_selected_particles()[0]
2359 p1 = sel1.get_selected_particles()[0]
2365 molecular_assignment[(molname,c0)]=(molname,c1)
2367 total_rmsd = math.sqrt(total_rmsd/total_N)
2369 self.pairwise_rmsd[(n0,n1)]=total_rmsd
2370 self.pairwise_molecular_assignment[(n0,n1)]=molecular_assignment
2371 self.pairwise_rmsd[(n1,n0)]=total_rmsd
2372 self.pairwise_molecular_assignment[(n1,n0)]=molecular_assignment
2374 return total_rmsd, molecular_assignment
2378 Fix the reference structure for structural alignment, rmsd and chain assignemnt
2379 @param reference can be either "Absolute" (cluster center of the first cluster)
2380 or Relative (cluster center of the current cluster)
2382 if reference==
"Absolute":
2384 elif reference==
"Relative":
2385 if cluster.center_index:
2386 n0=cluster.center_index
2388 n0=cluster.members[0]
2393 compute the molecular assignemnts between multiple copies
2394 of the same sequence. It changes the Copy index of Molecules
2397 _, molecular_assignment = self.
rmsd()
2398 for (m0, c0), (m1,c1)
in molecular_assignment.items():
2399 mol0 = self.molcopydict0[m0][c0]
2400 mol1 = self.molcopydict1[m1][c1]
2403 p1.set_value(cik0,c0)
2407 Undo the Copy index assignment
2410 _, molecular_assignment = self.
rmsd()
2412 for (m0, c0), (m1,c1)
in molecular_assignment.items():
2413 mol0 = self.molcopydict0[m0][c0]
2414 mol1 = self.molcopydict1[m1][c1]
2417 p1.set_value(cik0,c1)
2424 s=
"AnalysisReplicaExchange\n"
2425 s+=
"---- number of clusters %s \n"%str(len(self.clusters))
2426 s+=
"---- number of models %s \n"%str(len(self.stath0))
2429 def __getitem__(self,int_slice_adaptor):
2430 if isinstance(int_slice_adaptor, int):
2431 return self.clusters[int_slice_adaptor]
2432 elif isinstance(int_slice_adaptor, slice):
2433 return self.__iter__(int_slice_adaptor)
2435 raise TypeError(
"Unknown Type")
2438 return len(self.clusters)
2440 def __iter__(self,slice_key=None):
2441 if slice_key
is None:
2442 for i
in range(len(self)):
2445 for i
in range(len(self))[slice_key]:
Simplify creation of constraints and movers for an IMP Hierarchy.
def rmsd
Computes the RMSD.
def set_reference
Fix the reference structure for structural alignment, rmsd and chain assignemnt.
def load_clusters
Load the clusters from a pickle file.
def precision
Compute the precision of a cluster.
Extends the functionality of IMP.atom.Molecule.
A macro for running all the basic operations of analysis.
A container for models organized into clusters.
Sample using molecular dynamics.
def aggregate
initial filling of the clusters.
A class for reading stat files (either rmf or ascii v1 and v2)
A member of a rigid body, it has internal (local) coordinates.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Set of Python classes to create a multi-state, multi-resolution IMP hierarchy.
def prune_redundant_structures
remove structures that are similar append it to a new cluster
def rmsf
Compute the Root mean square fluctuations of a molecule in a cluster Returns an IMP.pmi.tools.OrderedDict() where the keys are the residue indexes and the value is the rmsf.
static XYZR setup_particle(Model *m, ParticleIndex pi)
Utility classes and functions for reading and storing PMI files.
def get_best_models
Given a list of stat files, read them all and find the best models.
A helper output for model evaluation.
def set_rmsd_selection
Setup the selection onto which the rmsd is computed.
def get_cluster_data
Return the model data from a cluster.
def __init__
Construction of the Class.
def get_molecules
Return list of all molecules grouped by state.
def set_data
Set the data from an external IMP.pmi.output.Data.
def undo_apply_molecular_assignments
Undo the Copy index assignment.
def set_alignment_selection
Setup the selection onto which the alignment is computed.
def rmsd_helper
a function that returns the permutation best_sel of sels0 that minimizes metric
def save_coordinates
Save the coordinates of the current cluster a single rmf file.
def clustering
Get the best scoring models, compute a distance matrix, cluster them, and create density maps...
def apply_molecular_assignments
compute the molecular assignemnts between multiple copies of the same sequence.
This class contains analysis utilities to investigate ReplicaExchange results.
Add uncertainty to a particle.
A macro to build a IMP::pmi::topology::System based on a TopologyReader object.
def merge_aggregates
merge the clusters that have close members
This class initializes the root node of the global IMP.atom.Hierarchy.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
A class to cluster structures.
def add_protocol_output
Capture details of the modeling protocol.
static Uncertainty setup_particle(Model *m, ParticleIndex pi, Float uncertainty)
Representation of the system.
def compute_cluster_center
Compute the cluster center for a given cluster.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def get_modeling_trajectory
Get a trajectory of the modeling run, for generating demonstrative movies.
Warning related to handling of structures.
A decorator for keeping track of copies of a molecule.
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
def get_trajectory_models
Given a list of stat files, read them all and find a trajectory of models.
The standard decorator for manipulating molecular structures.
Performs alignment and RMSD calculation for two sets of coordinates.
def update_seldicts
Update the seldicts.
def update_clusters
Update the cluster id numbers.
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
def refine
Refine the clusters by merging the ones whose centers are close.
A decorator for a particle with x,y,z coordinates.
Class for easy writing of PDBs, RMFs, and stat files.
def set_symmetric
Store names of symmetric molecules.
Warning for an expected, but missing, file.
Tools for clustering and cluster analysis.
bool get_is_canonical(atom::Hierarchy h)
Walk up a PMI2 hierarchy/representations and check if the root is named System.
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Classes for writing output files and processing them.
def deprecated_object
Python decorator to mark a class as deprecated.
Sample using Monte Carlo.
Create movers and set up constraints for PMI objects.
def merge
merge two clusters
def add_state
Add a state using the topology info in a IMP::pmi::topology::TopologyReader object.
The general base class for IMP exceptions.
static SampleProvenance setup_particle(Model *m, ParticleIndex pi, std::string method, int frames, int iterations, int replicas)
class to link stat files to several rmf files
def save_data
Save the data for the whole models into a pickle file.
Class to handle individual particles of a Model object.
def execute_macro
Builds representations and sets up degrees of freedom.
def bipartite_precision
Compute the bipartite precision (ie the cross-precision) between two clusters.
def read_coordinates_of_rmfs
Read in coordinates of a set of RMF tuples.
int get_copy_index(Hierarchy h)
Walk up the hierarchy to find the current copy index.
def cluster
Cluster the models based on RMSD.
static bool get_is_setup(Model *m, ParticleIndex pi)
def save_clusters
Save the clusters into a pickle file.
def have_close_members
returns true if c0 and c1 have members that are closer than rmsd_cutoff
void add_geometries(RMF::FileHandle file, const display::GeometriesTemp &r)
Add geometries to the file.
A macro to help setup and run replica exchange.
algebra::Transformation3D get_transformation_aligning_first_to_second(const Selection &s1, const Selection &s2)
Get the transformation to align two selections.
A dictionary-like wrapper for reading and storing sequence data.
void add_provenance(Model *m, ParticleIndex pi, Provenance p)
Add provenance to part of the model.
Hierarchies get_leaves(const Selection &h)
Select hierarchy particles identified by the biological name.
Compute mean density maps from structures.
def load_data
Load the data from an external pickled file.
Support for the RMF file format for storing hierarchical molecular data and markup.
Sample using replica exchange.
Warning for probably incorrect input parameters.
def add_provenance
Add provenance information in prov (a list of _TempProvenance objects) to each of the IMP hierarchies...
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A decorator for a particle with x,y,z coordinates and a radius.