3 """@namespace IMP.pmi.representation
4 Representation of the system.
7 from __future__
import print_function
20 from math
import pi, sqrt
21 from operator
import itemgetter
26 def __init__(self, doi, root):
30 def get_fname(self, fname):
31 """Return a path relative to the top of the repository"""
32 return os.path.relpath(fname, self._root)
34 class _StateInfo(object):
35 """Score state-specific information about this representation."""
40 "use IMP.pmi.topology.System instead")
41 class Representation(object):
45 Set up the representation of all proteins and nucleic acid macromolecules.
47 Create the molecular hierarchies, representation,
48 sequence connectivity for the various involved proteins and
49 nucleic acid macromolecules:
51 Create a protein, DNA or RNA, represent it as a set of connected balls of appropriate
52 radii and number of residues, PDB at given resolution(s), or ideal helices.
54 How to use the SimplifiedModel class (typical use):
56 see test/test_hierarchy_contruction.py
60 1) Create a chain of helices and flexible parts
62 c_1_119 =self.add_component_necklace("prot1",1,119,20)
63 c_120_131 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(120,131))
64 c_132_138 =self.add_component_beads("prot1",[(132,138)])
65 c_139_156 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(139,156))
66 c_157_174 =self.add_component_beads("prot1",[(157,174)])
67 c_175_182 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(175,182))
68 c_183_194 =self.add_component_beads("prot1",[(183,194)])
69 c_195_216 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(195,216))
70 c_217_250 =self.add_component_beads("prot1",[(217,250)])
73 self.set_rigid_body_from_hierarchies(c_120_131)
74 self.set_rigid_body_from_hierarchies(c_139_156)
75 self.set_rigid_body_from_hierarchies(c_175_182)
76 self.set_rigid_body_from_hierarchies(c_195_216)
78 clist=[c_1_119,c_120_131,c_132_138,c_139_156,c_157_174,c_175_182,c_183_194,c_195_216,
81 self.set_chain_of_super_rigid_bodies(clist,2,3)
83 self.set_super_rigid_bodies(["prot1"])
87 def __init__(self, model, upperharmonic=True, disorderedlength=True):
89 @param model the model
90 @param upperharmonic This flag uses either harmonic (False)
91 or upperharmonic (True) in the intra-pair
92 connectivity restraint.
93 @param disorderedlength This flag uses either disordered length
94 calculated for random coil peptides (True) or zero
95 surface-to-surface distance between beads (False)
96 as optimal distance for the sequence connectivity
100 self.state = _StateInfo()
102 self._file_dataset = {}
103 self._protocol_output = []
104 self._non_modeled_components = {}
113 self.upperharmonic = upperharmonic
114 self.disorderedlength = disorderedlength
115 self.rigid_bodies = []
116 self.fixed_rigid_bodies = []
117 self.fixed_floppy_bodies = []
118 self.floppy_bodies = []
124 self.super_rigid_bodies = []
125 self.rigid_body_symmetries = []
126 self.output_level =
"low"
129 self.maxtrans_rb = 2.0
130 self.maxrot_rb = 0.04
131 self.maxtrans_srb = 2.0
132 self.maxrot_srb = 0.2
133 self.rigidbodiesarefixed =
False
134 self.floppybodiesarefixed =
False
135 self.maxtrans_fb = 3.0
136 self.resolution = 10.0
137 self.bblenght = 100.0
141 self.representation_is_modified =
False
142 self.unmodeledregions_cr_dict = {}
143 self.sortedsegments_cr_dict = {}
145 self.connected_intra_pairs = []
148 self.sequence_dict = {}
149 self.hier_geometry_pairs = {}
155 self.hier_representation = {}
156 self.hier_resolution = {}
159 self.reference_structures = {}
162 self.linker_restraints.set_was_used(
True)
163 self.linker_restraints_dict = {}
164 self.threetoone = {
'ALA':
'A',
'ARG':
'R', 'ASN': 'N', 'ASP': 'D',
165 'CYS':
'C',
'GLU':
'E',
'GLN':
'Q',
'GLY':
'G',
166 'HIS':
'H',
'ILE':
'I',
'LEU':
'L',
'LYS':
'K',
167 'MET':
'M',
'PHE':
'F',
'PRO':
'P',
'SER':
'S',
168 'THR':
'T',
'TRP':
'W',
'TYR':
'Y',
'VAL':
'V',
'UNK':
'X'}
170 self.onetothree = dict((v, k)
for k, v
in self.threetoone.items())
180 """Associate some metadata with this modeling.
181 @param m an instance of an ihm metadata class, such as
182 ihm.Software, ihm.Citation, or ihm.location.Repository.
184 self._metadata.append(m)
187 """Associate a dataset with a filename.
188 This can be used to identify how the file was produced (in many
189 cases IMP can determine this automatically from a file header or
190 other metadata, but not always). For example, a manually-produced
191 PDB file (not from the PDB database or Modeller) can be
193 @param fname filename
194 @dataset the ihm.dataset.Dataset object to associate.
196 self._file_dataset[os.path.abspath(fname)] = dataset
199 """Get the dataset associated with a filename, or None.
200 @param fname filename
201 @return an ihm.dataset.Dataset, or None.
203 return self._file_dataset.get(os.path.abspath(fname),
None)
206 """Capture details of the modeling protocol.
207 @param p an instance of IMP.pmi.output.ProtocolOutput or a subclass.
209 state = p._add_state(self)
210 self._protocol_output.append((p, state))
211 p._each_metadata.append(self._metadata)
212 p._file_datasets.append(self._file_dataset)
213 state.model = self.model
214 state.prot = self.prot
215 protocol_output = property(
lambda self:
216 [x[0]
for x
in self._protocol_output])
218 def set_label(self, label):
221 def create_component(self, name, color=0.0):
228 protein_h.set_name(name)
229 self.hier_dict[name] = protein_h
230 self.hier_representation[name] = {}
231 self.hier_db.add_name(name)
232 self.prot.add_child(protein_h)
233 self.color_dict[name] = color
234 self.elements[name] = []
235 for p, state
in self._protocol_output:
236 p.create_component(state, name,
True)
239 """Create a transformed view of a component.
240 This does not copy the coordinates or representation of the
241 component, and the transformed component is not visible to IMP,
242 but when the model is written to a ProtocolOutput, the transformed
243 component is added. This can be used to easily add symmetry-related
244 'copies' of a component without copying the underlying
246 for p, state
in self._protocol_output:
247 p.create_transformed_component(state, name, original, transform)
250 """Create a component that isn't used in the modeling.
251 No coordinates or other structural data for this component will
252 be read or written, but a primary sequence can be assigned. This
253 is useful if the input experimental data is of a system larger
254 than that modeled. Any references to these non-modeled components
255 can then be correctly resolved later."""
256 self._non_modeled_components[name] =
None
257 self.elements[name] = []
258 for p, state
in self._protocol_output:
259 p.create_component(state, name,
False)
261 def get_component_names(self):
262 return list(self.hier_dict.keys())
267 Add the primary sequence for a single component.
269 @param name Human-readable name of the component
270 @param filename Name of the FASTA file
271 @param id Identifier of the sequence in the FASTA file header
272 (if not provided, use `name` instead)
277 if id
not in record_dict:
278 raise KeyError(
"id %s not found in fasta file" % id)
279 length = len(record_dict[id])
280 self.sequence_dict[name] = str(record_dict[id])
282 if name
not in self._non_modeled_components:
283 protein_h = self.hier_dict[name]
284 protein_h.set_sequence(self.sequence_dict[name])
287 self.sequence_dict[name]=offs_str+self.sequence_dict[name]
289 self.elements[name].append((length, length,
" ",
"end"))
290 for p, state
in self._protocol_output:
291 p.add_component_sequence(state, name, self.sequence_dict[name])
293 def autobuild_model(self, name, pdbname, chain,
294 resolutions=
None, resrange=
None,
296 color=
None, pdbresrange=
None, offset=0,
297 show=
False, isnucleicacid=
False,
300 self.representation_is_modified =
True
304 color = self.color_dict[name]
306 self.color_dict[name] = color
308 if resolutions
is None:
310 print(
"autobuild_model: constructing %s from pdb %s and chain %s" % (name, pdbname, str(chain)))
319 t.get_children()[0].get_children()[0]).
get_index()
321 t.get_children()[0].get_children()[-1]).
get_index()
327 if name
in self.sequence_dict:
328 resrange = (1, len(self.sequence_dict[name]))
330 resrange = (start + offset, end + offset)
332 if resrange[1]
in (-1,
'END'):
333 resrange = (resrange[0],end)
334 start = resrange[0] - offset
335 end = resrange[1] - offset
354 for n, g
in enumerate(gaps):
358 print(
"autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
360 chain, resolutions=resolutions,
361 color=color, cacenters=
True,
362 resrange=(first, last),
363 offset=offset, isnucleicacid=isnucleicacid)
364 elif g[2] ==
"gap" and n > 0:
365 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
366 parts = self.hier_db.get_particles_at_closest_resolution(name,
371 first+offset, last+offset, missingbeadsize, incoord=xyz)
373 elif g[2] ==
"gap" and n == 0:
375 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
377 first+offset, last+offset, missingbeadsize, incoord=xyznter)
382 resrange=
None, offset=0, cacenters=
True, show=
False,
383 isnucleicacid=
False, readnonwateratoms=
False,
384 read_ca_cb_only=
False):
386 Add a component that has an associated 3D structure in a PDB file.
388 Reads the PDB, and constructs the fragments corresponding to contiguous
391 @return a list of hierarchies.
393 @param name (string) the name of the component
394 @param pdbname (string) the name of the PDB file
395 @param chain (string or integer) can be either a string (eg, "A")
396 or an integer (eg, 0 or 1) in case you want
397 to get the corresponding chain number in the PDB.
398 @param resolutions (integers) a list of integers that corresponds
399 to the resolutions that have to be generated
400 @param color (float from 0 to 1) the color applied to the
401 hierarchies generated
402 @param resrange (tuple of integers): the residue range to extract
403 from the PDB. It is a tuple (beg,end). If not specified,
404 all residues belonging to the specified chain are read.
405 @param offset (integer) specifies the residue index offset to be
406 applied when reading the PDB (the FASTA sequence is
407 assumed to start from residue 1, so use this parameter
408 if the PDB file does not start at residue 1)
409 @param cacenters (boolean) if True generates resolution=1 beads
410 centered on C-alpha atoms.
411 @param show (boolean) print out the molecular hierarchy at the end.
412 @param isnucleicacid (boolean) use True if you're reading a PDB
414 @param readnonwateratoms (boolean) if True fixes some pathological PDB.
415 @param read_ca_cb_only (boolean) if True, only reads CA/CB
418 self.representation_is_modified =
True
421 color = self.color_dict[name]
422 protein_h = self.hier_dict[name]
432 if type(chain) == str:
440 t.get_children()[0].get_children()[0]).
get_index()
442 t.get_children()[0].get_children()[-1]).
get_index()
443 c =
IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
447 IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
454 if not resrange
is None:
455 if resrange[0] > start:
457 if resrange[1] < end:
460 if not isnucleicacid:
464 residue_indexes=list(range(
467 atom_type=IMP.atom.AT_CA)
474 residue_indexes=list(range(
477 atom_type=IMP.atom.AT_P)
479 ps = sel.get_selected_particles()
481 raise ValueError(
"%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
482 % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
490 if par.get_parent() != c0:
491 par.get_parent().remove_child(par)
493 start = start + offset
496 self.elements[name].append(
497 (start, end, pdbname.split(
"/")[-1] +
":" + chain,
"pdb"))
500 resolutions, isnucleicacid, c0, protein_h,
"pdb", color)
501 self._copy_pdb_provenance(t, hiers[0], offset)
503 for p, state
in self._protocol_output:
504 p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
518 for r
in residues.keys():
520 self.model.remove_particle(c0)
526 def _copy_pdb_provenance(self, pdb, h, offset):
527 """Copy the provenance information from the PDB to our hierarchy"""
528 c =
IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
535 def add_component_ideal_helix(
543 self.representation_is_modified =
True
544 from math
import pi, cos, sin
546 protein_h = self.hier_dict[name]
549 color = self.color_dict[name]
553 self.elements[name].append((start, end,
" ",
"helix"))
555 for n, res
in enumerate(range(start, end + 1)):
556 if name
in self.sequence_dict:
558 rtstr = self.onetothree[
559 self.sequence_dict[name][res-1]]
579 x = 2.3 * cos(n * 2 * pi / 3.6)
580 y = 2.3 * sin(n * 2 * pi / 3.6)
581 z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
590 resolutions,
False, c0, protein_h,
"helix", color)
599 """ add beads to the representation
600 @param name the component name
601 @param ds a list of tuples corresponding to the residue ranges
603 @param colors a list of colors associated to the beads
604 @param incoord the coordinate tuple corresponding to the position
609 self.representation_is_modified =
True
611 protein_h = self.hier_dict[name]
614 colors = [self.color_dict[name]]
617 for n, dss
in enumerate(ds):
618 ds_frag = (dss[0], dss[1])
619 self.elements[name].append((dss[0], dss[1],
" ",
"bead"))
621 if ds_frag[0] == ds_frag[1]:
623 if name
in self.sequence_dict:
625 rtstr = self.onetothree[
626 self.sequence_dict[name][ds_frag[0]-1]]
633 h.set_name(name +
'_%i_bead' % (ds_frag[0]))
634 prt.set_name(name +
'_%i_bead' % (ds_frag[0]))
638 h.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
639 prt.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
640 h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
641 resolution = len(h.get_residue_indexes())
642 if "Beads" not in self.hier_representation[name]:
644 root.set_name(
"Beads")
645 self.hier_representation[name][
"Beads"] = root
646 protein_h.add_child(root)
647 self.hier_representation[name][
"Beads"].add_child(h)
649 for kk
in range(ds_frag[0], ds_frag[1] + 1):
650 self.hier_db.add_particles(name, kk, resolution, [h])
673 ptem.set_radius(radius)
676 radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
678 ptem.set_radius(radius)
680 if not tuple(incoord)
is None:
681 ptem.set_coordinates(incoord)
686 self.floppy_bodies.append(prt)
690 for p, state
in self._protocol_output:
691 p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
699 Generates a string of beads with given length.
701 self.representation_is_modified =
True
707 for chunk
in IMP.pmi.tools.list_chunks_iterator(range(begin, end + 1), length):
709 [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
713 self, name, hierarchies=
None, selection_tuples=
None,
715 resolution=0.0, num_components=10,
716 inputfile=
None, outputfile=
None,
719 covariance_type=
'full', voxel_size=1.0,
721 sampled_points=1000000, num_iter=100,
723 multiply_by_total_mass=
True,
725 intermediate_map_fn=
None,
726 density_ps_to_copy=
None,
727 use_precomputed_gaussians=
False):
729 Sets up a Gaussian Mixture Model for this component.
730 Can specify input GMM file or it will be computed.
731 @param name component name
732 @param hierarchies set up GMM for some hierarchies
733 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
734 @param particles set up GMM for particles directly
735 @param resolution usual PMI resolution for selecting particles from the hierarchies
736 @param inputfile read the GMM from this file
737 @param outputfile fit and write the GMM to this file (text)
738 @param outputmap after fitting, create GMM density file (mrc)
739 @param kernel_type for creating the intermediate density (points are sampled to make GMM). Options are IMP.em.GAUSSIAN, IMP.em.SPHERE, and IMP.em.BINARIZED_SPHERE
740 @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
741 @param voxel_size for creating the intermediate density map and output map.
742 lower number increases accuracy but also rasterizing time grows
743 @param out_hier_name name of the output density hierarchy
744 @param sampled_points number of points to sample. more will increase accuracy and fitting time
745 @param num_iter num GMM iterations. more will increase accuracy and fitting time
746 @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
747 @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
748 @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
749 @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
750 @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
758 self.representation_is_modified =
True
760 protein_h = self.hier_dict[name]
761 if "Densities" not in self.hier_representation[name]:
763 root.set_name(
"Densities")
764 self.hier_representation[name][
"Densities"] = root
765 protein_h.add_child(root)
768 fragment_particles = []
769 if not particles
is None:
770 fragment_particles += particles
771 if not hierarchies
is None:
773 self, resolution=resolution,
774 hierarchies=hierarchies)
775 if not selection_tuples
is None:
776 for st
in selection_tuples:
777 fragment_particles += IMP.pmi.tools.select_by_tuple(
778 self, tupleselection=st,
779 resolution=resolution,
780 name_is_ambiguous=
False)
783 density_particles = []
786 inputfile, density_particles,
787 self.model, transform)
788 elif density_ps_to_copy:
789 for ip
in density_ps_to_copy:
795 density_particles.append(p)
796 elif use_precomputed_gaussians:
797 if len(fragment_particles) == 0:
798 print(
"add_component_density: no particle was selected")
800 for p
in fragment_particles:
803 raise Exception(
"The particles you selected must be Fragments and XYZs")
804 nres=len(
IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
805 pos=
IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
810 raise Exception(
"We haven't created a bead file for",nres,
"residues")
814 self.model, transform)
817 if len(fragment_particles) == 0:
818 print(
"add_component_density: no particle was selected")
821 density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
830 multiply_by_total_mass,
836 s0.set_name(out_hier_name)
837 self.hier_representation[name][
"Densities"].add_child(s0)
839 for nps, p
in enumerate(density_particles):
841 p.set_name(s0.get_name() +
'_gaussian_%i' % nps)
844 def get_component_density(self, name):
845 return self.hier_representation[name][
"Densities"]
848 selection_tuples=
None,
853 '''Decorates all specified particles as Gaussians directly.
854 @param name component name
855 @param hierarchies set up GMM for some hierarchies
856 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
857 @param particles set up GMM for particles directly
858 @param resolution usual PMI resolution for selecting particles from the hierarchies
864 from math
import sqrt
865 self.representation_is_modified =
True
867 if particles
is None:
868 fragment_particles = []
870 fragment_particles = particles
872 if not hierarchies
is None:
874 self, resolution=resolution,
875 hierarchies=hierarchies)
877 if not selection_tuples
is None:
878 for st
in selection_tuples:
879 fragment_particles += IMP.pmi.tools.select_by_tuple(
880 self, tupleselection=st,
881 resolution=resolution,
882 name_is_ambiguous=
False)
884 if len(fragment_particles) == 0:
885 print(
"add all atom densities: no particle was selected")
889 print(
'setting up all atom gaussians num_particles',len(fragment_particles))
890 for n,p
in enumerate(fragment_particles):
898 print(
'setting up particle',p.get_name(),
" as individual gaussian particle")
900 if not output_map
is None:
901 print(
'writing map to', output_map)
909 Make a copy of a hierarchy and append it to a component.
912 self.representation_is_modified =
True
913 protein_h = self.hier_dict[name]
914 hierclone = IMP.atom.create_clone(hierarchy)
915 hierclone.set_name(hierclone.get_name() +
"_clone")
916 protein_h.add_child(hierclone)
917 outhier.append(hierclone)
923 for n, pmain
in enumerate(psmain):
929 self.hier_db.add_particles(
945 def dump_particle_descriptors(self):
951 particles_attributes={}
952 floppy_body_attributes={}
958 hparent=h.get_parent()
959 while hparent != hroot:
960 hparent=h.get_parent()
961 name+=
"|"+hparent.get_name()
963 particles_attributes[name]={
"COORDINATES":numpy.array(
IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
969 rigid_body_attributes={}
970 for rb
in self.rigid_bodies:
972 rf=rb.get_reference_frame()
973 t=rf.get_transformation_to()
974 trans=t.get_translation()
976 rigid_body_attributes[name]={
"TRANSLATION":numpy.array(trans),
977 "ROTATION":numpy.array(rot.get_quaternion()),
978 "COORDINATES_NONRIGID_MEMBER":{},
979 "COORDINATES_RIGID_MEMBER":{}}
980 for mi
in rb.get_member_indexes():
981 rm=self.model.get_particle(mi)
983 name_part=rm.get_name()
985 rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
987 name_part=rm.get_name()
989 rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
993 pickle.dump(particles_attributes,
994 open(
"particles_attributes.pkl",
"wb"))
995 pickle.dump(rigid_body_attributes,
996 open(
"rigid_body_attributes.pkl",
"wb"))
1000 def load_particle_descriptors(self):
1006 particles_attributes = pickle.load(open(
"particles_attributes.pkl",
1008 rigid_body_attributes = pickle.load(open(
"rigid_body_attributes.pkl",
1018 hparent=h.get_parent()
1019 while hparent != hroot:
1020 hparent=h.get_parent()
1021 name+=
"|"+hparent.get_name()
1025 xyzr.set_coordinates(particles_attributes[name][
"COORDINATES"])
1031 for rb
in self.rigid_bodies:
1033 trans=rigid_body_attributes[name][
"TRANSLATION"]
1034 rot=rigid_body_attributes[name][
"ROTATION"]
1037 rb.set_reference_frame(rf)
1038 coor_nrm_ref=rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"]
1039 coor_rm_ref_dict=rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"]
1042 for mi
in rb.get_member_indexes():
1043 rm=self.model.get_particle(mi)
1045 name_part=rm.get_name()
1046 xyz=coor_nrm_ref[name_part]
1048 self.model.set_attribute(fk, rm,xyz[n])
1050 name_part=rm.get_name()
1052 coor_rm_model.append(
IMP.core.XYZ(rm).get_coordinates())
1053 if len(coor_rm_model)==0:
continue
1059 def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1060 """Print a warning if particle names in RMF and model don't match"""
1061 def match_any_suffix():
1063 suffixes = [
"pdb",
"bead_floppy_body_rigid_body_member_floppy_body",
1064 "bead_floppy_body_rigid_body_member",
1067 if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1069 if rmfname != reprname
and not match_any_suffix():
1070 print(
"set_coordinates_from_rmf: WARNING rmf particle and "
1071 "representation particle names don't match %s %s"
1072 % (rmfname, reprname))
1076 rmf_component_name=
None,
1077 check_number_particles=
True,
1078 representation_name_to_rmf_name_map=
None,
1080 skip_gaussian_in_rmf=
False,
1081 skip_gaussian_in_representation=
False,
1083 force_rigid_update=
False):
1084 '''Read and replace coordinates from an RMF file.
1085 Replace the coordinates of particles with the same name.
1086 It assumes that the RMF and the representation have the particles
1088 @param component_name Component name
1089 @param rmf_component_name Name of the component in the RMF file
1090 (if not specified, use `component_name`)
1091 @param representation_name_to_rmf_name_map a dictionary that map
1092 the original rmf particle name to the recipient particle component name
1093 @param save_file: save a file with the names of particles of the component
1094 @param force_rigid_update: update the coordinates of rigid bodies
1095 (normally this should be called before rigid bodies are set up)
1100 prots = IMP.pmi.analysis.get_hiers_from_rmf(
1106 raise ValueError(
"cannot read hierarchy from rmf")
1110 if force_rigid_update:
1121 p, self.hier_dict.keys())
1122 if (protname
is None)
and (rmf_component_name
is not None):
1124 p, rmf_component_name)
1125 if (skip_gaussian_in_rmf):
1128 if (rmf_component_name
is not None)
and (protname == rmf_component_name):
1130 elif (rmf_component_name
is None)
and (protname == component_name):
1134 if (skip_gaussian_in_representation):
1145 reprnames=[p.get_name()
for p
in psrepr]
1146 rmfnames=[p.get_name()
for p
in psrmf]
1149 fl=open(component_name+
".txt",
"w")
1150 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1153 if check_number_particles
and not representation_name_to_rmf_name_map:
1154 if len(psrmf) != len(psrepr):
1155 fl=open(component_name+
".txt",
"w")
1156 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1157 raise ValueError(
"%s cannot proceed the rmf and the representation don't have the same number of particles; "
1158 "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1161 if not representation_name_to_rmf_name_map:
1162 for n, prmf
in enumerate(psrmf):
1164 prmfname = prmf.get_name()
1165 preprname = psrepr[n].get_name()
1166 if force_rigid_update:
1172 raise ValueError(
"component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1174 self._compare_rmf_repr_names(prmfname, preprname,
1183 rbm.set_internal_coordinates(xyz)
1185 rb = rbm.get_rigid_body()
1187 raise ValueError(
"Cannot handle nested "
1189 rb.set_reference_frame_lazy(tr)
1196 g=gprmf.get_gaussian()
1197 grepr.set_gaussian(g)
1200 repr_name_particle_map={}
1201 rmf_name_particle_map={}
1203 rmf_name_particle_map[p.get_name()]=p
1207 for prepr
in psrepr:
1209 prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1211 print(
"set_coordinates_from_rmf: WARNING missing particle name in representation_name_to_rmf_name_map, skipping...")
1220 If the root hierarchy does not exist, construct it.
1222 if "Res:" + str(int(resolution))
not in self.hier_representation[name]:
1224 root.set_name(name +
"_Res:" + str(int(resolution)))
1225 self.hier_representation[name][
1226 "Res:" + str(int(resolution))] = root
1227 protein_h.add_child(root)
1230 input_hierarchy, protein_h, type, color):
1232 Generate all needed coarse grained layers.
1234 @param name name of the protein
1235 @param resolutions list of resolutions
1236 @param protein_h root hierarchy
1237 @param input_hierarchy hierarchy to coarse grain
1238 @param type a string, typically "pdb" or "helix"
1242 if (1
in resolutions)
or (0
in resolutions):
1245 if 1
in resolutions:
1248 s1.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1250 self.hier_representation[name][
"Res:1"].add_child(s1)
1252 if 0
in resolutions:
1255 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1257 self.hier_representation[name][
"Res:0"].add_child(s0)
1260 if not isnucleicacid:
1263 atom_type=IMP.atom.AT_CA)
1267 atom_type=IMP.atom.AT_P)
1269 for p
in sel.get_selected_particles():
1271 if 0
in resolutions:
1273 resclone0 = IMP.atom.create_clone(resobject)
1275 s0.add_child(resclone0)
1276 self.hier_db.add_particles(
1280 resclone0.get_children())
1282 chil = resclone0.get_children()
1291 if 1
in resolutions:
1293 resclone1 = IMP.atom.create_clone_one(resobject)
1295 s1.add_child(resclone1)
1296 self.hier_db.add_particles(name, resindex, 1, [resclone1])
1300 prt = resclone1.get_particle()
1301 prt.set_name(
'%s_%i_%s' % (name, resindex, type))
1323 for r
in resolutions:
1324 if r != 0
and r != 1:
1330 chil = s.get_children()
1332 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1337 self.hier_representation[name][
1338 "Res:" + str(int(r))].add_child(s0)
1347 prt.set_name(
'%s_%i_%s' % (name, first, type))
1349 prt.set_name(
'%s_%i_%i_%s' % (name, first, last, type))
1351 self.hier_db.add_particles(name, kk, r, [prt])
1372 Get the hierarchies at the given resolution.
1374 The map between resolution and hierarchies is cached to accelerate
1375 the selection algorithm. The cache is invalidated when the
1376 representation was changed by any add_component_xxx.
1379 if self.representation_is_modified:
1380 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1382 self.hier_resolution[resolution] = rhbr
1383 self.representation_is_modified =
False
1386 if resolution
in self.hier_resolution:
1387 return self.hier_resolution[resolution]
1389 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1391 self.hier_resolution[resolution] = rhbr
1395 self, max_translation=300., max_rotation=2.0 * pi,
1396 avoidcollision=
True, cutoff=10.0, niterations=100,
1398 excluded_rigid_bodies=
None,
1399 ignore_initial_coordinates=
False,
1400 hierarchies_excluded_from_collision=
None):
1402 Shuffle configuration; used to restart the optimization.
1404 The configuration of the system is initialized by placing each
1405 rigid body and each bead randomly in a box with a side of
1406 `max_translation` angstroms, and far enough from each other to
1407 prevent any steric clashes. The rigid bodies are also randomly rotated.
1409 @param avoidcollision check if the particle/rigid body was
1410 placed close to another particle; uses the optional
1411 arguments cutoff and niterations
1412 @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1415 if excluded_rigid_bodies
is None:
1416 excluded_rigid_bodies = []
1417 if hierarchies_excluded_from_collision
is None:
1418 hierarchies_excluded_from_collision = []
1420 if len(self.rigid_bodies) == 0:
1421 print(
"shuffle_configuration: rigid bodies were not intialized")
1424 gcpf.set_distance(cutoff)
1426 hierarchies_excluded_from_collision_indexes = []
1435 if bounding_box
is not None:
1436 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1441 for h
in hierarchies_excluded_from_collision:
1445 allparticleindexes = list(
1446 set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1448 print(hierarchies_excluded_from_collision)
1449 print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1451 print(
'shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies),
'rigid bodies')
1452 for rb
in self.rigid_bodies:
1453 if rb
not in excluded_rigid_bodies:
1455 rbindexes = rb.get_member_particle_indexes()
1457 set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1458 otherparticleindexes = list(
1459 set(allparticleindexes) - set(rbindexes))
1461 if len(otherparticleindexes)
is None:
1465 while niter < niterations:
1466 if (ignore_initial_coordinates):
1472 if bounding_box
is not None:
1488 gcpf.get_close_pairs(
1490 otherparticleindexes,
1494 if (ignore_initial_coordinates):
1495 print (rb.get_name(),
IMP.core.XYZ(rb).get_coordinates())
1498 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1499 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
1500 if niter == niterations:
1501 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1505 print(
'shuffling', len(self.floppy_bodies),
'floppy bodies')
1506 for fb
in self.floppy_bodies:
1507 if (avoidcollision):
1512 otherparticleindexes = list(
1513 set(allparticleindexes) - set(fbindexes))
1514 if len(otherparticleindexes)
is None:
1532 while niter < niterations:
1533 if (ignore_initial_coordinates):
1539 if bounding_box
is not None:
1551 if (avoidcollision):
1554 gcpf.get_close_pairs(
1556 otherparticleindexes,
1560 if (ignore_initial_coordinates):
1561 print (fb.get_name(),
IMP.core.XYZ(fb).get_coordinates())
1564 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1565 print(
"shuffle_configuration: floppy body name: " + fb.get_name())
1566 if niter == niterations:
1567 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1571 def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1575 self.reference_structures[label] = (
1579 def get_all_rmsds(self):
1582 for label
in self.reference_structures:
1584 current_coordinates = [
IMP.core.XYZ(p).get_coordinates()
1585 for p
in self.reference_structures[label][1]]
1586 reference_coordinates = self.reference_structures[label][0]
1587 if len(reference_coordinates) != len(current_coordinates):
1588 print(
"calculate_all_rmsds: reference and actual coordinates are not the same")
1591 current_coordinates,
1592 reference_coordinates)
1594 reference_coordinates,
1595 current_coordinates)
1599 reference_coordinates,
1600 current_coordinates)
1601 rmsds[label +
"_GlobalRMSD"] = rmsd_global
1602 rmsds[label +
"_RelativeDRMS"] = rmsd_relative
1605 def setup_component_geometry(self, name, color=None, resolution=1.0):
1607 color = self.color_dict[name]
1611 self.hier_geometry_pairs[name] = []
1612 protein_h = self.hier_dict[name]
1614 pbr = IMP.pmi.tools.sort_by_residues(pbr)
1616 for n
in range(len(pbr) - 1):
1617 self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1620 self, name, resolution=10, scale=1.0):
1622 Generate restraints between contiguous fragments in the hierarchy.
1623 The linkers are generated at resolution 10 by default.
1630 protein_h = self.hier_dict[name]
1632 frs = self.hier_db.get_preroot_fragments_by_resolution(
1638 start = fr.get_children()[0]
1643 end = fr.get_children()[-1]
1649 SortedSegments.append((start, end, startres))
1650 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1653 for x
in range(len(SortedSegments) - 1):
1654 last = SortedSegments[x][1]
1655 first = SortedSegments[x + 1][0]
1662 residuegap = firstresn - lastresn - 1
1663 if self.disorderedlength
and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1666 optdist = sqrt(5 / 3) * 1.93 * \
1667 (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1669 if self.upperharmonic:
1675 optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1676 if self.upperharmonic:
1682 pt0 = last.get_particle()
1683 pt1 = first.get_particle()
1686 print(
"Adding sequence connectivity restraint between", pt0.get_name(),
" and ", pt1.get_name(),
'of distance', optdist)
1687 sortedsegments_cr.add_restraint(r)
1688 self.linker_restraints_dict[
1689 "LinkerRestraint-" + pt0.get_name() +
"-" + pt1.get_name()] = r
1690 self.connected_intra_pairs.append((pt0, pt1))
1691 self.connected_intra_pairs.append((pt1, pt0))
1693 self.linker_restraints.add_restraint(sortedsegments_cr)
1694 self.linker_restraints.add_restraint(unmodeledregions_cr)
1696 self.sortedsegments_cr_dict[name] = sortedsegments_cr
1697 self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1699 def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1701 pts = IMP.pmi.tools.ParticleToSampleList()
1702 for n, fb
in enumerate(self.floppy_bodies):
1703 pts.add_particle(fb,
"Floppy_Bodies", 1.0,
"Floppy_Body_" + str(n))
1704 if len(pts.get_particles_to_sample()) > 0:
1706 print(
"optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1709 print(
"optimize_floppy_bodies: no particle to optimize")
1712 nSymmetry=
None, skip_gaussian_in_clones=
False):
1714 The copies must not contain rigid bodies.
1715 The symmetry restraints are applied at each leaf.
1719 self.representation_is_modified =
True
1720 ncopies = len(copies) + 1
1723 for k
in range(len(copies)):
1724 if (nSymmetry
is None):
1725 rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1728 rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1730 rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1738 for n, p
in enumerate(main_hiers):
1739 if (skip_gaussian_in_clones):
1748 self.model.add_score_state(c)
1749 print(
"Completed setting " + str(maincopy) +
" as a reference for "
1750 + str(copies[k]) +
" by rotating it by "
1751 + str(rotation_angle / 2.0 / pi * 360)
1752 +
" degrees around the " + str(rotational_axis) +
" axis.")
1755 def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1758 self.representation_is_modified =
True
1760 mainparticles = particles_reference
1762 t=initial_transformation
1764 p.set_name(
"RigidBody_Symmetry")
1770 copyparticles = particles_copy
1774 for n, p
in enumerate(mainparticles):
1776 pc = copyparticles[n]
1778 mainpurged.append(p)
1784 copypurged.append(pc)
1791 for n, p
in enumerate(mainpurged):
1794 print(
"setting " + p.get_name() +
" as reference for " + pc.get_name())
1797 lc.add(pc.get_index())
1800 self.model.add_score_state(c)
1803 self.rigid_bodies.append(rb)
1804 self.rigid_body_symmetries.append(rb)
1805 rb.set_name(label+
".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1808 def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1809 longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1812 self.representation_is_modified =
True
1815 protein_h = self.hier_dict[maincopy]
1818 for ilc
in range(-longitudinal_copies, longitudinal_copies + 1):
1819 for iac
in range(axial_copies):
1820 copyname = maincopy +
"_a" + str(ilc) +
"_l" + str(iac)
1821 self.create_component(copyname, 0.0)
1822 for hier
in protein_h.get_children():
1824 copyhier = self.hier_dict[copyname]
1825 outhiers.append(copyhier)
1829 2 * pi / axial_copies * (float(iac)))
1830 translation_vector = tuple(
1831 [translation_value * float(ilc) * x
for x
in axis])
1832 print(translation_vector)
1837 for n, p
in enumerate(mainparts):
1844 lc.add(pc.get_index())
1846 self.model.add_score_state(c)
1852 Load coordinates in the current representation.
1853 This should be done only if the hierarchy self.prot is identical
1854 to the one as stored in the rmf i.e. all components were added.
1859 rh = RMF.open_rmf_file_read_only(rmfname)
1867 create the representation (i.e. hierarchies) from the rmf file.
1868 it will be stored in self.prot, which will be overwritten.
1869 load the coordinates from the rmf file at frameindex.
1871 rh = RMF.open_rmf_file_read_only(rmfname)
1877 for p
in self.prot.get_children():
1878 self.create_component(p.get_name())
1879 self.hier_dict[p.get_name()] = p
1881 still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1882 save floppy bodies in self.floppy_bodies
1883 get the connectivity restraints
1888 Construct a rigid body from hierarchies (and optionally particles).
1890 @param hiers list of hierarchies to make rigid
1891 @param particles list of particles to add to the rigid body
1894 if particles
is None:
1897 rigid_parts = set(particles)
1900 print(
"set_rigid_body_from_hierarchies> setting up a new rigid body")
1907 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1910 name += hier.get_name() +
"-"
1911 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1913 if len(list(rigid_parts)) != 0:
1915 rb.set_coordinates_are_optimized(
True)
1916 rb.set_name(name +
"rigid_body")
1917 self.rigid_bodies.append(rb)
1921 print(
"set_rigid_body_from_hierarchies> rigid body could not be setup")
1925 Construct a rigid body from a list of subunits.
1927 Each subunit is a tuple that identifies the residue ranges and the
1928 component name (as used in create_component()).
1930 subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1932 [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1934 example: ["prot1","prot2",("prot3",(1,10))]
1936 sometimes, we know about structure of an interaction
1937 and here we make such PPIs rigid
1942 if type(s) == type(tuple())
and len(s) == 2:
1946 residue_indexes=list(range(s[1][0],
1948 if len(sel.get_selected_particles()) == 0:
1949 print(
"set_rigid_bodies: selected particle does not exist")
1950 for p
in sel.get_selected_particles():
1954 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1958 elif type(s) == type(str()):
1960 if len(sel.get_selected_particles()) == 0:
1961 print(
"set_rigid_bodies: selected particle does not exist")
1962 for p
in sel.get_selected_particles():
1966 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1971 rb.set_coordinates_are_optimized(
True)
1972 rb.set_name(
''.join(str(subunits)) +
"_rigid_body")
1973 self.rigid_bodies.append(rb)
1976 def set_super_rigid_body_from_hierarchies(
1982 super_rigid_xyzs = set()
1983 super_rigid_rbs = set()
1985 print(
"set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
1992 super_rigid_rbs.add(rb)
1994 super_rigid_xyzs.add(p)
1995 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1996 if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
1999 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2002 self.super_rigid_bodies.append(
2003 (super_rigid_xyzs, super_rigid_rbs, axis))
2005 def fix_rigid_bodies(self, rigid_bodies):
2006 self.fixed_rigid_bodies += rigid_bodies
2010 self, hiers, min_length=
None, max_length=
None, axis=
None):
2012 Make a chain of super rigid bodies from a list of hierarchies.
2014 Takes a linear list of hierarchies (they are supposed to be
2015 sequence-contiguous) and produces a chain of super rigid bodies
2016 with given length range, specified by min_length and max_length.
2019 hiers = IMP.pmi.tools.flatten_list(hiers)
2022 for hs
in IMP.pmi.tools.sublist_iterator(hiers, min_length, max_length):
2023 self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2025 def set_super_rigid_bodies(self, subunits, coords=None):
2026 super_rigid_xyzs = set()
2027 super_rigid_rbs = set()
2030 if type(s) == type(tuple())
and len(s) == 3:
2034 residue_indexes=list(range(s[0],
2036 if len(sel.get_selected_particles()) == 0:
2037 print(
"set_rigid_bodies: selected particle does not exist")
2038 for p
in sel.get_selected_particles():
2041 super_rigid_rbs.add(rb)
2043 super_rigid_xyzs.add(p)
2044 elif type(s) == type(str()):
2046 if len(sel.get_selected_particles()) == 0:
2047 print(
"set_rigid_bodies: selected particle does not exist")
2048 for p
in sel.get_selected_particles():
2052 super_rigid_rbs.add(rb)
2054 super_rigid_xyzs.add(p)
2055 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2059 Remove leaves of hierarchies from the floppy bodies list based
2060 on the component name
2063 ps=[h.get_particle()
for h
in hs]
2064 for p
in self.floppy_bodies:
2066 if p
in ps: self.floppy_bodies.remove(p)
2067 if p
in hs: self.floppy_bodies.remove(p)
2074 Remove leaves of hierarchies from the floppy bodies list.
2076 Given a list of hierarchies, get the leaves and remove the
2077 corresponding particles from the floppy bodies list. We need this
2078 function because sometimes
2079 we want to constrain the floppy bodies in a rigid body - for instance
2080 when you want to associate a bead with a density particle.
2082 for h
in hierarchies:
2085 if p
in self.floppy_bodies:
2086 print(
"remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2087 self.floppy_bodies.remove(p)
2090 def set_floppy_bodies(self):
2091 for p
in self.floppy_bodies:
2093 p.set_name(name +
"_floppy_body")
2095 print(
"I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2098 rb.set_is_rigid_member(p.get_index(),
False)
2101 rb.set_is_rigid_member(p.get_particle_index(),
False)
2102 p.set_name(p.get_name() +
"_rigid_body_member")
2104 def set_floppy_bodies_from_hierarchies(self, hiers):
2109 self.floppy_bodies.append(p)
2116 selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2117 @return the particles
2120 print(selection_tuples)
2121 for s
in selection_tuples:
2122 ps = IMP.pmi.tools.select_by_tuple(
2123 representation=self, tupleselection=tuple(s),
2124 resolution=
None, name_is_ambiguous=
False)
2128 def get_connected_intra_pairs(self):
2129 return self.connected_intra_pairs
2131 def set_rigid_bodies_max_trans(self, maxtrans):
2132 self.maxtrans_rb = maxtrans
2134 def set_rigid_bodies_max_rot(self, maxrot):
2135 self.maxrot_rb = maxrot
2137 def set_super_rigid_bodies_max_trans(self, maxtrans):
2138 self.maxtrans_srb = maxtrans
2140 def set_super_rigid_bodies_max_rot(self, maxrot):
2141 self.maxrot_srb = maxrot
2143 def set_floppy_bodies_max_trans(self, maxtrans):
2144 self.maxtrans_fb = maxtrans
2148 Fix rigid bodies in their actual position.
2149 The get_particles_to_sample() function will return
2150 just the floppy bodies (if they are not fixed).
2152 self.rigidbodiesarefixed = rigidbodiesarefixed
2156 Fix floppy bodies in their actual position.
2157 The get_particles_to_sample() function will return
2158 just the rigid bodies (if they are not fixed).
2160 self.floppybodiesarefixed=floppybodiesarefixed
2162 def draw_hierarchy_graph(self):
2164 print(
"Drawing hierarchy graph for " + c.get_name())
2165 IMP.pmi.output.get_graph_from_hierarchy(c)
2167 def get_geometries(self):
2170 for name
in self.hier_geometry_pairs:
2171 for pt
in self.hier_geometry_pairs[name]:
2185 def setup_bonds(self):
2188 for name
in self.hier_geometry_pairs:
2189 for pt
in self.hier_geometry_pairs[name]:
2204 def show_component_table(self, name):
2205 if name
in self.sequence_dict:
2206 lastresn = len(self.sequence_dict[name])
2209 residues = self.hier_db.get_residue_numbers(name)
2210 firstresn = min(residues)
2211 lastresn = max(residues)
2213 for nres
in range(firstresn, lastresn + 1):
2215 resolutions = self.hier_db.get_residue_resolutions(name, nres)
2217 for r
in resolutions:
2218 ps += self.hier_db.get_particles(name, nres, r)
2219 print(
"%20s %7s" % (name, nres),
" ".join([
"%20s %7s" % (str(p.get_name()),
2222 print(
"%20s %20s" % (name, nres),
"**** not represented ****")
2224 def draw_hierarchy_composition(self):
2226 ks = list(self.elements.keys())
2230 for k
in self.elements:
2231 for l
in self.elements[k]:
2236 self.draw_component_composition(k, max)
2238 def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2239 from matplotlib
import pyplot
2240 import matplotlib
as mpl
2242 tmplist = sorted(self.elements[k], key=itemgetter(0))
2244 endres = tmplist[-1][1]
2246 print(
"draw_component_composition: missing information for component %s" % name)
2248 fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2249 ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2254 norm = mpl.colors.Normalize(vmin=5, vmax=10)
2258 for n, l
in enumerate(tmplist):
2262 if bounds[-1] != l[0]:
2263 colors.append(
"white")
2266 colors.append(
"#99CCFF")
2268 colors.append(
"#FFFF99")
2270 colors.append(
"#33CCCC")
2272 bounds.append(l[1] + 1)
2275 colors.append(
"#99CCFF")
2277 colors.append(
"#FFFF99")
2279 colors.append(
"#33CCCC")
2281 bounds.append(l[1] + 1)
2283 if bounds[-1] - 1 == l[0]:
2287 colors.append(
"white")
2290 bounds.append(bounds[-1])
2291 colors.append(
"white")
2292 cmap = mpl.colors.ListedColormap(colors)
2293 cmap.set_over(
'0.25')
2294 cmap.set_under(
'0.75')
2296 norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2297 cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2303 spacing=
'proportional',
2304 orientation=
'horizontal')
2313 mid = 1.0 / endres * float(l[0])
2318 l[2], xy=(mid, 1), xycoords=
'axes fraction',
2319 xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords=
'axes fraction',
2320 arrowprops=dict(arrowstyle=
"->",
2321 connectionstyle=
"angle,angleA=0,angleB=90,rad=10"),
2323 extra_artists.append(t)
2326 title = ax.text(-0.005, 0.5, k, ha=
"right", va=
"center", rotation=90,
2329 extra_artists.append(title)
2331 labels = len(bounds) * [
" "]
2332 ax.set_xticklabels(labels)
2333 mid = 1.0 / endres * float(bounds[0])
2334 t = ax.annotate(bounds[0], xy=(mid, 0), xycoords=
'axes fraction',
2335 xytext=(mid - 0.01, -0.5), textcoords=
'axes fraction',)
2336 extra_artists.append(t)
2337 offsets = [0, -0.5, -1.0]
2339 for n
in range(1, len(bounds)):
2340 if bounds[n] == bounds[n - 1]:
2342 mid = 1.0 / endres * float(bounds[n])
2343 if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2345 offset = offsets[nclashes % 3]
2351 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2352 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction')
2355 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2356 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction', arrowprops=dict(arrowstyle=
"-"))
2357 extra_artists.append(t)
2359 cb2.add_lines(bounds, [
"black"] * len(bounds), [1] * len(bounds))
2363 k +
"structure.pdf",
2366 bbox_extra_artists=(extra_artists),
2367 bbox_inches=
'tight')
2370 def draw_coordinates_projection(self):
2371 import matplotlib.pyplot
as pp
2374 for name
in self.hier_geometry_pairs:
2375 for pt
in self.hier_geometry_pairs[name]:
2385 xpairs.append([x1, x2])
2386 ypairs.append([y1, y2])
2389 for xends, yends
in zip(xpairs, ypairs):
2394 pp.plot(xlist, ylist,
'b-', alpha=0.1)
2397 def get_prot_name_from_particle(self, particle):
2398 names = self.get_component_names()
2399 particle0 = particle
2401 while not name
in names:
2404 particle0 = h.get_particle()
2408 def get_particles_to_sample(self):
2418 if not self.rigidbodiesarefixed:
2419 for rb
in self.rigid_bodies:
2424 if rb
not in self.fixed_rigid_bodies:
2427 if not self.floppybodiesarefixed:
2428 for fb
in self.floppy_bodies:
2435 for srb
in self.super_rigid_bodies:
2438 rigid_bodies = list(srb[1])
2439 filtered_rigid_bodies = []
2440 for rb
in rigid_bodies:
2441 if rb
not in self.fixed_rigid_bodies:
2442 filtered_rigid_bodies.append(rb)
2443 srbtmp.append((srb[0], filtered_rigid_bodies))
2445 self.rigid_bodies = rbtmp
2446 self.floppy_bodies = fbtmp
2447 self.super_rigid_bodies = srbtmp
2449 ps[
"Rigid_Bodies_SimplifiedModel"] = (
2453 ps[
"Floppy_Bodies_SimplifiedModel"] = (
2456 ps[
"SR_Bodies_SimplifiedModel"] = (
2457 self.super_rigid_bodies,
2462 def set_output_level(self, level):
2463 self.output_level = level
2465 def _evaluate(self, deriv):
2466 """Evaluate the total score of all added restraints"""
2468 return r.evaluate(deriv)
2470 def get_output(self):
2474 output[
"SimplifiedModel_Total_Score_" +
2475 self.label] = str(self._evaluate(
False))
2476 output[
"SimplifiedModel_Linker_Score_" +
2477 self.label] = str(self.linker_restraints.unprotected_evaluate(
None))
2478 for name
in self.sortedsegments_cr_dict:
2479 partialscore = self.sortedsegments_cr_dict[name].evaluate(
False)
2480 score += partialscore
2482 "SimplifiedModel_Link_SortedSegments_" +
2487 partialscore = self.unmodeledregions_cr_dict[name].evaluate(
False)
2488 score += partialscore
2490 "SimplifiedModel_Link_UnmodeledRegions_" +
2495 for rb
in self.rigid_body_symmetries:
2497 output[name +
"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2498 for name
in self.linker_restraints_dict:
2503 self.linker_restraints_dict[
2504 name].unprotected_evaluate(
2507 if len(self.reference_structures.keys()) != 0:
2508 rmsds = self.get_all_rmsds()
2511 "SimplifiedModel_" +
2514 self.label] = rmsds[
2517 if self.output_level ==
"high":
2520 output[
"Coordinates_" +
2521 p.get_name() +
"_" + self.label] = str(d)
2523 output[
"_TotalScore"] = str(score)
2526 def get_test_output(self):
2528 output = self.get_output()
2529 for n, p
in enumerate(self.get_particles_to_sample()):
2530 output[
"Particle_to_sample_" + str(n)] = str(p)
2532 output[
"Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2533 output[
"Number_of_floppy_bodies"] = len(self.floppy_bodies)
2534 output[
"Number_of_rigid_bodies"] = len(self.rigid_bodies)
2535 output[
"Number_of_super_bodies"] = len(self.super_rigid_bodies)
2536 output[
"Selection_resolution_1"] = len(
2538 output[
"Selection_resolution_5"] = len(
2540 output[
"Selection_resolution_7"] = len(
2542 output[
"Selection_resolution_10"] = len(
2544 output[
"Selection_resolution_100"] = len(
2547 output[
"Selection_resolution=1"] = len(
2549 output[
"Selection_resolution=1,resid=10"] = len(
2551 for resolution
in self.hier_resolution:
2552 output[
"Hier_resolution_dictionary" +
2553 str(resolution)] = len(self.hier_resolution[resolution])
2554 for name
in self.hier_dict:
2556 "Selection_resolution=1,resid=10,name=" +
2564 "Selection_resolution=1,resid=10,name=" +
2566 ",ambiguous"] = len(
2571 name_is_ambiguous=
True,
2574 "Selection_resolution=1,resid=10,name=" +
2576 ",ambiguous"] = len(
2581 name_is_ambiguous=
True,
2584 "Selection_resolution=1,resrange=(10,20),name=" +
2593 "Selection_resolution=1,resrange=(10,20),name=" +
2595 ",ambiguous"] = len(
2600 name_is_ambiguous=
True,
2604 "Selection_resolution=10,resrange=(10,20),name=" +
2613 "Selection_resolution=10,resrange=(10,20),name=" +
2615 ",ambiguous"] = len(
2620 name_is_ambiguous=
True,
2624 "Selection_resolution=100,resrange=(10,20),name=" +
2633 "Selection_resolution=100,resrange=(10,20),name=" +
2635 ",ambiguous"] = len(
2640 name_is_ambiguous=
True,
def link_components_to_rmf
Load coordinates in the current representation.
Select non water and non hydrogen atoms.
double get_volume_from_residue_type(ResidueType rt)
Return an estimate for the volume of a given residue.
def add_metadata
Associate some metadata with this modeling.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A decorator to associate a particle with a part of a protein/DNA/RNA.
static Gaussian setup_particle(Model *m, ParticleIndex pi)
void show_molecular_hierarchy(Hierarchy h)
Print out the molecular hierarchy.
double get_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
Apply a SingletonFunction to a SingletonContainer to maintain an invariant.
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
def shuffle_configuration
Shuffle configuration; used to restart the optimization.
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.
static Provenanced setup_particle(Model *m, ParticleIndex pi, Provenance p)
static Atom setup_particle(Model *m, ParticleIndex pi, Atom other)
static Fragment setup_particle(Model *m, ParticleIndex pi)
Upper bound harmonic function (non-zero when feature > mean)
static Symmetric setup_particle(Model *m, ParticleIndex pi, Float symmetric)
def set_coordinates_from_rmf
Read and replace coordinates from an RMF file.
static XYZR setup_particle(Model *m, ParticleIndex pi)
A decorator for a particle which has bonds.
Select atoms which are selected by both selectors.
Rotation3D get_random_rotation_3d(const Rotation3D ¢er, double distance)
Pick a rotation at random near the provided one.
double get_mass(ResidueType c)
Get the mass from the residue type.
def create_non_modeled_component
Create a component that isn't used in the modeling.
Color get_rgb_color(double f)
Return the color for f from the RGB color map.
double get_mass_from_number_of_residues(unsigned int num_aa)
Estimate the mass of a protein from the number of amino acids.
Provenance create_clone(Provenance p)
Clone provenance (including previous provenance)
def set_rigid_bodies_as_fixed
Fix rigid bodies in their actual position.
def get_particles_from_selection_tuples
selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....
Add symmetric attribute to a particle.
double get_ball_radius_from_volume_3d(double volume)
Return the radius of a sphere with a given volume.
def setup_component_sequence_connectivity
Generate restraints between contiguous fragments in the hierarchy.
def create_rotational_symmetry
The copies must not contain rigid bodies.
Add uncertainty to a particle.
A score on the distance between the surfaces of two spheres.
static Residue setup_particle(Model *m, ParticleIndex pi, ResidueType t, int index, int insertion_code)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
static XYZ setup_particle(Model *m, ParticleIndex pi)
void read_pdb(TextInput input, int model, Hierarchy h)
def remove_floppy_bodies
Remove leaves of hierarchies from the floppy bodies list.
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
static Uncertainty setup_particle(Model *m, ParticleIndex pi, Float uncertainty)
def remove_floppy_bodies_from_component
Remove leaves of hierarchies from the floppy bodies list based on the component name.
Add resolution to a particle.
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
static Molecule setup_particle(Model *m, ParticleIndex pi)
Rotation3D get_rotation_about_axis(const Vector3D &axis, double angle)
Generate a Rotation3D object from a rotation around an axis.
Class for storing model, its restraints, constraints, and particles.
def create_transformed_component
Create a transformed view of a component.
Select all CB ATOM records.
A Gaussian distribution in 3D.
static bool get_is_setup(Model *m, ParticleIndex pi)
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
def deprecated_pmi1_object
Mark a PMI1 class as deprecated.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
def set_rigid_bodies
Construct a rigid body from a list of subunits.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Store a list of ParticleIndexes.
def deprecated_method
Python decorator to mark a method as deprecated.
A decorator for a particle representing an atom.
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
static Bonded setup_particle(Model *m, ParticleIndex pi)
def add_protocol_output
Capture details of the modeling protocol.
static bool get_is_setup(Model *m, ParticleIndex pi)
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def add_all_atom_densities
Decorates all specified particles as Gaussians directly.
def coarse_hierarchy
Generate all needed coarse grained layers.
def add_component_pdb
Add a component that has an associated 3D structure in a PDB file.
Track creation of a system fragment from a PDB file.
Basic utilities for handling cryo-electron microscopy 3D density maps.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
A decorator for a particle with x,y,z coordinates.
static bool get_is_setup(Model *m, ParticleIndex pi)
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
def set_chain_of_super_rigid_bodies
Make a chain of super rigid bodies from a list of hierarchies.
def add_component_sequence
Add the primary sequence for a single component.
Tools for clustering and cluster analysis.
static Resolution setup_particle(Model *m, ParticleIndex pi, Float resolution)
def create_components_from_rmf
still not working.
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Find all nearby pairs by testing all pairs.
Classes for writing output files and processing them.
Simple implementation of segments in 3D.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A decorator for a residue.
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...
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def add_component_density
Sets up a Gaussian Mixture Model for this component.
Sample using Monte Carlo.
The general base class for IMP exceptions.
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
def get_file_dataset
Get the dataset associated with a filename, or None.
IMP::core::RigidBody create_rigid_body(Hierarchy h)
static Reference setup_particle(Model *m, ParticleIndex pi, ParticleIndexAdaptor reference)
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual particles of a Model object.
Bond get_bond(Bonded a, Bonded b)
Get the bond between two particles.
std::string get_data_path(std::string file_name)
Return the full path to one of this module's data files.
Select atoms which are selected by either or both selectors.
def get_hierarchies_at_given_resolution
Get the hierarchies at the given resolution.
Store info for a chain of a protein.
Applies a PairScore to a Pair.
Select all CA ATOM records.
Python classes to represent, score, sample and analyze models.
static bool get_is_setup(Model *m, ParticleIndex pi)
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.
A decorator for a rigid body.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A dictionary-like wrapper for reading and storing sequence data.
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)
Output IMP model data in various file formats.
Functionality for loading, creating, manipulating and scoring atomic structures.
Tag part of the system to track how it was created.
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Hierarchies get_leaves(const Selection &h)
def add_component_necklace
Generates a string of beads with given length.
An exception for an invalid value being passed to IMP.
def set_file_dataset
Associate a dataset with a filename.
def set_floppy_bodies_as_fixed
Fix floppy bodies in their actual position.
Select hierarchy particles identified by the biological name.
static RigidBody setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor ps)
Select all ATOM and HETATM records with the given chain ids.
Support for the RMF file format for storing hierarchical molecular data and markup.
static bool get_is_setup(Model *m, ParticleIndex pi)
def set_rigid_body_from_hierarchies
Construct a rigid body from hierarchies (and optionally particles).
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Transformation3D get_random_local_transformation(Vector3D origin, double max_translation=5., double max_angle_in_rad=0.26)
Get a local transformation.
def check_root
If the root hierarchy does not exist, construct it.
def add_component_beads
add beads to the representation
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
def add_component_hierarchy_clone
Make a copy of a hierarchy and append it to a component.
Harmonic function (symmetric about the mean)
A decorator for a particle with x,y,z coordinates and a radius.