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
27 def __init__(self, doi, root):
31 def get_fname(self, fname):
32 """Return a path relative to the top of the repository"""
33 return os.path.relpath(fname, self._root)
35 class _StateInfo(object):
36 """Score state-specific information about this representation."""
41 "use IMP.pmi.topology.System instead")
42 class Representation(object):
46 Set up the representation of all proteins and nucleic acid macromolecules.
48 Create the molecular hierarchies, representation,
49 sequence connectivity for the various involved proteins and
50 nucleic acid macromolecules:
52 Create a protein, DNA or RNA, represent it as a set of connected balls of appropriate
53 radii and number of residues, PDB at given resolution(s), or ideal helices.
55 How to use the SimplifiedModel class (typical use):
57 see test/test_hierarchy_contruction.py
61 1) Create a chain of helices and flexible parts
63 c_1_119 =self.add_component_necklace("prot1",1,119,20)
64 c_120_131 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(120,131))
65 c_132_138 =self.add_component_beads("prot1",[(132,138)])
66 c_139_156 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(139,156))
67 c_157_174 =self.add_component_beads("prot1",[(157,174)])
68 c_175_182 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(175,182))
69 c_183_194 =self.add_component_beads("prot1",[(183,194)])
70 c_195_216 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(195,216))
71 c_217_250 =self.add_component_beads("prot1",[(217,250)])
74 self.set_rigid_body_from_hierarchies(c_120_131)
75 self.set_rigid_body_from_hierarchies(c_139_156)
76 self.set_rigid_body_from_hierarchies(c_175_182)
77 self.set_rigid_body_from_hierarchies(c_195_216)
79 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,
82 self.set_chain_of_super_rigid_bodies(clist,2,3)
84 self.set_super_rigid_bodies(["prot1"])
88 def __init__(self, model, upperharmonic=True, disorderedlength=True):
90 @param model the model
91 @param upperharmonic This flag uses either harmonic (False)
92 or upperharmonic (True) in the intra-pair
93 connectivity restraint.
94 @param disorderedlength This flag uses either disordered length
95 calculated for random coil peptides (True) or zero
96 surface-to-surface distance between beads (False)
97 as optimal distance for the sequence connectivity
101 self.state = _StateInfo()
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())
175 """Associate some metadata with this modeling.
176 @param m an instance of an ihm metadata class, such as
177 ihm.Software, ihm.Citation, or ihm.location.Repository.
179 self._metadata.append(m)
182 """Capture details of the modeling protocol.
183 @param p an instance of IMP.pmi.output.ProtocolOutput or a subclass.
185 state = p._add_state(self)
186 self._protocol_output.append((p, state))
187 p._each_metadata.append(self._metadata)
188 state.model = self.model
189 state.prot = self.prot
190 protocol_output = property(
lambda self:
191 [x[0]
for x
in self._protocol_output])
193 def set_label(self, label):
196 def create_component(self, name, color=0.0):
203 protein_h.set_name(name)
204 self.hier_dict[name] = protein_h
205 self.hier_representation[name] = {}
206 self.hier_db.add_name(name)
207 self.prot.add_child(protein_h)
208 self.color_dict[name] = color
209 self.elements[name] = []
210 for p, state
in self._protocol_output:
211 p.create_component(state, name,
True)
214 """Create a transformed view of a component.
215 This does not copy the coordinates or representation of the
216 component, and the transformed component is not visible to IMP,
217 but when the model is written to a ProtocolOutput, the transformed
218 component is added. This can be used to easily add symmetry-related
219 'copies' of a component without copying the underlying
221 for p, state
in self._protocol_output:
222 p.create_transformed_component(state, name, original, transform)
225 """Create a component that isn't used in the modeling.
226 No coordinates or other structural data for this component will
227 be read or written, but a primary sequence can be assigned. This
228 is useful if the input experimental data is of a system larger
229 than that modeled. Any references to these non-modeled components
230 can then be correctly resolved later."""
231 self._non_modeled_components[name] =
None
232 self.elements[name] = []
233 for p, state
in self._protocol_output:
234 p.create_component(state, name,
False)
236 def get_component_names(self):
237 return list(self.hier_dict.keys())
242 Add the primary sequence for a single component.
244 @param name Human-readable name of the component
245 @param filename Name of the FASTA file
246 @param id Identifier of the sequence in the FASTA file header
247 (if not provided, use `name` instead)
252 if id
not in record_dict:
253 raise KeyError(
"id %s not found in fasta file" % id)
254 length = len(record_dict[id])
255 self.sequence_dict[name] = str(record_dict[id])
257 if name
not in self._non_modeled_components:
258 protein_h = self.hier_dict[name]
259 protein_h.set_sequence(self.sequence_dict[name])
262 self.sequence_dict[name]=offs_str+self.sequence_dict[name]
264 self.elements[name].append((length, length,
" ",
"end"))
265 for p, state
in self._protocol_output:
266 p.add_component_sequence(state, name, self.sequence_dict[name])
268 def autobuild_model(self, name, pdbname, chain,
269 resolutions=
None, resrange=
None,
271 color=
None, pdbresrange=
None, offset=0,
272 show=
False, isnucleicacid=
False,
275 self.representation_is_modified =
True
279 color = self.color_dict[name]
281 self.color_dict[name] = color
283 if resolutions
is None:
285 print(
"autobuild_model: constructing %s from pdb %s and chain %s" % (name, pdbname, str(chain)))
294 t.get_children()[0].get_children()[0]).
get_index()
296 t.get_children()[0].get_children()[-1]).
get_index()
302 if name
in self.sequence_dict:
303 resrange = (1, len(self.sequence_dict[name]))
305 resrange = (start + offset, end + offset)
307 if resrange[1]
in (-1,
'END'):
308 resrange = (resrange[0],end)
309 start = resrange[0] - offset
310 end = resrange[1] - offset
329 for n, g
in enumerate(gaps):
333 print(
"autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
335 chain, resolutions=resolutions,
336 color=color, cacenters=
True,
337 resrange=(first, last),
338 offset=offset, isnucleicacid=isnucleicacid)
339 elif g[2] ==
"gap" and n > 0:
340 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
341 parts = self.hier_db.get_particles_at_closest_resolution(name,
346 first+offset, last+offset, missingbeadsize, incoord=xyz)
348 elif g[2] ==
"gap" and n == 0:
350 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
352 first+offset, last+offset, missingbeadsize, incoord=xyznter)
357 resrange=
None, offset=0, cacenters=
True, show=
False,
358 isnucleicacid=
False, readnonwateratoms=
False,
359 read_ca_cb_only=
False):
361 Add a component that has an associated 3D structure in a PDB file.
363 Reads the PDB, and constructs the fragments corresponding to contiguous
366 @return a list of hierarchies.
368 @param name (string) the name of the component
369 @param pdbname (string) the name of the PDB file
370 @param chain (string or integer) can be either a string (eg, "A")
371 or an integer (eg, 0 or 1) in case you want
372 to get the corresponding chain number in the PDB.
373 @param resolutions (integers) a list of integers that corresponds
374 to the resolutions that have to be generated
375 @param color (float from 0 to 1) the color applied to the
376 hierarchies generated
377 @param resrange (tuple of integers): the residue range to extract
378 from the PDB. It is a tuple (beg,end). If not specified,
379 all residues belonging to the specified chain are read.
380 @param offset (integer) specifies the residue index offset to be
381 applied when reading the PDB (the FASTA sequence is
382 assumed to start from residue 1, so use this parameter
383 if the PDB file does not start at residue 1)
384 @param cacenters (boolean) if True generates resolution=1 beads
385 centered on C-alpha atoms.
386 @param show (boolean) print out the molecular hierarchy at the end.
387 @param isnucleicacid (boolean) use True if you're reading a PDB
389 @param readnonwateratoms (boolean) if True fixes some pathological PDB.
390 @param read_ca_cb_only (boolean) if True, only reads CA/CB
393 self.representation_is_modified =
True
396 color = self.color_dict[name]
397 protein_h = self.hier_dict[name]
407 if isinstance(chain, str):
415 t.get_children()[0].get_children()[0]).
get_index()
417 t.get_children()[0].get_children()[-1]).
get_index()
418 c =
IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
422 IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
429 if not resrange
is None:
430 if resrange[0] > start:
432 if resrange[1] < end:
435 if not isnucleicacid:
439 residue_indexes=list(range(
442 atom_type=IMP.atom.AT_CA)
449 residue_indexes=list(range(
452 atom_type=IMP.atom.AT_P)
454 ps = sel.get_selected_particles()
456 raise ValueError(
"%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
457 % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
465 if par.get_parent() != c0:
466 par.get_parent().remove_child(par)
468 start = start + offset
471 self.elements[name].append(
472 (start, end, pdbname.split(
"/")[-1] +
":" + chain,
"pdb"))
475 resolutions, isnucleicacid, c0, protein_h,
"pdb", color)
476 self._copy_pdb_provenance(t, hiers[0], offset)
478 for p, state
in self._protocol_output:
479 p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
493 for r
in residues.keys():
495 self.model.remove_particle(c0)
501 def _copy_pdb_provenance(self, pdb, h, offset):
502 """Copy the provenance information from the PDB to our hierarchy"""
503 c =
IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
510 def add_component_ideal_helix(
518 self.representation_is_modified =
True
519 from math
import pi, cos, sin
521 protein_h = self.hier_dict[name]
524 color = self.color_dict[name]
528 self.elements[name].append((start, end,
" ",
"helix"))
530 for n, res
in enumerate(range(start, end + 1)):
531 if name
in self.sequence_dict:
533 rtstr = self.onetothree[
534 self.sequence_dict[name][res-1]]
554 x = 2.3 * cos(n * 2 * pi / 3.6)
555 y = 2.3 * sin(n * 2 * pi / 3.6)
556 z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
565 resolutions,
False, c0, protein_h,
"helix", color)
574 """ add beads to the representation
575 @param name the component name
576 @param ds a list of tuples corresponding to the residue ranges
578 @param colors a list of colors associated to the beads
579 @param incoord the coordinate tuple corresponding to the position
584 self.representation_is_modified =
True
586 protein_h = self.hier_dict[name]
589 colors = [self.color_dict[name]]
592 for n, dss
in enumerate(ds):
593 ds_frag = (dss[0], dss[1])
594 self.elements[name].append((dss[0], dss[1],
" ",
"bead"))
596 if ds_frag[0] == ds_frag[1]:
598 if name
in self.sequence_dict:
600 rtstr = self.onetothree[
601 self.sequence_dict[name][ds_frag[0]-1]]
608 h.set_name(name +
'_%i_bead' % (ds_frag[0]))
609 prt.set_name(name +
'_%i_bead' % (ds_frag[0]))
613 h.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
614 prt.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
615 h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
616 resolution = len(h.get_residue_indexes())
617 if "Beads" not in self.hier_representation[name]:
619 root.set_name(
"Beads")
620 self.hier_representation[name][
"Beads"] = root
621 protein_h.add_child(root)
622 self.hier_representation[name][
"Beads"].add_child(h)
624 for kk
in range(ds_frag[0], ds_frag[1] + 1):
625 self.hier_db.add_particles(name, kk, resolution, [h])
648 ptem.set_radius(radius)
651 radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
653 ptem.set_radius(radius)
655 if not tuple(incoord)
is None:
656 ptem.set_coordinates(incoord)
661 self.floppy_bodies.append(prt)
665 for p, state
in self._protocol_output:
666 p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
674 Generates a string of beads with given length.
676 self.representation_is_modified =
True
684 [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
688 self, name, hierarchies=
None, selection_tuples=
None,
690 resolution=0.0, num_components=10,
691 inputfile=
None, outputfile=
None,
694 covariance_type=
'full', voxel_size=1.0,
696 sampled_points=1000000, num_iter=100,
698 multiply_by_total_mass=
True,
700 intermediate_map_fn=
None,
701 density_ps_to_copy=
None,
702 use_precomputed_gaussians=
False):
704 Sets up a Gaussian Mixture Model for this component.
705 Can specify input GMM file or it will be computed.
706 @param name component name
707 @param hierarchies set up GMM for some hierarchies
708 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
709 @param particles set up GMM for particles directly
710 @param resolution usual PMI resolution for selecting particles from the hierarchies
711 @param inputfile read the GMM from this file
712 @param outputfile fit and write the GMM to this file (text)
713 @param outputmap after fitting, create GMM density file (mrc)
714 @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
715 @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
716 @param voxel_size for creating the intermediate density map and output map.
717 lower number increases accuracy but also rasterizing time grows
718 @param out_hier_name name of the output density hierarchy
719 @param sampled_points number of points to sample. more will increase accuracy and fitting time
720 @param num_iter num GMM iterations. more will increase accuracy and fitting time
721 @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
722 @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
723 @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
724 @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
725 @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
733 self.representation_is_modified =
True
735 protein_h = self.hier_dict[name]
736 if "Densities" not in self.hier_representation[name]:
738 root.set_name(
"Densities")
739 self.hier_representation[name][
"Densities"] = root
740 protein_h.add_child(root)
743 fragment_particles = []
744 if not particles
is None:
745 fragment_particles += particles
746 if not hierarchies
is None:
748 self, resolution=resolution,
749 hierarchies=hierarchies)
750 if not selection_tuples
is None:
751 for st
in selection_tuples:
752 fragment_particles += IMP.pmi.tools.select_by_tuple(
753 self, tupleselection=st,
754 resolution=resolution,
755 name_is_ambiguous=
False)
758 density_particles = []
761 inputfile, density_particles,
762 self.model, transform)
763 elif density_ps_to_copy:
764 for ip
in density_ps_to_copy:
770 density_particles.append(p)
771 elif use_precomputed_gaussians:
772 if len(fragment_particles) == 0:
773 print(
"add_component_density: no particle was selected")
775 for p
in fragment_particles:
778 raise Exception(
"The particles you selected must be Fragments and XYZs")
779 nres=len(
IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
780 pos=
IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
785 raise Exception(
"We haven't created a bead file for",nres,
"residues")
789 self.model, transform)
792 if len(fragment_particles) == 0:
793 print(
"add_component_density: no particle was selected")
796 density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
805 multiply_by_total_mass,
811 s0.set_name(out_hier_name)
812 self.hier_representation[name][
"Densities"].add_child(s0)
814 for nps, p
in enumerate(density_particles):
816 p.set_name(s0.get_name() +
'_gaussian_%i' % nps)
819 def get_component_density(self, name):
820 return self.hier_representation[name][
"Densities"]
823 selection_tuples=
None,
828 '''Decorates all specified particles as Gaussians directly.
829 @param name component name
830 @param hierarchies set up GMM for some hierarchies
831 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
832 @param particles set up GMM for particles directly
833 @param resolution usual PMI resolution for selecting particles from the hierarchies
839 from math
import sqrt
840 self.representation_is_modified =
True
842 if particles
is None:
843 fragment_particles = []
845 fragment_particles = particles
847 if not hierarchies
is None:
849 self, resolution=resolution,
850 hierarchies=hierarchies)
852 if not selection_tuples
is None:
853 for st
in selection_tuples:
854 fragment_particles += IMP.pmi.tools.select_by_tuple(
855 self, tupleselection=st,
856 resolution=resolution,
857 name_is_ambiguous=
False)
859 if len(fragment_particles) == 0:
860 print(
"add all atom densities: no particle was selected")
864 print(
'setting up all atom gaussians num_particles',len(fragment_particles))
865 for n,p
in enumerate(fragment_particles):
873 print(
'setting up particle',p.get_name(),
" as individual gaussian particle")
875 if not output_map
is None:
876 print(
'writing map to', output_map)
884 Make a copy of a hierarchy and append it to a component.
887 self.representation_is_modified =
True
888 protein_h = self.hier_dict[name]
889 hierclone = IMP.atom.create_clone(hierarchy)
890 hierclone.set_name(hierclone.get_name() +
"_clone")
891 protein_h.add_child(hierclone)
892 outhier.append(hierclone)
898 for n, pmain
in enumerate(psmain):
904 self.hier_db.add_particles(
920 def dump_particle_descriptors(self):
926 particles_attributes={}
927 floppy_body_attributes={}
933 hparent=h.get_parent()
934 while hparent != hroot:
935 hparent=h.get_parent()
936 name+=
"|"+hparent.get_name()
938 particles_attributes[name]={
"COORDINATES":numpy.array(
IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
944 rigid_body_attributes={}
945 for rb
in self.rigid_bodies:
947 rf=rb.get_reference_frame()
948 t=rf.get_transformation_to()
949 trans=t.get_translation()
951 rigid_body_attributes[name]={
"TRANSLATION":numpy.array(trans),
952 "ROTATION":numpy.array(rot.get_quaternion()),
953 "COORDINATES_NONRIGID_MEMBER":{},
954 "COORDINATES_RIGID_MEMBER":{}}
955 for mi
in rb.get_member_indexes():
956 rm=self.model.get_particle(mi)
958 name_part=rm.get_name()
960 rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
962 name_part=rm.get_name()
964 rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
968 pickle.dump(particles_attributes,
969 open(
"particles_attributes.pkl",
"wb"))
970 pickle.dump(rigid_body_attributes,
971 open(
"rigid_body_attributes.pkl",
"wb"))
975 def load_particle_descriptors(self):
981 particles_attributes = pickle.load(open(
"particles_attributes.pkl",
983 rigid_body_attributes = pickle.load(open(
"rigid_body_attributes.pkl",
993 hparent=h.get_parent()
994 while hparent != hroot:
995 hparent=h.get_parent()
996 name+=
"|"+hparent.get_name()
1000 xyzr.set_coordinates(particles_attributes[name][
"COORDINATES"])
1006 for rb
in self.rigid_bodies:
1008 trans=rigid_body_attributes[name][
"TRANSLATION"]
1009 rot=rigid_body_attributes[name][
"ROTATION"]
1012 rb.set_reference_frame(rf)
1013 coor_nrm_ref=rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"]
1014 coor_rm_ref_dict=rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"]
1017 for mi
in rb.get_member_indexes():
1018 rm=self.model.get_particle(mi)
1020 name_part=rm.get_name()
1021 xyz=coor_nrm_ref[name_part]
1023 self.model.set_attribute(fk, rm,xyz[n])
1025 name_part=rm.get_name()
1027 coor_rm_model.append(
IMP.core.XYZ(rm).get_coordinates())
1028 if len(coor_rm_model)==0:
continue
1034 def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1035 """Print a warning if particle names in RMF and model don't match"""
1036 def match_any_suffix():
1038 suffixes = [
"pdb",
"bead_floppy_body_rigid_body_member_floppy_body",
1039 "bead_floppy_body_rigid_body_member",
1042 if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1044 if rmfname != reprname
and not match_any_suffix():
1046 "set_coordinates_from_rmf: rmf particle and "
1047 "representation particle names don't match %s %s"
1052 rmf_component_name=
None,
1053 check_number_particles=
True,
1054 representation_name_to_rmf_name_map=
None,
1056 skip_gaussian_in_rmf=
False,
1057 skip_gaussian_in_representation=
False,
1059 force_rigid_update=
False):
1060 '''Read and replace coordinates from an RMF file.
1061 Replace the coordinates of particles with the same name.
1062 It assumes that the RMF and the representation have the particles
1064 @param component_name Component name
1065 @param rmf_component_name Name of the component in the RMF file
1066 (if not specified, use `component_name`)
1067 @param representation_name_to_rmf_name_map a dictionary that map
1068 the original rmf particle name to the recipient particle component name
1069 @param save_file: save a file with the names of particles of the component
1070 @param force_rigid_update: update the coordinates of rigid bodies
1071 (normally this should be called before rigid bodies are set up)
1076 prots = IMP.pmi.analysis.get_hiers_from_rmf(
1082 raise ValueError(
"cannot read hierarchy from rmf")
1086 if force_rigid_update:
1097 p, self.hier_dict.keys())
1098 if (protname
is None)
and (rmf_component_name
is not None):
1100 p, rmf_component_name)
1101 if (skip_gaussian_in_rmf):
1104 if (rmf_component_name
is not None)
and (protname == rmf_component_name):
1106 elif (rmf_component_name
is None)
and (protname == component_name):
1110 if (skip_gaussian_in_representation):
1121 reprnames=[p.get_name()
for p
in psrepr]
1122 rmfnames=[p.get_name()
for p
in psrmf]
1125 fl=open(component_name+
".txt",
"w")
1126 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1129 if check_number_particles
and not representation_name_to_rmf_name_map:
1130 if len(psrmf) != len(psrepr):
1131 fl=open(component_name+
".txt",
"w")
1132 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1133 raise ValueError(
"%s cannot proceed the rmf and the representation don't have the same number of particles; "
1134 "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1137 if not representation_name_to_rmf_name_map:
1138 for n, prmf
in enumerate(psrmf):
1140 prmfname = prmf.get_name()
1141 preprname = psrepr[n].get_name()
1142 if force_rigid_update:
1148 raise ValueError(
"component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1150 self._compare_rmf_repr_names(prmfname, preprname,
1159 rbm.set_internal_coordinates(xyz)
1161 rb = rbm.get_rigid_body()
1163 raise ValueError(
"Cannot handle nested "
1165 rb.set_reference_frame_lazy(tr)
1168 "set_coordinates_from_rmf: particles are not XYZ "
1177 g=gprmf.get_gaussian()
1178 grepr.set_gaussian(g)
1181 repr_name_particle_map={}
1182 rmf_name_particle_map={}
1184 rmf_name_particle_map[p.get_name()]=p
1188 for prepr
in psrepr:
1190 prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1193 "set_coordinates_from_rmf: missing particle name in "
1194 "representation_name_to_rmf_name_map, skipping...",
1204 If the root hierarchy does not exist, construct it.
1206 if "Res:" + str(int(resolution))
not in self.hier_representation[name]:
1208 root.set_name(name +
"_Res:" + str(int(resolution)))
1209 self.hier_representation[name][
1210 "Res:" + str(int(resolution))] = root
1211 protein_h.add_child(root)
1214 input_hierarchy, protein_h, type, color):
1216 Generate all needed coarse grained layers.
1218 @param name name of the protein
1219 @param resolutions list of resolutions
1220 @param protein_h root hierarchy
1221 @param input_hierarchy hierarchy to coarse grain
1222 @param type a string, typically "pdb" or "helix"
1226 if (1
in resolutions)
or (0
in resolutions):
1229 if 1
in resolutions:
1232 s1.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1234 self.hier_representation[name][
"Res:1"].add_child(s1)
1236 if 0
in resolutions:
1239 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1241 self.hier_representation[name][
"Res:0"].add_child(s0)
1244 if not isnucleicacid:
1247 atom_type=IMP.atom.AT_CA)
1251 atom_type=IMP.atom.AT_P)
1253 for p
in sel.get_selected_particles():
1255 if 0
in resolutions:
1257 resclone0 = IMP.atom.create_clone(resobject)
1259 s0.add_child(resclone0)
1260 self.hier_db.add_particles(
1264 resclone0.get_children())
1266 chil = resclone0.get_children()
1275 if 1
in resolutions:
1277 resclone1 = IMP.atom.create_clone_one(resobject)
1279 s1.add_child(resclone1)
1280 self.hier_db.add_particles(name, resindex, 1, [resclone1])
1284 prt = resclone1.get_particle()
1285 prt.set_name(
'%s_%i_%s' % (name, resindex, type))
1307 for r
in resolutions:
1308 if r != 0
and r != 1:
1314 chil = s.get_children()
1316 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1321 self.hier_representation[name][
1322 "Res:" + str(int(r))].add_child(s0)
1331 prt.set_name(
'%s_%i_%s' % (name, first, type))
1333 prt.set_name(
'%s_%i_%i_%s' % (name, first, last, type))
1335 self.hier_db.add_particles(name, kk, r, [prt])
1356 Get the hierarchies at the given resolution.
1358 The map between resolution and hierarchies is cached to accelerate
1359 the selection algorithm. The cache is invalidated when the
1360 representation was changed by any add_component_xxx.
1363 if self.representation_is_modified:
1364 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1366 self.hier_resolution[resolution] = rhbr
1367 self.representation_is_modified =
False
1370 if resolution
in self.hier_resolution:
1371 return self.hier_resolution[resolution]
1373 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1375 self.hier_resolution[resolution] = rhbr
1379 self, max_translation=300., max_rotation=2.0 * pi,
1380 avoidcollision=
True, cutoff=10.0, niterations=100,
1382 excluded_rigid_bodies=
None,
1383 ignore_initial_coordinates=
False,
1384 hierarchies_excluded_from_collision=
None):
1386 Shuffle configuration; used to restart the optimization.
1388 The configuration of the system is initialized by placing each
1389 rigid body and each bead randomly in a box with a side of
1390 `max_translation` angstroms, and far enough from each other to
1391 prevent any steric clashes. The rigid bodies are also randomly rotated.
1393 @param avoidcollision check if the particle/rigid body was
1394 placed close to another particle; uses the optional
1395 arguments cutoff and niterations
1396 @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1399 if excluded_rigid_bodies
is None:
1400 excluded_rigid_bodies = []
1401 if hierarchies_excluded_from_collision
is None:
1402 hierarchies_excluded_from_collision = []
1404 if len(self.rigid_bodies) == 0:
1405 print(
"shuffle_configuration: rigid bodies were not intialized")
1408 gcpf.set_distance(cutoff)
1410 hierarchies_excluded_from_collision_indexes = []
1419 if bounding_box
is not None:
1420 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1425 for h
in hierarchies_excluded_from_collision:
1429 allparticleindexes = list(
1430 set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1432 print(hierarchies_excluded_from_collision)
1433 print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1435 print(
'shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies),
'rigid bodies')
1436 for rb
in self.rigid_bodies:
1437 if rb
not in excluded_rigid_bodies:
1439 rbindexes = rb.get_member_particle_indexes()
1441 set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1442 otherparticleindexes = list(
1443 set(allparticleindexes) - set(rbindexes))
1445 if len(otherparticleindexes)
is None:
1449 while niter < niterations:
1450 if (ignore_initial_coordinates):
1456 if bounding_box
is not None:
1472 gcpf.get_close_pairs(
1474 otherparticleindexes,
1478 if (ignore_initial_coordinates):
1479 print (rb.get_name(),
IMP.core.XYZ(rb).get_coordinates())
1482 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1483 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
1484 if niter == niterations:
1485 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1489 print(
'shuffling', len(self.floppy_bodies),
'floppy bodies')
1490 for fb
in self.floppy_bodies:
1491 if (avoidcollision):
1496 otherparticleindexes = list(
1497 set(allparticleindexes) - set(fbindexes))
1498 if len(otherparticleindexes)
is None:
1516 while niter < niterations:
1517 if (ignore_initial_coordinates):
1523 if bounding_box
is not None:
1535 if (avoidcollision):
1538 gcpf.get_close_pairs(
1540 otherparticleindexes,
1544 if (ignore_initial_coordinates):
1545 print (fb.get_name(),
IMP.core.XYZ(fb).get_coordinates())
1548 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1549 print(
"shuffle_configuration: floppy body name: " + fb.get_name())
1550 if niter == niterations:
1551 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1555 def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1559 self.reference_structures[label] = (
1563 def get_all_rmsds(self):
1566 for label
in self.reference_structures:
1568 current_coordinates = [
IMP.core.XYZ(p).get_coordinates()
1569 for p
in self.reference_structures[label][1]]
1570 reference_coordinates = self.reference_structures[label][0]
1571 if len(reference_coordinates) != len(current_coordinates):
1572 print(
"calculate_all_rmsds: reference and actual coordinates are not the same")
1575 current_coordinates,
1576 reference_coordinates)
1578 reference_coordinates,
1579 current_coordinates)
1583 reference_coordinates,
1584 current_coordinates)
1585 rmsds[label +
"_GlobalRMSD"] = rmsd_global
1586 rmsds[label +
"_RelativeDRMS"] = rmsd_relative
1589 def setup_component_geometry(self, name, color=None, resolution=1.0):
1591 color = self.color_dict[name]
1595 self.hier_geometry_pairs[name] = []
1596 protein_h = self.hier_dict[name]
1598 pbr = IMP.pmi.tools.sort_by_residues(pbr)
1600 for n
in range(len(pbr) - 1):
1601 self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1604 self, name, resolution=10, scale=1.0):
1606 Generate restraints between contiguous fragments in the hierarchy.
1607 The linkers are generated at resolution 10 by default.
1614 protein_h = self.hier_dict[name]
1616 frs = self.hier_db.get_preroot_fragments_by_resolution(
1622 start = fr.get_children()[0]
1627 end = fr.get_children()[-1]
1633 SortedSegments.append((start, end, startres))
1634 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1637 for x
in range(len(SortedSegments) - 1):
1638 last = SortedSegments[x][1]
1639 first = SortedSegments[x + 1][0]
1646 residuegap = firstresn - lastresn - 1
1647 if self.disorderedlength
and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1650 optdist = sqrt(5 / 3) * 1.93 * \
1651 (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1653 if self.upperharmonic:
1659 optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1660 if self.upperharmonic:
1666 pt0 = last.get_particle()
1667 pt1 = first.get_particle()
1670 print(
"Adding sequence connectivity restraint between", pt0.get_name(),
" and ", pt1.get_name(),
'of distance', optdist)
1671 sortedsegments_cr.add_restraint(r)
1672 self.linker_restraints_dict[
1673 "LinkerRestraint-" + pt0.get_name() +
"-" + pt1.get_name()] = r
1674 self.connected_intra_pairs.append((pt0, pt1))
1675 self.connected_intra_pairs.append((pt1, pt0))
1677 self.linker_restraints.add_restraint(sortedsegments_cr)
1678 self.linker_restraints.add_restraint(unmodeledregions_cr)
1680 self.sortedsegments_cr_dict[name] = sortedsegments_cr
1681 self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1683 def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1685 pts = IMP.pmi.tools.ParticleToSampleList()
1686 for n, fb
in enumerate(self.floppy_bodies):
1687 pts.add_particle(fb,
"Floppy_Bodies", 1.0,
"Floppy_Body_" + str(n))
1688 if len(pts.get_particles_to_sample()) > 0:
1690 print(
"optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1693 print(
"optimize_floppy_bodies: no particle to optimize")
1696 nSymmetry=
None, skip_gaussian_in_clones=
False):
1698 The copies must not contain rigid bodies.
1699 The symmetry restraints are applied at each leaf.
1703 self.representation_is_modified =
True
1704 ncopies = len(copies) + 1
1707 for k
in range(len(copies)):
1708 if (nSymmetry
is None):
1709 rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1712 rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1714 rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1722 for n, p
in enumerate(main_hiers):
1723 if (skip_gaussian_in_clones):
1732 self.model.add_score_state(c)
1733 print(
"Completed setting " + str(maincopy) +
" as a reference for "
1734 + str(copies[k]) +
" by rotating it by "
1735 + str(rotation_angle / 2.0 / pi * 360)
1736 +
" degrees around the " + str(rotational_axis) +
" axis.")
1739 def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1742 self.representation_is_modified =
True
1744 mainparticles = particles_reference
1746 t=initial_transformation
1748 p.set_name(
"RigidBody_Symmetry")
1754 copyparticles = particles_copy
1758 for n, p
in enumerate(mainparticles):
1760 pc = copyparticles[n]
1762 mainpurged.append(p)
1768 copypurged.append(pc)
1775 for n, p
in enumerate(mainpurged):
1778 print(
"setting " + p.get_name() +
" as reference for " + pc.get_name())
1781 lc.add(pc.get_index())
1784 self.model.add_score_state(c)
1787 self.rigid_bodies.append(rb)
1788 self.rigid_body_symmetries.append(rb)
1789 rb.set_name(label+
".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1792 def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1793 longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1796 self.representation_is_modified =
True
1799 protein_h = self.hier_dict[maincopy]
1802 for ilc
in range(-longitudinal_copies, longitudinal_copies + 1):
1803 for iac
in range(axial_copies):
1804 copyname = maincopy +
"_a" + str(ilc) +
"_l" + str(iac)
1805 self.create_component(copyname, 0.0)
1806 for hier
in protein_h.get_children():
1808 copyhier = self.hier_dict[copyname]
1809 outhiers.append(copyhier)
1813 2 * pi / axial_copies * (float(iac)))
1814 translation_vector = tuple(
1815 [translation_value * float(ilc) * x
for x
in axis])
1816 print(translation_vector)
1821 for n, p
in enumerate(mainparts):
1828 lc.add(pc.get_index())
1830 self.model.add_score_state(c)
1836 Load coordinates in the current representation.
1837 This should be done only if the hierarchy self.prot is identical
1838 to the one as stored in the rmf i.e. all components were added.
1843 rh = RMF.open_rmf_file_read_only(rmfname)
1851 create the representation (i.e. hierarchies) from the rmf file.
1852 it will be stored in self.prot, which will be overwritten.
1853 load the coordinates from the rmf file at frameindex.
1855 rh = RMF.open_rmf_file_read_only(rmfname)
1861 for p
in self.prot.get_children():
1862 self.create_component(p.get_name())
1863 self.hier_dict[p.get_name()] = p
1865 still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1866 save floppy bodies in self.floppy_bodies
1867 get the connectivity restraints
1872 Construct a rigid body from hierarchies (and optionally particles).
1874 @param hiers list of hierarchies to make rigid
1875 @param particles list of particles to add to the rigid body
1878 if particles
is None:
1881 rigid_parts = set(particles)
1884 print(
"set_rigid_body_from_hierarchies> setting up a new rigid body")
1892 "set_rigid_body_from_hierarchies> particle %s "
1893 "already belongs to rigid body %s"
1894 % (p.get_name(), rb.get_name()),
1898 name += hier.get_name() +
"-"
1899 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1901 if len(list(rigid_parts)) != 0:
1903 rb.set_coordinates_are_optimized(
True)
1904 rb.set_name(name +
"rigid_body")
1905 self.rigid_bodies.append(rb)
1909 print(
"set_rigid_body_from_hierarchies> rigid body could not be setup")
1913 Construct a rigid body from a list of subunits.
1915 Each subunit is a tuple that identifies the residue ranges and the
1916 component name (as used in create_component()).
1918 subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1920 [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1922 example: ["prot1","prot2",("prot3",(1,10))]
1924 sometimes, we know about structure of an interaction
1925 and here we make such PPIs rigid
1930 if isinstance(s, tuple)
and len(s) == 2:
1934 residue_indexes=list(range(s[1][0],
1936 if len(sel.get_selected_particles()) == 0:
1937 print(
"set_rigid_bodies: selected particle does not exist")
1938 for p
in sel.get_selected_particles():
1943 "set_rigid_body_from_hierarchies> particle %s "
1944 "already belongs to rigid body %s"
1945 % (p.get_name(), rb.get_name()),
1950 elif isinstance(s, str):
1952 if len(sel.get_selected_particles()) == 0:
1953 print(
"set_rigid_bodies: selected particle does not exist")
1954 for p
in sel.get_selected_particles():
1959 "set_rigid_body_from_hierarchies> particle %s "
1960 "already belongs to rigid body %s"
1961 % (p.get_name(), rb.get_name()),
1967 rb.set_coordinates_are_optimized(
True)
1968 rb.set_name(
''.join(str(subunits)) +
"_rigid_body")
1969 self.rigid_bodies.append(rb)
1972 def set_super_rigid_body_from_hierarchies(
1978 super_rigid_xyzs = set()
1979 super_rigid_rbs = set()
1981 print(
"set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
1988 super_rigid_rbs.add(rb)
1990 super_rigid_xyzs.add(p)
1991 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1992 if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
1995 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
1998 self.super_rigid_bodies.append(
1999 (super_rigid_xyzs, super_rigid_rbs, axis))
2001 def fix_rigid_bodies(self, rigid_bodies):
2002 self.fixed_rigid_bodies += rigid_bodies
2006 self, hiers, min_length=
None, max_length=
None, axis=
None):
2008 Make a chain of super rigid bodies from a list of hierarchies.
2010 Takes a linear list of hierarchies (they are supposed to be
2011 sequence-contiguous) and produces a chain of super rigid bodies
2012 with given length range, specified by min_length and max_length.
2015 hiers = IMP.pmi.tools.flatten_list(hiers)
2019 self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2021 def set_super_rigid_bodies(self, subunits, coords=None):
2022 super_rigid_xyzs = set()
2023 super_rigid_rbs = set()
2026 if isinstance(s, tuple)
and len(s) == 3:
2030 residue_indexes=list(range(s[0],
2032 if len(sel.get_selected_particles()) == 0:
2033 print(
"set_rigid_bodies: selected particle does not exist")
2034 for p
in sel.get_selected_particles():
2037 super_rigid_rbs.add(rb)
2039 super_rigid_xyzs.add(p)
2040 elif isinstance(s, str):
2042 if len(sel.get_selected_particles()) == 0:
2043 print(
"set_rigid_bodies: selected particle does not exist")
2044 for p
in sel.get_selected_particles():
2048 super_rigid_rbs.add(rb)
2050 super_rigid_xyzs.add(p)
2051 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2055 Remove leaves of hierarchies from the floppy bodies list based
2056 on the component name
2059 ps=[h.get_particle()
for h
in hs]
2060 for p
in self.floppy_bodies:
2062 if p
in ps: self.floppy_bodies.remove(p)
2063 if p
in hs: self.floppy_bodies.remove(p)
2070 Remove leaves of hierarchies from the floppy bodies list.
2072 Given a list of hierarchies, get the leaves and remove the
2073 corresponding particles from the floppy bodies list. We need this
2074 function because sometimes
2075 we want to constrain the floppy bodies in a rigid body - for instance
2076 when you want to associate a bead with a density particle.
2078 for h
in hierarchies:
2081 if p
in self.floppy_bodies:
2082 print(
"remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2083 self.floppy_bodies.remove(p)
2086 def set_floppy_bodies(self):
2087 for p
in self.floppy_bodies:
2089 p.set_name(name +
"_floppy_body")
2091 print(
"I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2094 rb.set_is_rigid_member(p.get_index(),
False)
2097 rb.set_is_rigid_member(p.get_particle_index(),
False)
2098 p.set_name(p.get_name() +
"_rigid_body_member")
2100 def set_floppy_bodies_from_hierarchies(self, hiers):
2105 self.floppy_bodies.append(p)
2112 selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2113 @return the particles
2116 print(selection_tuples)
2117 for s
in selection_tuples:
2118 ps = IMP.pmi.tools.select_by_tuple(
2119 representation=self, tupleselection=tuple(s),
2120 resolution=
None, name_is_ambiguous=
False)
2124 def get_connected_intra_pairs(self):
2125 return self.connected_intra_pairs
2127 def set_rigid_bodies_max_trans(self, maxtrans):
2128 self.maxtrans_rb = maxtrans
2130 def set_rigid_bodies_max_rot(self, maxrot):
2131 self.maxrot_rb = maxrot
2133 def set_super_rigid_bodies_max_trans(self, maxtrans):
2134 self.maxtrans_srb = maxtrans
2136 def set_super_rigid_bodies_max_rot(self, maxrot):
2137 self.maxrot_srb = maxrot
2139 def set_floppy_bodies_max_trans(self, maxtrans):
2140 self.maxtrans_fb = maxtrans
2144 Fix rigid bodies in their actual position.
2145 The get_particles_to_sample() function will return
2146 just the floppy bodies (if they are not fixed).
2148 self.rigidbodiesarefixed = rigidbodiesarefixed
2152 Fix floppy bodies in their actual position.
2153 The get_particles_to_sample() function will return
2154 just the rigid bodies (if they are not fixed).
2156 self.floppybodiesarefixed=floppybodiesarefixed
2158 def draw_hierarchy_graph(self):
2160 print(
"Drawing hierarchy graph for " + c.get_name())
2161 IMP.pmi.output.get_graph_from_hierarchy(c)
2163 def get_geometries(self):
2166 for name
in self.hier_geometry_pairs:
2167 for pt
in self.hier_geometry_pairs[name]:
2181 def setup_bonds(self):
2184 for name
in self.hier_geometry_pairs:
2185 for pt
in self.hier_geometry_pairs[name]:
2200 def show_component_table(self, name):
2201 if name
in self.sequence_dict:
2202 lastresn = len(self.sequence_dict[name])
2205 residues = self.hier_db.get_residue_numbers(name)
2206 firstresn = min(residues)
2207 lastresn = max(residues)
2209 for nres
in range(firstresn, lastresn + 1):
2211 resolutions = self.hier_db.get_residue_resolutions(name, nres)
2213 for r
in resolutions:
2214 ps += self.hier_db.get_particles(name, nres, r)
2215 print(
"%20s %7s" % (name, nres),
" ".join([
"%20s %7s" % (str(p.get_name()),
2218 print(
"%20s %20s" % (name, nres),
"**** not represented ****")
2220 def draw_hierarchy_composition(self):
2222 ks = list(self.elements.keys())
2226 for k
in self.elements:
2227 for l
in self.elements[k]:
2232 self.draw_component_composition(k, max)
2234 def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2235 from matplotlib
import pyplot
2236 import matplotlib
as mpl
2238 tmplist = sorted(self.elements[k], key=itemgetter(0))
2240 endres = tmplist[-1][1]
2242 print(
"draw_component_composition: missing information for component %s" % name)
2244 fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2245 ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2250 norm = mpl.colors.Normalize(vmin=5, vmax=10)
2254 for n, l
in enumerate(tmplist):
2258 if bounds[-1] != l[0]:
2259 colors.append(
"white")
2262 colors.append(
"#99CCFF")
2264 colors.append(
"#FFFF99")
2266 colors.append(
"#33CCCC")
2268 bounds.append(l[1] + 1)
2271 colors.append(
"#99CCFF")
2273 colors.append(
"#FFFF99")
2275 colors.append(
"#33CCCC")
2277 bounds.append(l[1] + 1)
2279 if bounds[-1] - 1 == l[0]:
2283 colors.append(
"white")
2286 bounds.append(bounds[-1])
2287 colors.append(
"white")
2288 cmap = mpl.colors.ListedColormap(colors)
2289 cmap.set_over(
'0.25')
2290 cmap.set_under(
'0.75')
2292 norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2293 cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2299 spacing=
'proportional',
2300 orientation=
'horizontal')
2309 mid = 1.0 / endres * float(l[0])
2314 l[2], xy=(mid, 1), xycoords=
'axes fraction',
2315 xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords=
'axes fraction',
2316 arrowprops=dict(arrowstyle=
"->",
2317 connectionstyle=
"angle,angleA=0,angleB=90,rad=10"),
2319 extra_artists.append(t)
2322 title = ax.text(-0.005, 0.5, k, ha=
"right", va=
"center", rotation=90,
2325 extra_artists.append(title)
2327 labels = len(bounds) * [
" "]
2328 ax.set_xticklabels(labels)
2329 mid = 1.0 / endres * float(bounds[0])
2330 t = ax.annotate(bounds[0], xy=(mid, 0), xycoords=
'axes fraction',
2331 xytext=(mid - 0.01, -0.5), textcoords=
'axes fraction',)
2332 extra_artists.append(t)
2333 offsets = [0, -0.5, -1.0]
2335 for n
in range(1, len(bounds)):
2336 if bounds[n] == bounds[n - 1]:
2338 mid = 1.0 / endres * float(bounds[n])
2339 if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2341 offset = offsets[nclashes % 3]
2347 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2348 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction')
2351 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2352 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction', arrowprops=dict(arrowstyle=
"-"))
2353 extra_artists.append(t)
2355 cb2.add_lines(bounds, [
"black"] * len(bounds), [1] * len(bounds))
2359 k +
"structure.pdf",
2362 bbox_extra_artists=(extra_artists),
2363 bbox_inches=
'tight')
2366 def draw_coordinates_projection(self):
2367 import matplotlib.pyplot
as pp
2370 for name
in self.hier_geometry_pairs:
2371 for pt
in self.hier_geometry_pairs[name]:
2381 xpairs.append([x1, x2])
2382 ypairs.append([y1, y2])
2385 for xends, yends
in zip(xpairs, ypairs):
2390 pp.plot(xlist, ylist,
'b-', alpha=0.1)
2393 def get_prot_name_from_particle(self, particle):
2394 names = self.get_component_names()
2395 particle0 = particle
2397 while not name
in names:
2400 particle0 = h.get_particle()
2404 def get_particles_to_sample(self):
2414 if not self.rigidbodiesarefixed:
2415 for rb
in self.rigid_bodies:
2420 if rb
not in self.fixed_rigid_bodies:
2423 if not self.floppybodiesarefixed:
2424 for fb
in self.floppy_bodies:
2431 for srb
in self.super_rigid_bodies:
2434 rigid_bodies = list(srb[1])
2435 filtered_rigid_bodies = []
2436 for rb
in rigid_bodies:
2437 if rb
not in self.fixed_rigid_bodies:
2438 filtered_rigid_bodies.append(rb)
2439 srbtmp.append((srb[0], filtered_rigid_bodies))
2441 self.rigid_bodies = rbtmp
2442 self.floppy_bodies = fbtmp
2443 self.super_rigid_bodies = srbtmp
2445 ps[
"Rigid_Bodies_SimplifiedModel"] = (
2449 ps[
"Floppy_Bodies_SimplifiedModel"] = (
2452 ps[
"SR_Bodies_SimplifiedModel"] = (
2453 self.super_rigid_bodies,
2458 def set_output_level(self, level):
2459 self.output_level = level
2461 def _evaluate(self, deriv):
2462 """Evaluate the total score of all added restraints"""
2464 return r.evaluate(deriv)
2466 def get_output(self):
2470 output[
"SimplifiedModel_Total_Score_" +
2471 self.label] = str(self._evaluate(
False))
2472 output[
"SimplifiedModel_Linker_Score_" +
2473 self.label] = str(self.linker_restraints.unprotected_evaluate(
None))
2474 for name
in self.sortedsegments_cr_dict:
2475 partialscore = self.sortedsegments_cr_dict[name].evaluate(
False)
2476 score += partialscore
2478 "SimplifiedModel_Link_SortedSegments_" +
2483 partialscore = self.unmodeledregions_cr_dict[name].evaluate(
False)
2484 score += partialscore
2486 "SimplifiedModel_Link_UnmodeledRegions_" +
2491 for rb
in self.rigid_body_symmetries:
2493 output[name +
"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2494 for name
in self.linker_restraints_dict:
2499 self.linker_restraints_dict[
2500 name].unprotected_evaluate(
2503 if len(self.reference_structures.keys()) != 0:
2504 rmsds = self.get_all_rmsds()
2507 "SimplifiedModel_" +
2510 self.label] = rmsds[
2513 if self.output_level ==
"high":
2516 output[
"Coordinates_" +
2517 p.get_name() +
"_" + self.label] = str(d)
2519 output[
"_TotalScore"] = str(score)
2522 def get_test_output(self):
2524 output = self.get_output()
2525 for n, p
in enumerate(self.get_particles_to_sample()):
2526 output[
"Particle_to_sample_" + str(n)] = str(p)
2528 output[
"Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2529 output[
"Number_of_floppy_bodies"] = len(self.floppy_bodies)
2530 output[
"Number_of_rigid_bodies"] = len(self.rigid_bodies)
2531 output[
"Number_of_super_bodies"] = len(self.super_rigid_bodies)
2532 output[
"Selection_resolution_1"] = len(
2534 output[
"Selection_resolution_5"] = len(
2536 output[
"Selection_resolution_7"] = len(
2538 output[
"Selection_resolution_10"] = len(
2540 output[
"Selection_resolution_100"] = len(
2543 output[
"Selection_resolution=1"] = len(
2545 output[
"Selection_resolution=1,resid=10"] = len(
2547 for resolution
in self.hier_resolution:
2548 output[
"Hier_resolution_dictionary" +
2549 str(resolution)] = len(self.hier_resolution[resolution])
2550 for name
in self.hier_dict:
2552 "Selection_resolution=1,resid=10,name=" +
2560 "Selection_resolution=1,resid=10,name=" +
2562 ",ambiguous"] = len(
2567 name_is_ambiguous=
True,
2570 "Selection_resolution=1,resid=10,name=" +
2572 ",ambiguous"] = len(
2577 name_is_ambiguous=
True,
2580 "Selection_resolution=1,resrange=(10,20),name=" +
2589 "Selection_resolution=1,resrange=(10,20),name=" +
2591 ",ambiguous"] = len(
2596 name_is_ambiguous=
True,
2600 "Selection_resolution=10,resrange=(10,20),name=" +
2609 "Selection_resolution=10,resrange=(10,20),name=" +
2611 ",ambiguous"] = len(
2616 name_is_ambiguous=
True,
2620 "Selection_resolution=100,resrange=(10,20),name=" +
2629 "Selection_resolution=100,resrange=(10,20),name=" +
2631 ",ambiguous"] = len(
2636 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.
Warning related to handling of structures.
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.
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)
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_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.