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)))
322 t.get_children()[0].get_children()[0]).
get_index()
324 t.get_children()[0].get_children()[-1]).
get_index()
330 if name
in self.sequence_dict:
331 resrange = (1, len(self.sequence_dict[name]))
333 resrange = (start + offset, end + offset)
335 if resrange[1]
in (-1,
'END'):
336 resrange = (resrange[0],end)
337 start = resrange[0] - offset
338 end = resrange[1] - offset
357 for n, g
in enumerate(gaps):
361 print(
"autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
363 chain, resolutions=resolutions,
364 color=color, cacenters=
True,
365 resrange=(first, last),
366 offset=offset, isnucleicacid=isnucleicacid)
367 elif g[2] ==
"gap" and n > 0:
368 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
369 parts = self.hier_db.get_particles_at_closest_resolution(name,
374 first+offset, last+offset, missingbeadsize, incoord=xyz)
376 elif g[2] ==
"gap" and n == 0:
378 print(
"autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
380 first+offset, last+offset, missingbeadsize, incoord=xyznter)
387 def autobuild_pdb_and_intervening_beads(self, *args, **kwargs):
388 r = self.autobuild_model(*args, **kwargs)
392 resrange=
None, offset=0, cacenters=
True, show=
False,
393 isnucleicacid=
False, readnonwateratoms=
False,
394 read_ca_cb_only=
False):
396 Add a component that has an associated 3D structure in a PDB file.
398 Reads the PDB, and constructs the fragments corresponding to contiguous
401 @return a list of hierarchies.
403 @param name (string) the name of the component
404 @param pdbname (string) the name of the PDB file
405 @param chain (string or integer) can be either a string (eg, "A")
406 or an integer (eg, 0 or 1) in case you want
407 to get the corresponding chain number in the PDB.
408 @param resolutions (integers) a list of integers that corresponds
409 to the resolutions that have to be generated
410 @param color (float from 0 to 1) the color applied to the
411 hierarchies generated
412 @param resrange (tuple of integers): the residue range to extract
413 from the PDB. It is a tuple (beg,end). If not specified,
414 all residues belonging to the specified chain are read.
415 @param offset (integer) specifies the residue index offset to be
416 applied when reading the PDB (the FASTA sequence is
417 assumed to start from residue 1, so use this parameter
418 if the PDB file does not start at residue 1)
419 @param cacenters (boolean) if True generates resolution=1 beads
420 centered on C-alpha atoms.
421 @param show (boolean) print out the molecular hierarchy at the end.
422 @param isnucleicacid (boolean) use True if you're reading a PDB
424 @param readnonwateratoms (boolean) if True fixes some pathological PDB.
425 @param read_ca_cb_only (boolean) if True, only reads CA/CB
428 self.representation_is_modified =
True
431 color = self.color_dict[name]
432 protein_h = self.hier_dict[name]
442 if type(chain) == str:
450 t.get_children()[0].get_children()[0]).
get_index()
452 t.get_children()[0].get_children()[-1]).
get_index()
453 c =
IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
457 IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
464 if not resrange
is None:
465 if resrange[0] > start:
467 if resrange[1] < end:
470 if not isnucleicacid:
474 residue_indexes=list(range(
477 atom_type=IMP.atom.AT_CA)
484 residue_indexes=list(range(
487 atom_type=IMP.atom.AT_P)
489 ps = sel.get_selected_particles()
491 raise ValueError(
"%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
492 % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
500 if par.get_parent() != c0:
501 par.get_parent().remove_child(par)
503 start = start + offset
506 self.elements[name].append(
507 (start, end, pdbname.split(
"/")[-1] +
":" + chain,
"pdb"))
510 resolutions, isnucleicacid, c0, protein_h,
"pdb", color)
511 self._copy_pdb_provenance(t, hiers[0], offset)
513 for p, state
in self._protocol_output:
514 p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
528 for r
in residues.keys():
530 self.model.remove_particle(c0)
536 def _copy_pdb_provenance(self, pdb, h, offset):
537 """Copy the provenance information from the PDB to our hierarchy"""
538 c =
IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
545 def add_component_ideal_helix(
553 self.representation_is_modified =
True
554 from math
import pi, cos, sin
556 protein_h = self.hier_dict[name]
559 color = self.color_dict[name]
563 self.elements[name].append((start, end,
" ",
"helix"))
565 for n, res
in enumerate(range(start, end + 1)):
566 if name
in self.sequence_dict:
568 rtstr = self.onetothree[
569 self.sequence_dict[name][res-1]]
589 x = 2.3 * cos(n * 2 * pi / 3.6)
590 y = 2.3 * sin(n * 2 * pi / 3.6)
591 z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
600 resolutions,
False, c0, protein_h,
"helix", color)
609 """ add beads to the representation
610 @param name the component name
611 @param ds a list of tuples corresponding to the residue ranges
613 @param colors a list of colors associated to the beads
614 @param incoord the coordinate tuple corresponding to the position
619 self.representation_is_modified =
True
621 protein_h = self.hier_dict[name]
624 colors = [self.color_dict[name]]
627 for n, dss
in enumerate(ds):
628 ds_frag = (dss[0], dss[1])
629 self.elements[name].append((dss[0], dss[1],
" ",
"bead"))
631 if ds_frag[0] == ds_frag[1]:
633 if name
in self.sequence_dict:
635 rtstr = self.onetothree[
636 self.sequence_dict[name][ds_frag[0]-1]]
643 h.set_name(name +
'_%i_bead' % (ds_frag[0]))
644 prt.set_name(name +
'_%i_bead' % (ds_frag[0]))
648 h.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
649 prt.set_name(name +
'_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
650 h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
651 resolution = len(h.get_residue_indexes())
652 if "Beads" not in self.hier_representation[name]:
654 root.set_name(
"Beads")
655 self.hier_representation[name][
"Beads"] = root
656 protein_h.add_child(root)
657 self.hier_representation[name][
"Beads"].add_child(h)
659 for kk
in range(ds_frag[0], ds_frag[1] + 1):
660 self.hier_db.add_particles(name, kk, resolution, [h])
683 ptem.set_radius(radius)
686 radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
688 ptem.set_radius(radius)
690 if not tuple(incoord)
is None:
691 ptem.set_coordinates(incoord)
696 self.floppy_bodies.append(prt)
700 for p, state
in self._protocol_output:
701 p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
709 Generates a string of beads with given length.
711 self.representation_is_modified =
True
719 [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
723 self, name, hierarchies=
None, selection_tuples=
None,
725 resolution=0.0, num_components=10,
726 inputfile=
None, outputfile=
None,
729 covariance_type=
'full', voxel_size=1.0,
731 sampled_points=1000000, num_iter=100,
733 multiply_by_total_mass=
True,
735 intermediate_map_fn=
None,
736 density_ps_to_copy=
None,
737 use_precomputed_gaussians=
False):
739 Sets up a Gaussian Mixture Model for this component.
740 Can specify input GMM file or it will be computed.
741 @param name component name
742 @param hierarchies set up GMM for some hierarchies
743 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
744 @param particles set up GMM for particles directly
745 @param resolution usual PMI resolution for selecting particles from the hierarchies
746 @param inputfile read the GMM from this file
747 @param outputfile fit and write the GMM to this file (text)
748 @param outputmap after fitting, create GMM density file (mrc)
749 @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
750 @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
751 @param voxel_size for creating the intermediate density map and output map.
752 lower number increases accuracy but also rasterizing time grows
753 @param out_hier_name name of the output density hierarchy
754 @param sampled_points number of points to sample. more will increase accuracy and fitting time
755 @param num_iter num GMM iterations. more will increase accuracy and fitting time
756 @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
757 @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
758 @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
759 @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
760 @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
768 self.representation_is_modified =
True
770 protein_h = self.hier_dict[name]
771 if "Densities" not in self.hier_representation[name]:
773 root.set_name(
"Densities")
774 self.hier_representation[name][
"Densities"] = root
775 protein_h.add_child(root)
778 fragment_particles = []
779 if not particles
is None:
780 fragment_particles += particles
781 if not hierarchies
is None:
783 self, resolution=resolution,
784 hierarchies=hierarchies)
785 if not selection_tuples
is None:
786 for st
in selection_tuples:
787 fragment_particles += IMP.pmi1.tools.select_by_tuple(
788 self, tupleselection=st,
789 resolution=resolution,
790 name_is_ambiguous=
False)
793 density_particles = []
796 inputfile, density_particles,
797 self.model, transform)
798 elif density_ps_to_copy:
799 for ip
in density_ps_to_copy:
805 density_particles.append(p)
806 elif use_precomputed_gaussians:
807 if len(fragment_particles) == 0:
808 print(
"add_component_density: no particle was selected")
810 for p
in fragment_particles:
813 raise Exception(
"The particles you selected must be Fragments and XYZs")
814 nres=len(
IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
815 pos=
IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
820 raise Exception(
"We haven't created a bead file for",nres,
"residues")
824 self.model, transform)
827 if len(fragment_particles) == 0:
828 print(
"add_component_density: no particle was selected")
831 density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
840 multiply_by_total_mass,
846 s0.set_name(out_hier_name)
847 self.hier_representation[name][
"Densities"].add_child(s0)
849 for nps, p
in enumerate(density_particles):
851 p.set_name(s0.get_name() +
'_gaussian_%i' % nps)
854 def get_component_density(self, name):
855 return self.hier_representation[name][
"Densities"]
858 selection_tuples=
None,
863 '''Decorates all specified particles as Gaussians directly.
864 @param name component name
865 @param hierarchies set up GMM for some hierarchies
866 @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
867 @param particles set up GMM for particles directly
868 @param resolution usual PMI resolution for selecting particles from the hierarchies
874 from math
import sqrt
875 self.representation_is_modified =
True
877 if particles
is None:
878 fragment_particles = []
880 fragment_particles = particles
882 if not hierarchies
is None:
884 self, resolution=resolution,
885 hierarchies=hierarchies)
887 if not selection_tuples
is None:
888 for st
in selection_tuples:
889 fragment_particles += IMP.pmi1.tools.select_by_tuple(
890 self, tupleselection=st,
891 resolution=resolution,
892 name_is_ambiguous=
False)
894 if len(fragment_particles) == 0:
895 print(
"add all atom densities: no particle was selected")
899 print(
'setting up all atom gaussians num_particles',len(fragment_particles))
900 for n,p
in enumerate(fragment_particles):
908 print(
'setting up particle',p.get_name(),
" as individual gaussian particle")
910 if not output_map
is None:
911 print(
'writing map to', output_map)
919 Make a copy of a hierarchy and append it to a component.
922 self.representation_is_modified =
True
923 protein_h = self.hier_dict[name]
924 hierclone = IMP.atom.create_clone(hierarchy)
925 hierclone.set_name(hierclone.get_name() +
"_clone")
926 protein_h.add_child(hierclone)
927 outhier.append(hierclone)
933 for n, pmain
in enumerate(psmain):
939 self.hier_db.add_particles(
955 def dump_particle_descriptors(self):
961 particles_attributes={}
962 floppy_body_attributes={}
968 hparent=h.get_parent()
969 while hparent != hroot:
970 hparent=h.get_parent()
971 name+=
"|"+hparent.get_name()
973 particles_attributes[name]={
"COORDINATES":numpy.array(
IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
979 rigid_body_attributes={}
980 for rb
in self.rigid_bodies:
982 rf=rb.get_reference_frame()
983 t=rf.get_transformation_to()
984 trans=t.get_translation()
986 rigid_body_attributes[name]={
"TRANSLATION":numpy.array(trans),
987 "ROTATION":numpy.array(rot.get_quaternion()),
988 "COORDINATES_NONRIGID_MEMBER":{},
989 "COORDINATES_RIGID_MEMBER":{}}
990 for mi
in rb.get_member_indexes():
991 rm=self.model.get_particle(mi)
993 name_part=rm.get_name()
995 rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
997 name_part=rm.get_name()
999 rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
1003 pickle.dump(particles_attributes,
1004 open(
"particles_attributes.pkl",
"wb"))
1005 pickle.dump(rigid_body_attributes,
1006 open(
"rigid_body_attributes.pkl",
"wb"))
1010 def load_particle_descriptors(self):
1016 particles_attributes = pickle.load(open(
"particles_attributes.pkl",
1018 rigid_body_attributes = pickle.load(open(
"rigid_body_attributes.pkl",
1028 hparent=h.get_parent()
1029 while hparent != hroot:
1030 hparent=h.get_parent()
1031 name+=
"|"+hparent.get_name()
1035 xyzr.set_coordinates(particles_attributes[name][
"COORDINATES"])
1041 for rb
in self.rigid_bodies:
1043 trans=rigid_body_attributes[name][
"TRANSLATION"]
1044 rot=rigid_body_attributes[name][
"ROTATION"]
1047 rb.set_reference_frame(rf)
1048 coor_nrm_ref=rigid_body_attributes[name][
"COORDINATES_NONRIGID_MEMBER"]
1049 coor_rm_ref_dict=rigid_body_attributes[name][
"COORDINATES_RIGID_MEMBER"]
1052 for mi
in rb.get_member_indexes():
1053 rm=self.model.get_particle(mi)
1055 name_part=rm.get_name()
1056 xyz=coor_nrm_ref[name_part]
1058 self.model.set_attribute(fk, rm,xyz[n])
1060 name_part=rm.get_name()
1062 coor_rm_model.append(
IMP.core.XYZ(rm).get_coordinates())
1063 if len(coor_rm_model)==0:
continue
1069 def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1070 """Print a warning if particle names in RMF and model don't match"""
1071 def match_any_suffix():
1073 suffixes = [
"pdb",
"bead_floppy_body_rigid_body_member_floppy_body",
1074 "bead_floppy_body_rigid_body_member",
1077 if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1079 if rmfname != reprname
and not match_any_suffix():
1080 print(
"set_coordinates_from_rmf: WARNING rmf particle and "
1081 "representation particle names don't match %s %s"
1082 % (rmfname, reprname))
1086 rmf_component_name=
None,
1087 check_number_particles=
True,
1088 representation_name_to_rmf_name_map=
None,
1090 skip_gaussian_in_rmf=
False,
1091 skip_gaussian_in_representation=
False,
1093 force_rigid_update=
False):
1094 '''Read and replace coordinates from an RMF file.
1095 Replace the coordinates of particles with the same name.
1096 It assumes that the RMF and the representation have the particles
1098 @param component_name Component name
1099 @param rmf_component_name Name of the component in the RMF file
1100 (if not specified, use `component_name`)
1101 @param representation_name_to_rmf_name_map a dictionary that map
1102 the original rmf particle name to the recipient particle component name
1103 @param save_file: save a file with the names of particles of the component
1104 @param force_rigid_update: update the coordinates of rigid bodies
1105 (normally this should be called before rigid bodies are set up)
1110 prots = IMP.pmi1.analysis.get_hiers_from_rmf(
1116 raise ValueError(
"cannot read hierarchy from rmf")
1120 if force_rigid_update:
1131 p, self.hier_dict.keys())
1132 if (protname
is None)
and (rmf_component_name
is not None):
1134 p, rmf_component_name)
1135 if (skip_gaussian_in_rmf):
1138 if (rmf_component_name
is not None)
and (protname == rmf_component_name):
1140 elif (rmf_component_name
is None)
and (protname == component_name):
1144 if (skip_gaussian_in_representation):
1155 reprnames=[p.get_name()
for p
in psrepr]
1156 rmfnames=[p.get_name()
for p
in psrmf]
1159 fl=open(component_name+
".txt",
"w")
1160 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1163 if check_number_particles
and not representation_name_to_rmf_name_map:
1164 if len(psrmf) != len(psrepr):
1165 fl=open(component_name+
".txt",
"w")
1166 for i
in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+
","+str(i[1])+
"\n")
1167 raise ValueError(
"%s cannot proceed the rmf and the representation don't have the same number of particles; "
1168 "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1171 if not representation_name_to_rmf_name_map:
1172 for n, prmf
in enumerate(psrmf):
1174 prmfname = prmf.get_name()
1175 preprname = psrepr[n].get_name()
1176 if force_rigid_update:
1182 raise ValueError(
"component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1184 self._compare_rmf_repr_names(prmfname, preprname,
1193 rbm.set_internal_coordinates(xyz)
1195 rb = rbm.get_rigid_body()
1197 raise ValueError(
"Cannot handle nested "
1199 rb.set_reference_frame_lazy(tr)
1206 g=gprmf.get_gaussian()
1207 grepr.set_gaussian(g)
1210 repr_name_particle_map={}
1211 rmf_name_particle_map={}
1213 rmf_name_particle_map[p.get_name()]=p
1217 for prepr
in psrepr:
1219 prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1221 print(
"set_coordinates_from_rmf: WARNING missing particle name in representation_name_to_rmf_name_map, skipping...")
1230 If the root hierarchy does not exist, construct it.
1232 if "Res:" + str(int(resolution))
not in self.hier_representation[name]:
1234 root.set_name(name +
"_Res:" + str(int(resolution)))
1235 self.hier_representation[name][
1236 "Res:" + str(int(resolution))] = root
1237 protein_h.add_child(root)
1240 input_hierarchy, protein_h, type, color):
1242 Generate all needed coarse grained layers.
1244 @param name name of the protein
1245 @param resolutions list of resolutions
1246 @param protein_h root hierarchy
1247 @param input_hierarchy hierarchy to coarse grain
1248 @param type a string, typically "pdb" or "helix"
1252 if (1
in resolutions)
or (0
in resolutions):
1255 if 1
in resolutions:
1258 s1.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1260 self.hier_representation[name][
"Res:1"].add_child(s1)
1262 if 0
in resolutions:
1265 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1267 self.hier_representation[name][
"Res:0"].add_child(s0)
1270 if not isnucleicacid:
1273 atom_type=IMP.atom.AT_CA)
1277 atom_type=IMP.atom.AT_P)
1279 for p
in sel.get_selected_particles():
1281 if 0
in resolutions:
1283 resclone0 = IMP.atom.create_clone(resobject)
1285 s0.add_child(resclone0)
1286 self.hier_db.add_particles(
1290 resclone0.get_children())
1292 chil = resclone0.get_children()
1301 if 1
in resolutions:
1303 resclone1 = IMP.atom.create_clone_one(resobject)
1305 s1.add_child(resclone1)
1306 self.hier_db.add_particles(name, resindex, 1, [resclone1])
1310 prt = resclone1.get_particle()
1311 prt.set_name(
'%s_%i_%s' % (name, resindex, type))
1333 for r
in resolutions:
1334 if r != 0
and r != 1:
1340 chil = s.get_children()
1342 s0.set_name(
'%s_%i-%i_%s' % (name, start, end, type))
1347 self.hier_representation[name][
1348 "Res:" + str(int(r))].add_child(s0)
1357 prt.set_name(
'%s_%i_%s' % (name, first, type))
1359 prt.set_name(
'%s_%i_%i_%s' % (name, first, last, type))
1361 self.hier_db.add_particles(name, kk, r, [prt])
1382 Get the hierarchies at the given resolution.
1384 The map between resolution and hierarchies is cached to accelerate
1385 the selection algorithm. The cache is invalidated when the
1386 representation was changed by any add_component_xxx.
1389 if self.representation_is_modified:
1390 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1392 self.hier_resolution[resolution] = rhbr
1393 self.representation_is_modified =
False
1396 if resolution
in self.hier_resolution:
1397 return self.hier_resolution[resolution]
1399 rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1401 self.hier_resolution[resolution] = rhbr
1405 self, max_translation=300., max_rotation=2.0 * pi,
1406 avoidcollision=
True, cutoff=10.0, niterations=100,
1408 excluded_rigid_bodies=
None,
1409 ignore_initial_coordinates=
False,
1410 hierarchies_excluded_from_collision=
None):
1412 Shuffle configuration; used to restart the optimization.
1414 The configuration of the system is initialized by placing each
1415 rigid body and each bead randomly in a box with a side of
1416 `max_translation` angstroms, and far enough from each other to
1417 prevent any steric clashes. The rigid bodies are also randomly rotated.
1419 @param avoidcollision check if the particle/rigid body was
1420 placed close to another particle; uses the optional
1421 arguments cutoff and niterations
1422 @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1425 if excluded_rigid_bodies
is None:
1426 excluded_rigid_bodies = []
1427 if hierarchies_excluded_from_collision
is None:
1428 hierarchies_excluded_from_collision = []
1430 if len(self.rigid_bodies) == 0:
1431 print(
"shuffle_configuration: rigid bodies were not intialized")
1434 gcpf.set_distance(cutoff)
1436 hierarchies_excluded_from_collision_indexes = []
1442 hierarchies_excluded_from_collision_indexes.extend(
1446 if bounding_box
is not None:
1447 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1452 for h
in hierarchies_excluded_from_collision:
1453 hierarchies_excluded_from_collision_indexes.extend(
1457 allparticleindexes = list(
1458 set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1460 print(hierarchies_excluded_from_collision)
1461 print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1463 print(
'shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies),
'rigid bodies')
1464 for rb
in self.rigid_bodies:
1465 if rb
not in excluded_rigid_bodies:
1467 rbindexes = rb.get_member_particle_indexes()
1469 set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1470 otherparticleindexes = list(
1471 set(allparticleindexes) - set(rbindexes))
1473 if len(otherparticleindexes)
is None:
1477 while niter < niterations:
1478 if (ignore_initial_coordinates):
1484 if bounding_box
is not None:
1500 gcpf.get_close_pairs(
1502 otherparticleindexes,
1506 if (ignore_initial_coordinates):
1507 print (rb.get_name(),
IMP.core.XYZ(rb).get_coordinates())
1510 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1511 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
1512 if niter == niterations:
1513 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1517 print(
'shuffling', len(self.floppy_bodies),
'floppy bodies')
1518 for fb
in self.floppy_bodies:
1519 if (avoidcollision):
1524 otherparticleindexes = list(
1525 set(allparticleindexes) - set(fbindexes))
1526 if len(otherparticleindexes)
is None:
1544 while niter < niterations:
1545 if (ignore_initial_coordinates):
1551 if bounding_box
is not None:
1563 if (avoidcollision):
1566 gcpf.get_close_pairs(
1568 otherparticleindexes,
1572 if (ignore_initial_coordinates):
1573 print (fb.get_name(),
IMP.core.XYZ(fb).get_coordinates())
1576 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1577 print(
"shuffle_configuration: floppy body name: " + fb.get_name())
1578 if niter == niterations:
1579 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1583 def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1587 self.reference_structures[label] = (
1591 def get_all_rmsds(self):
1594 for label
in self.reference_structures:
1596 current_coordinates = [
IMP.core.XYZ(p).get_coordinates()
1597 for p
in self.reference_structures[label][1]]
1598 reference_coordinates = self.reference_structures[label][0]
1599 if len(reference_coordinates) != len(current_coordinates):
1600 print(
"calculate_all_rmsds: reference and actual coordinates are not the same")
1603 current_coordinates,
1604 reference_coordinates)
1606 reference_coordinates,
1607 current_coordinates)
1611 reference_coordinates,
1612 current_coordinates)
1613 rmsds[label +
"_GlobalRMSD"] = rmsd_global
1614 rmsds[label +
"_RelativeDRMS"] = rmsd_relative
1617 def setup_component_geometry(self, name, color=None, resolution=1.0):
1619 color = self.color_dict[name]
1623 self.hier_geometry_pairs[name] = []
1624 protein_h = self.hier_dict[name]
1626 pbr = IMP.pmi1.tools.sort_by_residues(pbr)
1628 for n
in range(len(pbr) - 1):
1629 self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1632 self, name, resolution=10, scale=1.0):
1634 Generate restraints between contiguous fragments in the hierarchy.
1635 The linkers are generated at resolution 10 by default.
1642 protein_h = self.hier_dict[name]
1644 frs = self.hier_db.get_preroot_fragments_by_resolution(
1650 start = fr.get_children()[0]
1655 end = fr.get_children()[-1]
1661 SortedSegments.append((start, end, startres))
1662 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1665 for x
in range(len(SortedSegments) - 1):
1666 last = SortedSegments[x][1]
1667 first = SortedSegments[x + 1][0]
1674 residuegap = firstresn - lastresn - 1
1675 if self.disorderedlength
and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1678 optdist = sqrt(5 / 3) * 1.93 * \
1679 (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1681 if self.upperharmonic:
1687 optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1688 if self.upperharmonic:
1694 pt0 = last.get_particle()
1695 pt1 = first.get_particle()
1698 print(
"Adding sequence connectivity restraint between", pt0.get_name(),
" and ", pt1.get_name(),
'of distance', optdist)
1699 sortedsegments_cr.add_restraint(r)
1700 self.linker_restraints_dict[
1701 "LinkerRestraint-" + pt0.get_name() +
"-" + pt1.get_name()] = r
1702 self.connected_intra_pairs.append((pt0, pt1))
1703 self.connected_intra_pairs.append((pt1, pt0))
1705 self.linker_restraints.add_restraint(sortedsegments_cr)
1706 self.linker_restraints.add_restraint(unmodeledregions_cr)
1708 self.sortedsegments_cr_dict[name] = sortedsegments_cr
1709 self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1711 def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1713 pts = IMP.pmi1.tools.ParticleToSampleList()
1714 for n, fb
in enumerate(self.floppy_bodies):
1715 pts.add_particle(fb,
"Floppy_Bodies", 1.0,
"Floppy_Body_" + str(n))
1716 if len(pts.get_particles_to_sample()) > 0:
1718 print(
"optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1721 print(
"optimize_floppy_bodies: no particle to optimize")
1724 nSymmetry=
None, skip_gaussian_in_clones=
False):
1726 The copies must not contain rigid bodies.
1727 The symmetry restraints are applied at each leaf.
1731 self.representation_is_modified =
True
1732 ncopies = len(copies) + 1
1735 for k
in range(len(copies)):
1736 if (nSymmetry
is None):
1737 rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1740 rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1742 rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1750 for n, p
in enumerate(main_hiers):
1751 if (skip_gaussian_in_clones):
1760 self.model.add_score_state(c)
1761 print(
"Completed setting " + str(maincopy) +
" as a reference for "
1762 + str(copies[k]) +
" by rotating it by "
1763 + str(rotation_angle / 2.0 / pi * 360)
1764 +
" degrees around the " + str(rotational_axis) +
" axis.")
1767 def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1770 self.representation_is_modified =
True
1772 mainparticles = particles_reference
1774 t=initial_transformation
1776 p.set_name(
"RigidBody_Symmetry")
1782 copyparticles = particles_copy
1786 for n, p
in enumerate(mainparticles):
1788 pc = copyparticles[n]
1790 mainpurged.append(p)
1796 copypurged.append(pc)
1803 for n, p
in enumerate(mainpurged):
1806 print(
"setting " + p.get_name() +
" as reference for " + pc.get_name())
1809 lc.add(pc.get_index())
1812 self.model.add_score_state(c)
1815 self.rigid_bodies.append(rb)
1816 self.rigid_body_symmetries.append(rb)
1817 rb.set_name(label+
".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1820 def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1821 longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1824 self.representation_is_modified =
True
1827 protein_h = self.hier_dict[maincopy]
1830 for ilc
in range(-longitudinal_copies, longitudinal_copies + 1):
1831 for iac
in range(axial_copies):
1832 copyname = maincopy +
"_a" + str(ilc) +
"_l" + str(iac)
1833 self.create_component(copyname, 0.0)
1834 for hier
in protein_h.get_children():
1836 copyhier = self.hier_dict[copyname]
1837 outhiers.append(copyhier)
1841 2 * pi / axial_copies * (float(iac)))
1842 translation_vector = tuple(
1843 [translation_value * float(ilc) * x
for x
in axis])
1844 print(translation_vector)
1849 for n, p
in enumerate(mainparts):
1856 lc.add(pc.get_index())
1858 self.model.add_score_state(c)
1864 Load coordinates in the current representation.
1865 This should be done only if the hierarchy self.prot is identical
1866 to the one as stored in the rmf i.e. all components were added.
1871 rh = RMF.open_rmf_file_read_only(rmfname)
1879 create the representation (i.e. hierarchies) from the rmf file.
1880 it will be stored in self.prot, which will be overwritten.
1881 load the coordinates from the rmf file at frameindex.
1883 rh = RMF.open_rmf_file_read_only(rmfname)
1889 for p
in self.prot.get_children():
1890 self.create_component(p.get_name())
1891 self.hier_dict[p.get_name()] = p
1893 still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1894 save floppy bodies in self.floppy_bodies
1895 get the connectivity restraints
1900 Construct a rigid body from hierarchies (and optionally particles).
1902 @param hiers list of hierarchies to make rigid
1903 @param particles list of particles to add to the rigid body
1906 if particles
is None:
1909 rigid_parts = set(particles)
1912 print(
"set_rigid_body_from_hierarchies> setting up a new rigid body")
1919 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1922 name += hier.get_name() +
"-"
1923 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1925 if len(list(rigid_parts)) != 0:
1927 rb.set_coordinates_are_optimized(
True)
1928 rb.set_name(name +
"rigid_body")
1929 self.rigid_bodies.append(rb)
1933 print(
"set_rigid_body_from_hierarchies> rigid body could not be setup")
1937 Construct a rigid body from a list of subunits.
1939 Each subunit is a tuple that identifies the residue ranges and the
1940 component name (as used in create_component()).
1942 subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1944 [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1946 example: ["prot1","prot2",("prot3",(1,10))]
1948 sometimes, we know about structure of an interaction
1949 and here we make such PPIs rigid
1954 if type(s) == type(tuple())
and len(s) == 2:
1958 residue_indexes=list(range(s[1][0],
1960 if len(sel.get_selected_particles()) == 0:
1961 print(
"set_rigid_bodies: selected particle does not exist")
1962 for p
in sel.get_selected_particles():
1966 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1970 elif type(s) == type(str()):
1972 if len(sel.get_selected_particles()) == 0:
1973 print(
"set_rigid_bodies: selected particle does not exist")
1974 for p
in sel.get_selected_particles():
1978 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1983 rb.set_coordinates_are_optimized(
True)
1984 rb.set_name(
''.join(str(subunits)) +
"_rigid_body")
1985 self.rigid_bodies.append(rb)
1988 def set_super_rigid_body_from_hierarchies(
1994 super_rigid_xyzs = set()
1995 super_rigid_rbs = set()
1997 print(
"set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
2004 super_rigid_rbs.add(rb)
2006 super_rigid_xyzs.add(p)
2007 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
2008 if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
2011 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2014 self.super_rigid_bodies.append(
2015 (super_rigid_xyzs, super_rigid_rbs, axis))
2017 def fix_rigid_bodies(self, rigid_bodies):
2018 self.fixed_rigid_bodies += rigid_bodies
2022 self, hiers, min_length=
None, max_length=
None, axis=
None):
2024 Make a chain of super rigid bodies from a list of hierarchies.
2026 Takes a linear list of hierarchies (they are supposed to be
2027 sequence-contiguous) and produces a chain of super rigid bodies
2028 with given length range, specified by min_length and max_length.
2031 hiers = IMP.pmi1.tools.flatten_list(hiers)
2035 self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2037 def set_super_rigid_bodies(self, subunits, coords=None):
2038 super_rigid_xyzs = set()
2039 super_rigid_rbs = set()
2042 if type(s) == type(tuple())
and len(s) == 3:
2046 residue_indexes=list(range(s[0],
2048 if len(sel.get_selected_particles()) == 0:
2049 print(
"set_rigid_bodies: selected particle does not exist")
2050 for p
in sel.get_selected_particles():
2053 super_rigid_rbs.add(rb)
2055 super_rigid_xyzs.add(p)
2056 elif type(s) == type(str()):
2058 if len(sel.get_selected_particles()) == 0:
2059 print(
"set_rigid_bodies: selected particle does not exist")
2060 for p
in sel.get_selected_particles():
2064 super_rigid_rbs.add(rb)
2066 super_rigid_xyzs.add(p)
2067 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2071 Remove leaves of hierarchies from the floppy bodies list based
2072 on the component name
2075 ps=[h.get_particle()
for h
in hs]
2076 for p
in self.floppy_bodies:
2078 if p
in ps: self.floppy_bodies.remove(p)
2079 if p
in hs: self.floppy_bodies.remove(p)
2086 Remove leaves of hierarchies from the floppy bodies list.
2088 Given a list of hierarchies, get the leaves and remove the
2089 corresponding particles from the floppy bodies list. We need this
2090 function because sometimes
2091 we want to constrain the floppy bodies in a rigid body - for instance
2092 when you want to associate a bead with a density particle.
2094 for h
in hierarchies:
2097 if p
in self.floppy_bodies:
2098 print(
"remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2099 self.floppy_bodies.remove(p)
2102 def set_floppy_bodies(self):
2103 for p
in self.floppy_bodies:
2105 p.set_name(name +
"_floppy_body")
2107 print(
"I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2110 rb.set_is_rigid_member(p.get_index(),
False)
2113 rb.set_is_rigid_member(p.get_particle_index(),
False)
2114 p.set_name(p.get_name() +
"_rigid_body_member")
2116 def set_floppy_bodies_from_hierarchies(self, hiers):
2121 self.floppy_bodies.append(p)
2128 selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2129 @return the particles
2132 print(selection_tuples)
2133 for s
in selection_tuples:
2134 ps = IMP.pmi1.tools.select_by_tuple(
2135 representation=self, tupleselection=tuple(s),
2136 resolution=
None, name_is_ambiguous=
False)
2140 def get_connected_intra_pairs(self):
2141 return self.connected_intra_pairs
2143 def set_rigid_bodies_max_trans(self, maxtrans):
2144 self.maxtrans_rb = maxtrans
2146 def set_rigid_bodies_max_rot(self, maxrot):
2147 self.maxrot_rb = maxrot
2149 def set_super_rigid_bodies_max_trans(self, maxtrans):
2150 self.maxtrans_srb = maxtrans
2152 def set_super_rigid_bodies_max_rot(self, maxrot):
2153 self.maxrot_srb = maxrot
2155 def set_floppy_bodies_max_trans(self, maxtrans):
2156 self.maxtrans_fb = maxtrans
2160 Fix rigid bodies in their actual position.
2161 The get_particles_to_sample() function will return
2162 just the floppy bodies (if they are not fixed).
2164 self.rigidbodiesarefixed = rigidbodiesarefixed
2168 Fix floppy bodies in their actual position.
2169 The get_particles_to_sample() function will return
2170 just the rigid bodies (if they are not fixed).
2172 self.floppybodiesarefixed=floppybodiesarefixed
2174 def draw_hierarchy_graph(self):
2176 print(
"Drawing hierarchy graph for " + c.get_name())
2177 IMP.pmi1.output.get_graph_from_hierarchy(c)
2179 def get_geometries(self):
2182 for name
in self.hier_geometry_pairs:
2183 for pt
in self.hier_geometry_pairs[name]:
2197 def setup_bonds(self):
2200 for name
in self.hier_geometry_pairs:
2201 for pt
in self.hier_geometry_pairs[name]:
2216 def show_component_table(self, name):
2217 if name
in self.sequence_dict:
2218 lastresn = len(self.sequence_dict[name])
2221 residues = self.hier_db.get_residue_numbers(name)
2222 firstresn = min(residues)
2223 lastresn = max(residues)
2225 for nres
in range(firstresn, lastresn + 1):
2227 resolutions = self.hier_db.get_residue_resolutions(name, nres)
2229 for r
in resolutions:
2230 ps += self.hier_db.get_particles(name, nres, r)
2231 print(
"%20s %7s" % (name, nres),
" ".join([
"%20s %7s" % (str(p.get_name()),
2234 print(
"%20s %20s" % (name, nres),
"**** not represented ****")
2236 def draw_hierarchy_composition(self):
2238 ks = list(self.elements.keys())
2242 for k
in self.elements:
2243 for l
in self.elements[k]:
2248 self.draw_component_composition(k, max)
2250 def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2251 from matplotlib
import pyplot
2252 import matplotlib
as mpl
2254 tmplist = sorted(self.elements[k], key=itemgetter(0))
2256 endres = tmplist[-1][1]
2258 print(
"draw_component_composition: missing information for component %s" % name)
2260 fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2261 ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2266 norm = mpl.colors.Normalize(vmin=5, vmax=10)
2270 for n, l
in enumerate(tmplist):
2274 if bounds[-1] != l[0]:
2275 colors.append(
"white")
2278 colors.append(
"#99CCFF")
2280 colors.append(
"#FFFF99")
2282 colors.append(
"#33CCCC")
2284 bounds.append(l[1] + 1)
2287 colors.append(
"#99CCFF")
2289 colors.append(
"#FFFF99")
2291 colors.append(
"#33CCCC")
2293 bounds.append(l[1] + 1)
2295 if bounds[-1] - 1 == l[0]:
2299 colors.append(
"white")
2302 bounds.append(bounds[-1])
2303 colors.append(
"white")
2304 cmap = mpl.colors.ListedColormap(colors)
2305 cmap.set_over(
'0.25')
2306 cmap.set_under(
'0.75')
2308 norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2309 cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2315 spacing=
'proportional',
2316 orientation=
'horizontal')
2325 mid = 1.0 / endres * float(l[0])
2330 l[2], xy=(mid, 1), xycoords=
'axes fraction',
2331 xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords=
'axes fraction',
2332 arrowprops=dict(arrowstyle=
"->",
2333 connectionstyle=
"angle,angleA=0,angleB=90,rad=10"),
2335 extra_artists.append(t)
2338 title = ax.text(-0.005, 0.5, k, ha=
"right", va=
"center", rotation=90,
2341 extra_artists.append(title)
2343 labels = len(bounds) * [
" "]
2344 ax.set_xticklabels(labels)
2345 mid = 1.0 / endres * float(bounds[0])
2346 t = ax.annotate(bounds[0], xy=(mid, 0), xycoords=
'axes fraction',
2347 xytext=(mid - 0.01, -0.5), textcoords=
'axes fraction',)
2348 extra_artists.append(t)
2349 offsets = [0, -0.5, -1.0]
2351 for n
in range(1, len(bounds)):
2352 if bounds[n] == bounds[n - 1]:
2354 mid = 1.0 / endres * float(bounds[n])
2355 if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2357 offset = offsets[nclashes % 3]
2363 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2364 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction')
2367 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2368 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction', arrowprops=dict(arrowstyle=
"-"))
2369 extra_artists.append(t)
2371 cb2.add_lines(bounds, [
"black"] * len(bounds), [1] * len(bounds))
2375 k +
"structure.pdf",
2378 bbox_extra_artists=(extra_artists),
2379 bbox_inches=
'tight')
2382 def draw_coordinates_projection(self):
2383 import matplotlib.pyplot
as pp
2386 for name
in self.hier_geometry_pairs:
2387 for pt
in self.hier_geometry_pairs[name]:
2397 xpairs.append([x1, x2])
2398 ypairs.append([y1, y2])
2401 for xends, yends
in zip(xpairs, ypairs):
2406 pp.plot(xlist, ylist,
'b-', alpha=0.1)
2409 def get_prot_name_from_particle(self, particle):
2410 names = self.get_component_names()
2411 particle0 = particle
2413 while not name
in names:
2416 particle0 = h.get_particle()
2420 def get_particles_to_sample(self):
2430 if not self.rigidbodiesarefixed:
2431 for rb
in self.rigid_bodies:
2436 if rb
not in self.fixed_rigid_bodies:
2439 if not self.floppybodiesarefixed:
2440 for fb
in self.floppy_bodies:
2447 for srb
in self.super_rigid_bodies:
2450 rigid_bodies = list(srb[1])
2451 filtered_rigid_bodies = []
2452 for rb
in rigid_bodies:
2453 if rb
not in self.fixed_rigid_bodies:
2454 filtered_rigid_bodies.append(rb)
2455 srbtmp.append((srb[0], filtered_rigid_bodies))
2457 self.rigid_bodies = rbtmp
2458 self.floppy_bodies = fbtmp
2459 self.super_rigid_bodies = srbtmp
2461 ps[
"Rigid_Bodies_SimplifiedModel"] = (
2465 ps[
"Floppy_Bodies_SimplifiedModel"] = (
2468 ps[
"SR_Bodies_SimplifiedModel"] = (
2469 self.super_rigid_bodies,
2474 def set_output_level(self, level):
2475 self.output_level = level
2477 def _evaluate(self, deriv):
2478 """Evaluate the total score of all added restraints"""
2480 return r.evaluate(deriv)
2482 def get_output(self):
2486 output[
"SimplifiedModel_Total_Score_" +
2487 self.label] = str(self._evaluate(
False))
2488 output[
"SimplifiedModel_Linker_Score_" +
2489 self.label] = str(self.linker_restraints.unprotected_evaluate(
None))
2490 for name
in self.sortedsegments_cr_dict:
2491 partialscore = self.sortedsegments_cr_dict[name].evaluate(
False)
2492 score += partialscore
2494 "SimplifiedModel_Link_SortedSegments_" +
2499 partialscore = self.unmodeledregions_cr_dict[name].evaluate(
False)
2500 score += partialscore
2502 "SimplifiedModel_Link_UnmodeledRegions_" +
2507 for rb
in self.rigid_body_symmetries:
2509 output[name +
"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2510 for name
in self.linker_restraints_dict:
2515 self.linker_restraints_dict[
2516 name].unprotected_evaluate(
2519 if len(self.reference_structures.keys()) != 0:
2520 rmsds = self.get_all_rmsds()
2523 "SimplifiedModel_" +
2526 self.label] = rmsds[
2529 if self.output_level ==
"high":
2532 output[
"Coordinates_" +
2533 p.get_name() +
"_" + self.label] = str(d)
2535 output[
"_TotalScore"] = str(score)
2538 def get_test_output(self):
2540 output = self.get_output()
2541 for n, p
in enumerate(self.get_particles_to_sample()):
2542 output[
"Particle_to_sample_" + str(n)] = str(p)
2544 output[
"Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2545 output[
"Number_of_floppy_bodies"] = len(self.floppy_bodies)
2546 output[
"Number_of_rigid_bodies"] = len(self.rigid_bodies)
2547 output[
"Number_of_super_bodies"] = len(self.super_rigid_bodies)
2548 output[
"Selection_resolution_1"] = len(
2550 output[
"Selection_resolution_5"] = len(
2552 output[
"Selection_resolution_7"] = len(
2554 output[
"Selection_resolution_10"] = len(
2556 output[
"Selection_resolution_100"] = len(
2559 output[
"Selection_resolution=1"] = len(
2561 output[
"Selection_resolution=1,resid=10"] = len(
2563 for resolution
in self.hier_resolution:
2564 output[
"Hier_resolution_dictionary" +
2565 str(resolution)] = len(self.hier_resolution[resolution])
2566 for name
in self.hier_dict:
2568 "Selection_resolution=1,resid=10,name=" +
2576 "Selection_resolution=1,resid=10,name=" +
2578 ",ambiguous"] = len(
2583 name_is_ambiguous=
True,
2586 "Selection_resolution=1,resid=10,name=" +
2588 ",ambiguous"] = len(
2593 name_is_ambiguous=
True,
2596 "Selection_resolution=1,resrange=(10,20),name=" +
2605 "Selection_resolution=1,resrange=(10,20),name=" +
2607 ",ambiguous"] = len(
2612 name_is_ambiguous=
True,
2616 "Selection_resolution=10,resrange=(10,20),name=" +
2625 "Selection_resolution=10,resrange=(10,20),name=" +
2627 ",ambiguous"] = len(
2632 name_is_ambiguous=
True,
2636 "Selection_resolution=100,resrange=(10,20),name=" +
2645 "Selection_resolution=100,resrange=(10,20),name=" +
2647 ",ambiguous"] = len(
2652 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.
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.