3 """@namespace IMP.pmi1.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."""
39 class Representation(object):
43 Set up the representation of all proteins and nucleic acid macromolecules.
45 Create the molecular hierarchies, representation,
46 sequence connectivity for the various involved proteins and
47 nucleic acid macromolecules:
49 Create a protein, DNA or RNA, represent it as a set of connected balls of appropriate
50 radii and number of residues, PDB at given resolution(s), or ideal helices.
52 How to use the SimplifiedModel class (typical use):
54 see test/test_hierarchy_contruction.py
58 1) Create a chain of helices and flexible parts
60 c_1_119 =self.add_component_necklace("prot1",1,119,20)
61 c_120_131 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(120,131))
62 c_132_138 =self.add_component_beads("prot1",[(132,138)])
63 c_139_156 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(139,156))
64 c_157_174 =self.add_component_beads("prot1",[(157,174)])
65 c_175_182 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(175,182))
66 c_183_194 =self.add_component_beads("prot1",[(183,194)])
67 c_195_216 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(195,216))
68 c_217_250 =self.add_component_beads("prot1",[(217,250)])
71 self.set_rigid_body_from_hierarchies(c_120_131)
72 self.set_rigid_body_from_hierarchies(c_139_156)
73 self.set_rigid_body_from_hierarchies(c_175_182)
74 self.set_rigid_body_from_hierarchies(c_195_216)
76 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,
79 self.set_chain_of_super_rigid_bodies(clist,2,3)
81 self.set_super_rigid_bodies(["prot1"])
85 def __init__(self, model, upperharmonic=True, disorderedlength=True):
87 @param model the model
88 @param upperharmonic This flag uses either harmonic (False)
89 or upperharmonic (True) in the intra-pair
90 connectivity restraint.
91 @param disorderedlength This flag uses either disordered length
92 calculated for random coil peptides (True) or zero
93 surface-to-surface distance between beads (False)
94 as optimal distance for the sequence connectivity
98 self.state = _StateInfo()
100 self._file_dataset = {}
101 self._protocol_output = []
102 self._non_modeled_components = {}
111 self.upperharmonic = upperharmonic
112 self.disorderedlength = disorderedlength
113 self.rigid_bodies = []
114 self.fixed_rigid_bodies = []
115 self.fixed_floppy_bodies = []
116 self.floppy_bodies = []
122 self.super_rigid_bodies = []
123 self.rigid_body_symmetries = []
124 self.output_level =
"low"
127 self.maxtrans_rb = 2.0
128 self.maxrot_rb = 0.04
129 self.maxtrans_srb = 2.0
130 self.maxrot_srb = 0.2
131 self.rigidbodiesarefixed =
False
132 self.floppybodiesarefixed =
False
133 self.maxtrans_fb = 3.0
134 self.resolution = 10.0
135 self.bblenght = 100.0
139 self.representation_is_modified =
False
140 self.unmodeledregions_cr_dict = {}
141 self.sortedsegments_cr_dict = {}
143 self.connected_intra_pairs = []
146 self.sequence_dict = {}
147 self.hier_geometry_pairs = {}
153 self.hier_representation = {}
154 self.hier_resolution = {}
157 self.reference_structures = {}
160 self.linker_restraints.set_was_used(
True)
161 self.linker_restraints_dict = {}
162 self.threetoone = {
'ALA':
'A',
'ARG':
'R', 'ASN': 'N', 'ASP': 'D',
163 'CYS':
'C',
'GLU':
'E',
'GLN':
'Q',
'GLY':
'G',
164 'HIS':
'H',
'ILE':
'I',
'LEU':
'L',
'LYS':
'K',
165 'MET':
'M',
'PHE':
'F',
'PRO':
'P',
'SER':
'S',
166 'THR':
'T',
'TRP':
'W',
'TYR':
'Y',
'VAL':
'V',
'UNK':
'X'}
168 self.onetothree = dict((v, k)
for k, v
in self.threetoone.items())
178 """Associate some metadata with this modeling.
179 @param m an instance of IMP.pmi1.metadata.Metadata or a subclass.
181 self._metadata.append(m)
184 """Associate a dataset with a filename.
185 This can be used to identify how the file was produced (in many
186 cases IMP can determine this automatically from a file header or
187 other metadata, but not always). For example, a manually-produced
188 PDB file (not from the PDB database or Modeller) can be
190 @param fname filename
191 @dataset the ihm.dataset.Dataset object to associate.
193 self._file_dataset[os.path.abspath(fname)] = dataset
196 """Get the dataset associated with a filename, or None.
197 @param fname filename
198 @return an ihm.dataset.Dataset, or None.
200 return self._file_dataset.get(os.path.abspath(fname),
None)
203 """Capture details of the modeling protocol.
204 @param p an instance of IMP.pmi1.output.ProtocolOutput or a subclass.
206 state = p._add_state(self)
207 self._protocol_output.append((p, state))
208 p._each_metadata.append(self._metadata)
209 p._file_datasets.append(self._file_dataset)
210 state.model = self.model
211 state.prot = self.prot
212 protocol_output = property(
lambda self:
213 [x[0]
for x
in self._protocol_output])
215 def set_label(self, label):
218 def create_component(self, name, color=0.0):
225 protein_h.set_name(name)
226 self.hier_dict[name] = protein_h
227 self.hier_representation[name] = {}
228 self.hier_db.add_name(name)
229 self.prot.add_child(protein_h)
230 self.color_dict[name] = color
231 self.elements[name] = []
232 for p, state
in self._protocol_output:
233 p.create_component(state, name,
True)
236 """Create a transformed view of a component.
237 This does not copy the coordinates or representation of the
238 component, and the transformed component is not visible to IMP,
239 but when the model is written to a ProtocolOutput, the transformed
240 component is added. This can be used to easily add symmetry-related
241 'copies' of a component without copying the underlying
243 for p, state
in self._protocol_output:
244 p.create_transformed_component(state, name, original, transform)
247 """Create a component that isn't used in the modeling.
248 No coordinates or other structural data for this component will
249 be read or written, but a primary sequence can be assigned. This
250 is useful if the input experimental data is of a system larger
251 than that modeled. Any references to these non-modeled components
252 can then be correctly resolved later."""
253 self._non_modeled_components[name] =
None
254 self.elements[name] = []
255 for p, state
in self._protocol_output:
256 p.create_component(state, name,
False)
261 def add_component_name(self, *args, **kwargs):
262 self.create_component(*args, **kwargs)
264 def get_component_names(self):
265 return list(self.hier_dict.keys())
270 Add the primary sequence for a single component.
272 @param name Human-readable name of the component
273 @param filename Name of the FASTA file
274 @param id Identifier of the sequence in the FASTA file header
275 (if not provided, use `name` instead)
280 if id
not in record_dict:
281 raise KeyError(
"id %s not found in fasta file" % id)
282 length = len(record_dict[id])
283 self.sequence_dict[name] = str(record_dict[id])
285 if name
not in self._non_modeled_components:
286 protein_h = self.hier_dict[name]
287 protein_h.set_sequence(self.sequence_dict[name])
290 self.sequence_dict[name]=offs_str+self.sequence_dict[name]
292 self.elements[name].append((length, length,
" ",
"end"))
293 for p, state
in self._protocol_output:
294 p.add_component_sequence(state, name, self.sequence_dict[name])
296 def autobuild_model(self, name, pdbname, chain,
297 resolutions=
None, resrange=
None,
299 color=
None, pdbresrange=
None, offset=0,
300 show=
False, isnucleicacid=
False,
303 self.representation_is_modified =
True
307 color = self.color_dict[name]
309 self.color_dict[name] = color
311 if resolutions
is None:
313 print(
"autobuild_model: constructing %s from pdb %s and chain %s" % (name, pdbname, str(chain)))
323 t.get_children()[0].get_children()[0]).
get_index()
325 t.get_children()[0].get_children()[-1]).
get_index()
331 if name
in self.sequence_dict:
332 resrange = (1, len(self.sequence_dict[name]))
334 resrange = (start + offset, end + offset)
336 if resrange[1]
in (-1,
'END'):
337 resrange = (resrange[0],end)
338 start = resrange[0] - offset
339 end = resrange[1] - offset
358 for n, g
in enumerate(gaps):
362 print(
"autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
364 chain, resolutions=resolutions,
365 color=color, cacenters=
True,
366 resrange=(first, last),
367 offset=offset, isnucleicacid=isnucleicacid)
368 elif g[2] ==
"gap" and n > 0:
369 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
370 parts = self.hier_db.get_particles_at_closest_resolution(name,
375 first+offset, last+offset, missingbeadsize, incoord=xyz)
377 elif g[2] ==
"gap" and n == 0:
379 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
381 first+offset, last+offset, missingbeadsize, incoord=xyznter)
388 def autobuild_pdb_and_intervening_beads(self, *args, **kwargs):
389 r = self.autobuild_model(*args, **kwargs)
393 resrange=
None, offset=0, cacenters=
True, show=
False,
394 isnucleicacid=
False, readnonwateratoms=
False,
395 read_ca_cb_only=
False):
397 Add a component that has an associated 3D structure in a PDB file.
399 Reads the PDB, and constructs the fragments corresponding to contiguous
402 @return a list of hierarchies.
404 @param name (string) the name of the component
405 @param pdbname (string) the name of the PDB file
406 @param chain (string or integer) can be either a string (eg, "A")
407 or an integer (eg, 0 or 1) in case you want
408 to get the corresponding chain number in the PDB.
409 @param resolutions (integers) a list of integers that corresponds
410 to the resolutions that have to be generated
411 @param color (float from 0 to 1) the color applied to the
412 hierarchies generated
413 @param resrange (tuple of integers): the residue range to extract
414 from the PDB. It is a tuple (beg,end). If not specified,
415 all residues belonging to the specified chain are read.
416 @param offset (integer) specifies the residue index offset to be
417 applied when reading the PDB (the FASTA sequence is
418 assumed to start from residue 1, so use this parameter
419 if the PDB file does not start at residue 1)
420 @param cacenters (boolean) if True generates resolution=1 beads
421 centered on C-alpha atoms.
422 @param show (boolean) print out the molecular hierarchy at the end.
423 @param isnucleicacid (boolean) use True if you're reading a PDB
425 @param readnonwateratoms (boolean) if True fixes some pathological PDB.
426 @param read_ca_cb_only (boolean) if True, only reads CA/CB
429 self.representation_is_modified =
True
432 color = self.color_dict[name]
433 protein_h = self.hier_dict[name]
443 if type(chain) == str:
449 t.get_children()[0].get_children()[0]).
get_index()
451 t.get_children()[0].get_children()[-1]).
get_index()
452 c =
IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
456 IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
463 if not resrange
is None:
464 if resrange[0] > start:
466 if resrange[1] < end:
469 if not isnucleicacid:
473 residue_indexes=list(range(
476 atom_type=IMP.atom.AT_CA)
483 residue_indexes=list(range(
486 atom_type=IMP.atom.AT_P)
488 ps = sel.get_selected_particles()
490 raise ValueError(
"%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
491 % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
499 if par.get_parent() != c0:
500 par.get_parent().remove_child(par)
502 start = start + offset
505 self.elements[name].append(
506 (start, end, pdbname.split(
"/")[-1] +
":" + chain,
"pdb"))
509 resolutions, isnucleicacid, c0, protein_h,
"pdb", color)
510 self._copy_pdb_provenance(t, hiers[0], offset)
512 for p, state
in self._protocol_output:
513 p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
527 for r
in residues.keys():
529 self.model.remove_particle(c0)
535 def _copy_pdb_provenance(self, pdb, h, offset):
536 """Copy the provenance information from the PDB to our hierarchy"""
537 c =
IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
544 def add_component_ideal_helix(
552 self.representation_is_modified =
True
553 from math
import pi, cos, sin
555 protein_h = self.hier_dict[name]
558 color = self.color_dict[name]
562 self.elements[name].append((start, end,
" ",
"helix"))
564 for n, res
in enumerate(range(start, end + 1)):
565 if name
in self.sequence_dict:
567 rtstr = self.onetothree[
568 self.sequence_dict[name][res-1]]
588 x = 2.3 * cos(n * 2 * pi / 3.6)
589 y = 2.3 * sin(n * 2 * pi / 3.6)
590 z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
599 resolutions,
False, c0, protein_h,
"helix", color)
608 """ add beads to the representation
609 @param name the component name
610 @param ds a list of tuples corresponding to the residue ranges
612 @param colors a list of colors associated to the beads
613 @param incoord the coordinate tuple corresponding to the position
618 self.representation_is_modified =
True
620 protein_h = self.hier_dict[name]
623 colors = [self.color_dict[name]]
626 for n, dss
in enumerate(ds):
627 ds_frag = (dss[0], dss[1])
628 self.elements[name].append((dss[0], dss[1],
" ",
"bead"))
630 if ds_frag[0] == ds_frag[1]:
632 if name
in self.sequence_dict:
634 rtstr = self.onetothree[
635 self.sequence_dict[name][ds_frag[0]-1]]
642 h.set_name(name +
'_%i_bead' % (ds_frag[0]))
643 prt.set_name(name +
'_%i_bead' % (ds_frag[0]))
647 h.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
648 prt.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
649 h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
650 resolution = len(h.get_residue_indexes())
651 if "Beads" not in self.hier_representation[name]:
653 root.set_name(
"Beads")
654 self.hier_representation[name][
"Beads"] = root
655 protein_h.add_child(root)
656 self.hier_representation[name][
"Beads"].add_child(h)
658 for kk
in range(ds_frag[0], ds_frag[1] + 1):
659 self.hier_db.add_particles(name, kk, resolution, [h])
682 ptem.set_radius(radius)
685 radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
687 ptem.set_radius(radius)
689 if not tuple(incoord)
is None:
690 ptem.set_coordinates(incoord)
695 self.floppy_bodies.append(prt)
699 for p, state
in self._protocol_output:
700 p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
708 Generates a string of beads with given length.
710 self.representation_is_modified =
True
718 [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
722 self, name, hierarchies=
None, selection_tuples=
None,
724 resolution=0.0, num_components=10,
725 inputfile=
None, outputfile=
None,
728 covariance_type=
'full', voxel_size=1.0,
730 sampled_points=1000000, num_iter=100,
732 multiply_by_total_mass=
True,
734 intermediate_map_fn=
None,
735 density_ps_to_copy=
None,
736 use_precomputed_gaussians=
False):
738 Sets up a Gaussian Mixture Model for this component.
739 Can specify input GMM file or it will be computed.
740 @param name component name
741 @param hierarchies set up GMM for some hierarchies
742 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
743 @param particles set up GMM for particles directly
744 @param resolution usual PMI resolution for selecting particles from the hierarchies
745 @param inputfile read the GMM from this file
746 @param outputfile fit and write the GMM to this file (text)
747 @param outputmap after fitting, create GMM density file (mrc)
748 @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
749 @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
750 @param voxel_size for creating the intermediate density map and output map.
751 lower number increases accuracy but also rasterizing time grows
752 @param out_hier_name name of the output density hierarchy
753 @param sampled_points number of points to sample. more will increase accuracy and fitting time
754 @param num_iter num GMM iterations. more will increase accuracy and fitting time
755 @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
756 @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
757 @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
758 @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
759 @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
767 self.representation_is_modified =
True
769 protein_h = self.hier_dict[name]
770 if "Densities" not in self.hier_representation[name]:
772 root.set_name(
"Densities")
773 self.hier_representation[name][
"Densities"] = root
774 protein_h.add_child(root)
777 fragment_particles = []
778 if not particles
is None:
779 fragment_particles += particles
780 if not hierarchies
is None:
782 self, resolution=resolution,
783 hierarchies=hierarchies)
784 if not selection_tuples
is None:
785 for st
in selection_tuples:
786 fragment_particles += IMP.pmi1.tools.select_by_tuple(
787 self, tupleselection=st,
788 resolution=resolution,
789 name_is_ambiguous=
False)
792 density_particles = []
795 inputfile, density_particles,
796 self.model, transform)
797 elif density_ps_to_copy:
798 for ip
in density_ps_to_copy:
804 density_particles.append(p)
805 elif use_precomputed_gaussians:
806 if len(fragment_particles) == 0:
807 print(
"add_component_density: no particle was selected")
809 for p
in fragment_particles:
812 raise Exception(
"The particles you selected must be Fragments and XYZs")
813 nres=len(
IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
814 pos=
IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
819 raise Exception(
"We haven't created a bead file for",nres,
"residues")
823 self.model, transform)
826 if len(fragment_particles) == 0:
827 print(
"add_component_density: no particle was selected")
830 density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
839 multiply_by_total_mass,
845 s0.set_name(out_hier_name)
846 self.hier_representation[name][
"Densities"].add_child(s0)
848 for nps, p
in enumerate(density_particles):
850 p.set_name(s0.get_name() +
'_gaussian_%i' % nps)
853 def get_component_density(self, name):
854 return self.hier_representation[name][
"Densities"]
857 selection_tuples=
None,
862 '''Decorates all specified particles as Gaussians directly.
863 @param name component name
864 @param hierarchies set up GMM for some hierarchies
865 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
866 @param particles set up GMM for particles directly
867 @param resolution usual PMI resolution for selecting particles from the hierarchies
873 from math
import sqrt
874 self.representation_is_modified =
True
876 if particles
is None:
877 fragment_particles = []
879 fragment_particles = particles
881 if not hierarchies
is None:
883 self, resolution=resolution,
884 hierarchies=hierarchies)
886 if not selection_tuples
is None:
887 for st
in selection_tuples:
888 fragment_particles += IMP.pmi1.tools.select_by_tuple(
889 self, tupleselection=st,
890 resolution=resolution,
891 name_is_ambiguous=
False)
893 if len(fragment_particles) == 0:
894 print(
"add all atom densities: no particle was selected")
898 print(
'setting up all atom gaussians num_particles',len(fragment_particles))
899 for n,p
in enumerate(fragment_particles):
907 print(
'setting up particle',p.get_name(),
" as individual gaussian particle")
909 if not output_map
is None:
910 print(
'writing map to', output_map)
918 Make a copy of a hierarchy and append it to a component.
921 self.representation_is_modified =
True
922 protein_h = self.hier_dict[name]
923 hierclone = IMP.atom.create_clone(hierarchy)
924 hierclone.set_name(hierclone.get_name() +
"_clone")
925 protein_h.add_child(hierclone)
926 outhier.append(hierclone)
932 for n, pmain
in enumerate(psmain):
938 self.hier_db.add_particles(
954 def dump_particle_descriptors(self):
960 particles_attributes={}
961 floppy_body_attributes={}
967 hparent=h.get_parent()
968 while hparent != hroot:
969 hparent=h.get_parent()
970 name+=
"|"+hparent.get_name()
972 particles_attributes[name]={
"COORDINATES":numpy.array(
IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
978 rigid_body_attributes={}
979 for rb
in self.rigid_bodies:
981 rf=rb.get_reference_frame()
982 t=rf.get_transformation_to()
983 trans=t.get_translation()
985 rigid_body_attributes[name]={
"TRANSLATION":numpy.array(trans),
986 "ROTATION":numpy.array(rot.get_quaternion()),
987 "COORDINATES_NONRIGID_MEMBER":{},
988 "COORDINATES_RIGID_MEMBER":{}}
989 for mi
in rb.get_member_indexes():
990 rm=self.model.get_particle(mi)
992 name_part=rm.get_name()
994 rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
996 name_part=rm.get_name()
998 rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
1002 pickle.dump(particles_attributes,
1003 open(
"particles_attributes.pkl",
"wb"))
1004 pickle.dump(rigid_body_attributes,
1005 open(
"rigid_body_attributes.pkl",
"wb"))
1009 def load_particle_descriptors(self):
1015 particles_attributes = pickle.load(open(
"particles_attributes.pkl",
1017 rigid_body_attributes = pickle.load(open(
"rigid_body_attributes.pkl",
1027 hparent=h.get_parent()
1028 while hparent != hroot:
1029 hparent=h.get_parent()
1030 name+=
"|"+hparent.get_name()
1034 xyzr.set_coordinates(particles_attributes[name][
"COORDINATES"])
1040 for rb
in self.rigid_bodies:
1042 trans=rigid_body_attributes[name][
"TRANSLATION"]
1043 rot=rigid_body_attributes[name][
"ROTATION"]
1046 rb.set_reference_frame(rf)
1047 coor_nrm_ref=rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"]
1048 coor_rm_ref_dict=rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"]
1051 for mi
in rb.get_member_indexes():
1052 rm=self.model.get_particle(mi)
1054 name_part=rm.get_name()
1055 xyz=coor_nrm_ref[name_part]
1057 self.model.set_attribute(fk, rm,xyz[n])
1059 name_part=rm.get_name()
1061 coor_rm_model.append(
IMP.core.XYZ(rm).get_coordinates())
1062 if len(coor_rm_model)==0:
continue
1068 def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1069 """Print a warning if particle names in RMF and model don't match"""
1070 def match_any_suffix():
1072 suffixes = [
"pdb",
"bead_floppy_body_rigid_body_member_floppy_body",
1073 "bead_floppy_body_rigid_body_member",
1076 if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1078 if rmfname != reprname
and not match_any_suffix():
1079 print(
"set_coordinates_from_rmf: WARNING rmf particle and "
1080 "representation particle names don't match %s %s"
1081 % (rmfname, reprname))
1085 rmf_component_name=
None,
1086 check_number_particles=
True,
1087 representation_name_to_rmf_name_map=
None,
1089 skip_gaussian_in_rmf=
False,
1090 skip_gaussian_in_representation=
False,
1092 force_rigid_update=
False):
1093 '''Read and replace coordinates from an RMF file.
1094 Replace the coordinates of particles with the same name.
1095 It assumes that the RMF and the representation have the particles
1097 @param component_name Component name
1098 @param rmf_component_name Name of the component in the RMF file
1099 (if not specified, use `component_name`)
1100 @param representation_name_to_rmf_name_map a dictionary that map
1101 the original rmf particle name to the recipient particle component name
1102 @param save_file: save a file with the names of particles of the component
1103 @param force_rigid_update: update the coordinates of rigid bodies
1104 (normally this should be called before rigid bodies are set up)
1109 prots = IMP.pmi1.analysis.get_hiers_from_rmf(
1115 raise ValueError(
"cannot read hierarchy from rmf")
1119 if force_rigid_update:
1130 p, self.hier_dict.keys())
1131 if (protname
is None)
and (rmf_component_name
is not None):
1133 p, rmf_component_name)
1134 if (skip_gaussian_in_rmf):
1137 if (rmf_component_name
is not None)
and (protname == rmf_component_name):
1139 elif (rmf_component_name
is None)
and (protname == component_name):
1143 if (skip_gaussian_in_representation):
1154 reprnames=[p.get_name()
for p
in psrepr]
1155 rmfnames=[p.get_name()
for p
in psrmf]
1158 fl=open(component_name+
".txt",
"w")
1159 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1162 if check_number_particles
and not representation_name_to_rmf_name_map:
1163 if len(psrmf) != len(psrepr):
1164 fl=open(component_name+
".txt",
"w")
1165 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1166 raise ValueError(
"%s cannot proceed the rmf and the representation don't have the same number of particles; "
1167 "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1170 if not representation_name_to_rmf_name_map:
1171 for n, prmf
in enumerate(psrmf):
1173 prmfname = prmf.get_name()
1174 preprname = psrepr[n].get_name()
1175 if force_rigid_update:
1181 raise ValueError(
"component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1183 self._compare_rmf_repr_names(prmfname, preprname,
1192 rbm.set_internal_coordinates(xyz)
1194 rb = rbm.get_rigid_body()
1196 raise ValueError(
"Cannot handle nested "
1198 rb.set_reference_frame_lazy(tr)
1205 g=gprmf.get_gaussian()
1206 grepr.set_gaussian(g)
1209 repr_name_particle_map={}
1210 rmf_name_particle_map={}
1212 rmf_name_particle_map[p.get_name()]=p
1216 for prepr
in psrepr:
1218 prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1220 print(
"set_coordinates_from_rmf: WARNING missing particle name in representation_name_to_rmf_name_map, skipping...")
1229 If the root hierarchy does not exist, construct it.
1231 if "Res:" + str(int(resolution))
not in self.hier_representation[name]:
1233 root.set_name(name +
"_Res:" + str(int(resolution)))
1234 self.hier_representation[name][
1235 "Res:" + str(int(resolution))] = root
1236 protein_h.add_child(root)
1239 input_hierarchy, protein_h, type, color):
1241 Generate all needed coarse grained layers.
1243 @param name name of the protein
1244 @param resolutions list of resolutions
1245 @param protein_h root hierarchy
1246 @param input_hierarchy hierarchy to coarse grain
1247 @param type a string, typically "pdb" or "helix"
1251 if (1
in resolutions)
or (0
in resolutions):
1254 if 1
in resolutions:
1257 s1.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1259 self.hier_representation[name][
"Res:1"].add_child(s1)
1261 if 0
in resolutions:
1264 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1266 self.hier_representation[name][
"Res:0"].add_child(s0)
1269 if not isnucleicacid:
1272 atom_type=IMP.atom.AT_CA)
1276 atom_type=IMP.atom.AT_P)
1278 for p
in sel.get_selected_particles():
1280 if 0
in resolutions:
1282 resclone0 = IMP.atom.create_clone(resobject)
1284 s0.add_child(resclone0)
1285 self.hier_db.add_particles(
1289 resclone0.get_children())
1291 chil = resclone0.get_children()
1300 if 1
in resolutions:
1302 resclone1 = IMP.atom.create_clone_one(resobject)
1304 s1.add_child(resclone1)
1305 self.hier_db.add_particles(name, resindex, 1, [resclone1])
1309 prt = resclone1.get_particle()
1310 prt.set_name(
'%s_%i_%s' % (name, resindex, type))
1332 for r
in resolutions:
1333 if r != 0
and r != 1:
1339 chil = s.get_children()
1341 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1346 self.hier_representation[name][
1347 "Res:" + str(int(r))].add_child(s0)
1356 prt.set_name(
'%s_%i_%s' % (name, first, type))
1358 prt.set_name(
'%s_%i_%i_%s' % (name, first, last, type))
1360 self.hier_db.add_particles(name, kk, r, [prt])
1381 Get the hierarchies at the given resolution.
1383 The map between resolution and hierarchies is cached to accelerate
1384 the selection algorithm. The cache is invalidated when the
1385 representation was changed by any add_component_xxx.
1388 if self.representation_is_modified:
1389 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1391 self.hier_resolution[resolution] = rhbr
1392 self.representation_is_modified =
False
1395 if resolution
in self.hier_resolution:
1396 return self.hier_resolution[resolution]
1398 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1400 self.hier_resolution[resolution] = rhbr
1404 self, max_translation=300., max_rotation=2.0 * pi,
1405 avoidcollision=
True, cutoff=10.0, niterations=100,
1407 excluded_rigid_bodies=
None,
1408 ignore_initial_coordinates=
False,
1409 hierarchies_excluded_from_collision=
None):
1411 Shuffle configuration; used to restart the optimization.
1413 The configuration of the system is initialized by placing each
1414 rigid body and each bead randomly in a box with a side of
1415 `max_translation` angstroms, and far enough from each other to
1416 prevent any steric clashes. The rigid bodies are also randomly rotated.
1418 @param avoidcollision check if the particle/rigid body was
1419 placed close to another particle; uses the optional
1420 arguments cutoff and niterations
1421 @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1424 if excluded_rigid_bodies
is None:
1425 excluded_rigid_bodies = []
1426 if hierarchies_excluded_from_collision
is None:
1427 hierarchies_excluded_from_collision = []
1429 if len(self.rigid_bodies) == 0:
1430 print(
"shuffle_configuration: rigid bodies were not initialized")
1433 gcpf.set_distance(cutoff)
1435 hierarchies_excluded_from_collision_indexes = []
1441 hierarchies_excluded_from_collision_indexes.extend(
1445 if bounding_box
is not None:
1446 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1451 for h
in hierarchies_excluded_from_collision:
1452 hierarchies_excluded_from_collision_indexes.extend(
1456 allparticleindexes = list(
1457 set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1459 print(hierarchies_excluded_from_collision)
1460 print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1462 print(
'shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies),
'rigid bodies')
1463 for rb
in self.rigid_bodies:
1464 if rb
not in excluded_rigid_bodies:
1466 rbindexes = rb.get_member_particle_indexes()
1468 set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1469 otherparticleindexes = list(
1470 set(allparticleindexes) - set(rbindexes))
1472 if len(otherparticleindexes)
is None:
1476 while niter < niterations:
1477 if (ignore_initial_coordinates):
1483 if bounding_box
is not None:
1499 gcpf.get_close_pairs(
1501 otherparticleindexes,
1505 if (ignore_initial_coordinates):
1506 print (rb.get_name(),
IMP.core.XYZ(rb).get_coordinates())
1509 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1510 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
1511 if niter == niterations:
1512 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1516 print(
'shuffling', len(self.floppy_bodies),
'floppy bodies')
1517 for fb
in self.floppy_bodies:
1518 if (avoidcollision):
1523 otherparticleindexes = list(
1524 set(allparticleindexes) - set(fbindexes))
1525 if len(otherparticleindexes)
is None:
1543 while niter < niterations:
1544 if (ignore_initial_coordinates):
1550 if bounding_box
is not None:
1562 if (avoidcollision):
1565 gcpf.get_close_pairs(
1567 otherparticleindexes,
1571 if (ignore_initial_coordinates):
1572 print (fb.get_name(),
IMP.core.XYZ(fb).get_coordinates())
1575 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1576 print(
"shuffle_configuration: floppy body name: " + fb.get_name())
1577 if niter == niterations:
1578 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1582 def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1586 self.reference_structures[label] = (
1590 def get_all_rmsds(self):
1593 for label
in self.reference_structures:
1595 current_coordinates = [
IMP.core.XYZ(p).get_coordinates()
1596 for p
in self.reference_structures[label][1]]
1597 reference_coordinates = self.reference_structures[label][0]
1598 if len(reference_coordinates) != len(current_coordinates):
1599 print(
"calculate_all_rmsds: reference and actual coordinates are not the same")
1602 current_coordinates,
1603 reference_coordinates)
1605 reference_coordinates,
1606 current_coordinates)
1610 reference_coordinates,
1611 current_coordinates)
1612 rmsds[label +
"_GlobalRMSD"] = rmsd_global
1613 rmsds[label +
"_RelativeDRMS"] = rmsd_relative
1616 def setup_component_geometry(self, name, color=None, resolution=1.0):
1618 color = self.color_dict[name]
1622 self.hier_geometry_pairs[name] = []
1623 protein_h = self.hier_dict[name]
1625 pbr = IMP.pmi1.tools.sort_by_residues(pbr)
1627 for n
in range(len(pbr) - 1):
1628 self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1631 self, name, resolution=10, scale=1.0):
1633 Generate restraints between contiguous fragments in the hierarchy.
1634 The linkers are generated at resolution 10 by default.
1641 protein_h = self.hier_dict[name]
1643 frs = self.hier_db.get_preroot_fragments_by_resolution(
1649 start = fr.get_children()[0]
1654 end = fr.get_children()[-1]
1660 SortedSegments.append((start, end, startres))
1661 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1664 for x
in range(len(SortedSegments) - 1):
1665 last = SortedSegments[x][1]
1666 first = SortedSegments[x + 1][0]
1673 residuegap = firstresn - lastresn - 1
1674 if self.disorderedlength
and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1677 optdist = sqrt(5 / 3) * 1.93 * \
1678 (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1680 if self.upperharmonic:
1686 optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1687 if self.upperharmonic:
1693 pt0 = last.get_particle()
1694 pt1 = first.get_particle()
1697 print(
"Adding sequence connectivity restraint between", pt0.get_name(),
" and ", pt1.get_name(),
'of distance', optdist)
1698 sortedsegments_cr.add_restraint(r)
1699 self.linker_restraints_dict[
1700 "LinkerRestraint-" + pt0.get_name() +
"-" + pt1.get_name()] = r
1701 self.connected_intra_pairs.append((pt0, pt1))
1702 self.connected_intra_pairs.append((pt1, pt0))
1704 self.linker_restraints.add_restraint(sortedsegments_cr)
1705 self.linker_restraints.add_restraint(unmodeledregions_cr)
1707 self.sortedsegments_cr_dict[name] = sortedsegments_cr
1708 self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1710 def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1712 pts = IMP.pmi1.tools.ParticleToSampleList()
1713 for n, fb
in enumerate(self.floppy_bodies):
1714 pts.add_particle(fb,
"Floppy_Bodies", 1.0,
"Floppy_Body_" + str(n))
1715 if len(pts.get_particles_to_sample()) > 0:
1717 print(
"optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1720 print(
"optimize_floppy_bodies: no particle to optimize")
1723 nSymmetry=
None, skip_gaussian_in_clones=
False):
1725 The copies must not contain rigid bodies.
1726 The symmetry restraints are applied at each leaf.
1730 self.representation_is_modified =
True
1731 ncopies = len(copies) + 1
1734 for k
in range(len(copies)):
1735 if (nSymmetry
is None):
1736 rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1739 rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1741 rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1749 for n, p
in enumerate(main_hiers):
1750 if (skip_gaussian_in_clones):
1759 self.model.add_score_state(c)
1760 print(
"Completed setting " + str(maincopy) +
" as a reference for "
1761 + str(copies[k]) +
" by rotating it by "
1762 + str(rotation_angle / 2.0 / pi * 360)
1763 +
" degrees around the " + str(rotational_axis) +
" axis.")
1766 def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1769 self.representation_is_modified =
True
1771 mainparticles = particles_reference
1773 t=initial_transformation
1775 p.set_name(
"RigidBody_Symmetry")
1781 copyparticles = particles_copy
1785 for n, p
in enumerate(mainparticles):
1787 pc = copyparticles[n]
1789 mainpurged.append(p)
1795 copypurged.append(pc)
1802 for n, p
in enumerate(mainpurged):
1805 print(
"setting " + p.get_name() +
" as reference for " + pc.get_name())
1808 lc.add(pc.get_index())
1811 self.model.add_score_state(c)
1814 self.rigid_bodies.append(rb)
1815 self.rigid_body_symmetries.append(rb)
1816 rb.set_name(label+
".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1819 def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1820 longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1823 self.representation_is_modified =
True
1826 protein_h = self.hier_dict[maincopy]
1829 for ilc
in range(-longitudinal_copies, longitudinal_copies + 1):
1830 for iac
in range(axial_copies):
1831 copyname = maincopy +
"_a" + str(ilc) +
"_l" + str(iac)
1832 self.create_component(copyname, 0.0)
1833 for hier
in protein_h.get_children():
1835 copyhier = self.hier_dict[copyname]
1836 outhiers.append(copyhier)
1840 2 * pi / axial_copies * (float(iac)))
1841 translation_vector = tuple(
1842 [translation_value * float(ilc) * x
for x
in axis])
1843 print(translation_vector)
1848 for n, p
in enumerate(mainparts):
1855 lc.add(pc.get_index())
1857 self.model.add_score_state(c)
1863 Load coordinates in the current representation.
1864 This should be done only if the hierarchy self.prot is identical
1865 to the one as stored in the rmf i.e. all components were added.
1870 rh = RMF.open_rmf_file_read_only(rmfname)
1878 create the representation (i.e. hierarchies) from the rmf file.
1879 it will be stored in self.prot, which will be overwritten.
1880 load the coordinates from the rmf file at frameindex.
1882 rh = RMF.open_rmf_file_read_only(rmfname)
1888 for p
in self.prot.get_children():
1889 self.create_component(p.get_name())
1890 self.hier_dict[p.get_name()] = p
1892 still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1893 save floppy bodies in self.floppy_bodies
1894 get the connectivity restraints
1899 Construct a rigid body from hierarchies (and optionally particles).
1901 @param hiers list of hierarchies to make rigid
1902 @param particles list of particles to add to the rigid body
1905 if particles
is None:
1908 rigid_parts = set(particles)
1911 print(
"set_rigid_body_from_hierarchies> setting up a new rigid body")
1918 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1921 name += hier.get_name() +
"-"
1922 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1924 if len(list(rigid_parts)) != 0:
1926 rb.set_coordinates_are_optimized(
True)
1927 rb.set_name(name +
"rigid_body")
1928 self.rigid_bodies.append(rb)
1932 print(
"set_rigid_body_from_hierarchies> rigid body could not be setup")
1936 Construct a rigid body from a list of subunits.
1938 Each subunit is a tuple that identifies the residue ranges and the
1939 component name (as used in create_component()).
1941 subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1943 [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1945 example: ["prot1","prot2",("prot3",(1,10))]
1947 sometimes, we know about structure of an interaction
1948 and here we make such PPIs rigid
1953 if type(s) == type(tuple())
and len(s) == 2:
1957 residue_indexes=list(range(s[1][0],
1959 if len(sel.get_selected_particles()) == 0:
1960 print(
"set_rigid_bodies: selected particle does not exist")
1961 for p
in sel.get_selected_particles():
1965 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1969 elif type(s) == type(str()):
1971 if len(sel.get_selected_particles()) == 0:
1972 print(
"set_rigid_bodies: selected particle does not exist")
1973 for p
in sel.get_selected_particles():
1977 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1982 rb.set_coordinates_are_optimized(
True)
1983 rb.set_name(
''.join(str(subunits)) +
"_rigid_body")
1984 self.rigid_bodies.append(rb)
1987 def set_super_rigid_body_from_hierarchies(
1993 super_rigid_xyzs = set()
1994 super_rigid_rbs = set()
1996 print(
"set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
2003 super_rigid_rbs.add(rb)
2005 super_rigid_xyzs.add(p)
2006 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
2007 if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
2010 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2013 self.super_rigid_bodies.append(
2014 (super_rigid_xyzs, super_rigid_rbs, axis))
2016 def fix_rigid_bodies(self, rigid_bodies):
2017 self.fixed_rigid_bodies += rigid_bodies
2021 self, hiers, min_length=
None, max_length=
None, axis=
None):
2023 Make a chain of super rigid bodies from a list of hierarchies.
2025 Takes a linear list of hierarchies (they are supposed to be
2026 sequence-contiguous) and produces a chain of super rigid bodies
2027 with given length range, specified by min_length and max_length.
2030 hiers = IMP.pmi1.tools.flatten_list(hiers)
2034 self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2036 def set_super_rigid_bodies(self, subunits, coords=None):
2037 super_rigid_xyzs = set()
2038 super_rigid_rbs = set()
2041 if type(s) == type(tuple())
and len(s) == 3:
2045 residue_indexes=list(range(s[0],
2047 if len(sel.get_selected_particles()) == 0:
2048 print(
"set_rigid_bodies: selected particle does not exist")
2049 for p
in sel.get_selected_particles():
2052 super_rigid_rbs.add(rb)
2054 super_rigid_xyzs.add(p)
2055 elif type(s) == type(str()):
2057 if len(sel.get_selected_particles()) == 0:
2058 print(
"set_rigid_bodies: selected particle does not exist")
2059 for p
in sel.get_selected_particles():
2063 super_rigid_rbs.add(rb)
2065 super_rigid_xyzs.add(p)
2066 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2070 Remove leaves of hierarchies from the floppy bodies list based
2071 on the component name
2074 ps=[h.get_particle()
for h
in hs]
2075 for p
in self.floppy_bodies:
2077 if p
in ps: self.floppy_bodies.remove(p)
2078 if p
in hs: self.floppy_bodies.remove(p)
2085 Remove leaves of hierarchies from the floppy bodies list.
2087 Given a list of hierarchies, get the leaves and remove the
2088 corresponding particles from the floppy bodies list. We need this
2089 function because sometimes
2090 we want to constrain the floppy bodies in a rigid body - for instance
2091 when you want to associate a bead with a density particle.
2093 for h
in hierarchies:
2096 if p
in self.floppy_bodies:
2097 print(
"remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2098 self.floppy_bodies.remove(p)
2101 def set_floppy_bodies(self):
2102 for p
in self.floppy_bodies:
2104 p.set_name(name +
"_floppy_body")
2106 print(
"I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2109 rb.set_is_rigid_member(p.get_index(),
False)
2112 rb.set_is_rigid_member(p.get_particle_index(),
False)
2113 p.set_name(p.get_name() +
"_rigid_body_member")
2115 def set_floppy_bodies_from_hierarchies(self, hiers):
2120 self.floppy_bodies.append(p)
2127 selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2128 @return the particles
2131 print(selection_tuples)
2132 for s
in selection_tuples:
2133 ps = IMP.pmi1.tools.select_by_tuple(
2134 representation=self, tupleselection=tuple(s),
2135 resolution=
None, name_is_ambiguous=
False)
2139 def get_connected_intra_pairs(self):
2140 return self.connected_intra_pairs
2142 def set_rigid_bodies_max_trans(self, maxtrans):
2143 self.maxtrans_rb = maxtrans
2145 def set_rigid_bodies_max_rot(self, maxrot):
2146 self.maxrot_rb = maxrot
2148 def set_super_rigid_bodies_max_trans(self, maxtrans):
2149 self.maxtrans_srb = maxtrans
2151 def set_super_rigid_bodies_max_rot(self, maxrot):
2152 self.maxrot_srb = maxrot
2154 def set_floppy_bodies_max_trans(self, maxtrans):
2155 self.maxtrans_fb = maxtrans
2159 Fix rigid bodies in their actual position.
2160 The get_particles_to_sample() function will return
2161 just the floppy bodies (if they are not fixed).
2163 self.rigidbodiesarefixed = rigidbodiesarefixed
2167 Fix floppy bodies in their actual position.
2168 The get_particles_to_sample() function will return
2169 just the rigid bodies (if they are not fixed).
2171 self.floppybodiesarefixed=floppybodiesarefixed
2173 def draw_hierarchy_graph(self):
2175 print(
"Drawing hierarchy graph for " + c.get_name())
2176 IMP.pmi1.output.get_graph_from_hierarchy(c)
2178 def get_geometries(self):
2181 for name
in self.hier_geometry_pairs:
2182 for pt
in self.hier_geometry_pairs[name]:
2196 def setup_bonds(self):
2199 for name
in self.hier_geometry_pairs:
2200 for pt
in self.hier_geometry_pairs[name]:
2215 def show_component_table(self, name):
2216 if name
in self.sequence_dict:
2217 lastresn = len(self.sequence_dict[name])
2220 residues = self.hier_db.get_residue_numbers(name)
2221 firstresn = min(residues)
2222 lastresn = max(residues)
2224 for nres
in range(firstresn, lastresn + 1):
2226 resolutions = self.hier_db.get_residue_resolutions(name, nres)
2228 for r
in resolutions:
2229 ps += self.hier_db.get_particles(name, nres, r)
2230 print(
"%20s %7s" % (name, nres),
" ".join([
"%20s %7s" % (str(p.get_name()),
2233 print(
"%20s %20s" % (name, nres),
"**** not represented ****")
2235 def draw_hierarchy_composition(self):
2237 ks = list(self.elements.keys())
2241 for k
in self.elements:
2242 for l
in self.elements[k]:
2247 self.draw_component_composition(k, max)
2249 def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2250 from matplotlib
import pyplot
2251 import matplotlib
as mpl
2253 tmplist = sorted(self.elements[k], key=itemgetter(0))
2255 endres = tmplist[-1][1]
2257 print(
"draw_component_composition: missing information for component %s" % name)
2259 fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2260 ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2265 norm = mpl.colors.Normalize(vmin=5, vmax=10)
2269 for n, l
in enumerate(tmplist):
2273 if bounds[-1] != l[0]:
2274 colors.append(
"white")
2277 colors.append(
"#99CCFF")
2279 colors.append(
"#FFFF99")
2281 colors.append(
"#33CCCC")
2283 bounds.append(l[1] + 1)
2286 colors.append(
"#99CCFF")
2288 colors.append(
"#FFFF99")
2290 colors.append(
"#33CCCC")
2292 bounds.append(l[1] + 1)
2294 if bounds[-1] - 1 == l[0]:
2298 colors.append(
"white")
2301 bounds.append(bounds[-1])
2302 colors.append(
"white")
2303 cmap = mpl.colors.ListedColormap(colors)
2304 cmap.set_over(
'0.25')
2305 cmap.set_under(
'0.75')
2307 norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2308 cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2314 spacing=
'proportional',
2315 orientation=
'horizontal')
2324 mid = 1.0 / endres * float(l[0])
2329 l[2], xy=(mid, 1), xycoords=
'axes fraction',
2330 xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords=
'axes fraction',
2331 arrowprops=dict(arrowstyle=
"->",
2332 connectionstyle=
"angle,angleA=0,angleB=90,rad=10"),
2334 extra_artists.append(t)
2337 title = ax.text(-0.005, 0.5, k, ha=
"right", va=
"center", rotation=90,
2340 extra_artists.append(title)
2342 labels = len(bounds) * [
" "]
2343 ax.set_xticklabels(labels)
2344 mid = 1.0 / endres * float(bounds[0])
2345 t = ax.annotate(bounds[0], xy=(mid, 0), xycoords=
'axes fraction',
2346 xytext=(mid - 0.01, -0.5), textcoords=
'axes fraction',)
2347 extra_artists.append(t)
2348 offsets = [0, -0.5, -1.0]
2350 for n
in range(1, len(bounds)):
2351 if bounds[n] == bounds[n - 1]:
2353 mid = 1.0 / endres * float(bounds[n])
2354 if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2356 offset = offsets[nclashes % 3]
2362 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2363 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction')
2366 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2367 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction', arrowprops=dict(arrowstyle=
"-"))
2368 extra_artists.append(t)
2370 cb2.add_lines(bounds, [
"black"] * len(bounds), [1] * len(bounds))
2374 k +
"structure.pdf",
2377 bbox_extra_artists=(extra_artists),
2378 bbox_inches=
'tight')
2381 def draw_coordinates_projection(self):
2382 import matplotlib.pyplot
as pp
2385 for name
in self.hier_geometry_pairs:
2386 for pt
in self.hier_geometry_pairs[name]:
2396 xpairs.append([x1, x2])
2397 ypairs.append([y1, y2])
2400 for xends, yends
in zip(xpairs, ypairs):
2405 pp.plot(xlist, ylist,
'b-', alpha=0.1)
2408 def get_prot_name_from_particle(self, particle):
2409 names = self.get_component_names()
2410 particle0 = particle
2412 while not name
in names:
2415 particle0 = h.get_particle()
2419 def get_particles_to_sample(self):
2429 if not self.rigidbodiesarefixed:
2430 for rb
in self.rigid_bodies:
2435 if rb
not in self.fixed_rigid_bodies:
2438 if not self.floppybodiesarefixed:
2439 for fb
in self.floppy_bodies:
2446 for srb
in self.super_rigid_bodies:
2449 rigid_bodies = list(srb[1])
2450 filtered_rigid_bodies = []
2451 for rb
in rigid_bodies:
2452 if rb
not in self.fixed_rigid_bodies:
2453 filtered_rigid_bodies.append(rb)
2454 srbtmp.append((srb[0], filtered_rigid_bodies))
2456 self.rigid_bodies = rbtmp
2457 self.floppy_bodies = fbtmp
2458 self.super_rigid_bodies = srbtmp
2460 ps[
"Rigid_Bodies_SimplifiedModel"] = (
2464 ps[
"Floppy_Bodies_SimplifiedModel"] = (
2467 ps[
"SR_Bodies_SimplifiedModel"] = (
2468 self.super_rigid_bodies,
2473 def set_output_level(self, level):
2474 self.output_level = level
2476 def _evaluate(self, deriv):
2477 """Evaluate the total score of all added restraints"""
2479 return r.evaluate(deriv)
2481 def get_output(self):
2485 output[
"SimplifiedModel_Total_Score_" +
2486 self.label] = str(self._evaluate(
False))
2487 output[
"SimplifiedModel_Linker_Score_" +
2488 self.label] = str(self.linker_restraints.unprotected_evaluate(
None))
2489 for name
in self.sortedsegments_cr_dict:
2490 partialscore = self.sortedsegments_cr_dict[name].evaluate(
False)
2491 score += partialscore
2493 "SimplifiedModel_Link_SortedSegments_" +
2498 partialscore = self.unmodeledregions_cr_dict[name].evaluate(
False)
2499 score += partialscore
2501 "SimplifiedModel_Link_UnmodeledRegions_" +
2506 for rb
in self.rigid_body_symmetries:
2508 output[name +
"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2509 for name
in self.linker_restraints_dict:
2514 self.linker_restraints_dict[
2515 name].unprotected_evaluate(
2518 if len(self.reference_structures.keys()) != 0:
2519 rmsds = self.get_all_rmsds()
2522 "SimplifiedModel_" +
2525 self.label] = rmsds[
2528 if self.output_level ==
"high":
2531 output[
"Coordinates_" +
2532 p.get_name() +
"_" + self.label] = str(d)
2534 output[
"_TotalScore"] = str(score)
2537 def get_test_output(self):
2539 output = self.get_output()
2540 for n, p
in enumerate(self.get_particles_to_sample()):
2541 output[
"Particle_to_sample_" + str(n)] = str(p)
2543 output[
"Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2544 output[
"Number_of_floppy_bodies"] = len(self.floppy_bodies)
2545 output[
"Number_of_rigid_bodies"] = len(self.rigid_bodies)
2546 output[
"Number_of_super_bodies"] = len(self.super_rigid_bodies)
2547 output[
"Selection_resolution_1"] = len(
2549 output[
"Selection_resolution_5"] = len(
2551 output[
"Selection_resolution_7"] = len(
2553 output[
"Selection_resolution_10"] = len(
2555 output[
"Selection_resolution_100"] = len(
2558 output[
"Selection_resolution=1"] = len(
2560 output[
"Selection_resolution=1,resid=10"] = len(
2562 for resolution
in self.hier_resolution:
2563 output[
"Hier_resolution_dictionary" +
2564 str(resolution)] = len(self.hier_resolution[resolution])
2565 for name
in self.hier_dict:
2567 "Selection_resolution=1,resid=10,name=" +
2575 "Selection_resolution=1,resid=10,name=" +
2577 ",ambiguous"] = len(
2582 name_is_ambiguous=
True,
2585 "Selection_resolution=1,resid=10,name=" +
2587 ",ambiguous"] = len(
2592 name_is_ambiguous=
True,
2595 "Selection_resolution=1,resrange=(10,20),name=" +
2604 "Selection_resolution=1,resrange=(10,20),name=" +
2606 ",ambiguous"] = len(
2611 name_is_ambiguous=
True,
2615 "Selection_resolution=10,resrange=(10,20),name=" +
2624 "Selection_resolution=10,resrange=(10,20),name=" +
2626 ",ambiguous"] = len(
2631 name_is_ambiguous=
True,
2635 "Selection_resolution=100,resrange=(10,20),name=" +
2644 "Selection_resolution=100,resrange=(10,20),name=" +
2646 ",ambiguous"] = len(
2651 name_is_ambiguous=
True,
static Symmetric setup_particle(Model *m, ParticleIndex pi, Float symmetric)
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.
static Resolution setup_particle(Model *m, ParticleIndex pi, Float resolution)
def add_protocol_output
Capture details of the modeling protocol.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A decorator to associate a particle with a part of a protein/DNA/RNA.
def add_component_pdb
Add a component that has an associated 3D structure in a PDB file.
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.
static bool get_is_setup(Model *m, ParticleIndex pi)
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
A member of a rigid body, it has internal (local) coordinates.
Add resolution to a particle.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
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)
def coarse_hierarchy
Generate all needed coarse grained layers.
static XYZR setup_particle(Model *m, ParticleIndex pi)
A decorator for a particle which has bonds.
def create_components_from_rmf
still not working.
def check_root
If the root hierarchy does not exist, construct it.
Select atoms which are selected by both selectors.
static bool get_is_setup(Model *m, ParticleIndex pi)
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 add_component_hierarchy_clone
Make a copy of a hierarchy and append it to a component.
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.
def link_components_to_rmf
Load coordinates in the current representation.
Legacy PMI1 module to represent, score, sample and analyze models.
Provenance create_clone(Provenance p)
Clone provenance (including previous provenance)
def create_transformed_component
Create a transformed view of a component.
double get_ball_radius_from_volume_3d(double volume)
Return the radius of a sphere with a given volume.
Set of python classes to create a multi-state, multi-resolution IMP hierarchy.
def get_particles_from_selection_tuples
selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....
def get_hierarchies_at_given_resolution
Get the hierarchies at the given resolution.
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)
def shuffle_configuration
Shuffle configuration; used to restart the optimization.
def setup_component_sequence_connectivity
Generate restraints between contiguous fragments in the hierarchy.
static XYZ setup_particle(Model *m, ParticleIndex pi)
void read_pdb(TextInput input, int model, Hierarchy h)
def add_component_beads
add beads to the representation
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
Sample using Monte Carlo.
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.
Tools for clustering and cluster analysis.
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 add_component_density
Sets up a Gaussian Mixture Model for this component.
Select all CB ATOM records.
Add uncertainty to a particle.
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.
static bool get_is_setup(Model *m, ParticleIndex pi)
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
A dictionary-like wrapper for reading and storing sequence data.
def add_metadata
Associate some metadata with this modeling.
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)
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def get_file_dataset
Get the dataset associated with a filename, or None.
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.
def set_coordinates_from_rmf
Read and replace coordinates from an RMF file.
A decorator for a particle with x,y,z coordinates.
def set_file_dataset
Associate a dataset with a filename.
def add_all_atom_densities
Decorates all specified particles as Gaussians directly.
def add_component_sequence
Add the primary sequence for a single component.
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
def set_rigid_body_from_hierarchies
Construct a rigid body from hierarchies (and optionally particles).
static Uncertainty setup_particle(Model *m, ParticleIndex pi, Float uncertainty)
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Find all nearby pairs by testing all pairs.
Score a pair of particles based on the distance between their centers.
Simple implementation of segments in 3D.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def remove_floppy_bodies
Remove leaves of hierarchies from the floppy bodies list.
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...
def set_chain_of_super_rigid_bodies
Make a chain of super rigid bodies from a list of hierarchies.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def set_rigid_bodies
Construct a rigid body from a list of subunits.
def add_component_necklace
Generates a string of beads with given length.
The general base class for IMP exceptions.
def create_rotational_symmetry
The copies must not contain rigid bodies.
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.
Add symmetric attribute to a particle.
Store info for a chain of a protein.
Applies a PairScore to a Pair.
Classes for writing output files and processing them.
Select all CA ATOM records.
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)
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)
An exception for an invalid value being passed to IMP.
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.
def create_non_modeled_component
Create a component that isn't used in the modeling.
Support for the RMF file format for storing hierarchical molecular data and markup.
static bool get_is_setup(Model *m, ParticleIndex pi)
def remove_floppy_bodies_from_component
Remove leaves of hierarchies from the floppy bodies list based on the component name.
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 set_rigid_bodies_as_fixed
Fix rigid bodies in their actual position.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
def set_floppy_bodies_as_fixed
Fix floppy bodies in their actual position.
Harmonic function (symmetric about the mean)
A decorator for a particle with x,y,z coordinates and a radius.