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 = []
1445 if bounding_box
is not None:
1446 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1451 for h
in hierarchies_excluded_from_collision:
1455 allparticleindexes = list(
1456 set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1458 print(hierarchies_excluded_from_collision)
1459 print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1461 print(
'shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies),
'rigid bodies')
1462 for rb
in self.rigid_bodies:
1463 if rb
not in excluded_rigid_bodies:
1465 rbindexes = rb.get_member_particle_indexes()
1467 set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1468 otherparticleindexes = list(
1469 set(allparticleindexes) - set(rbindexes))
1471 if len(otherparticleindexes)
is None:
1475 while niter < niterations:
1476 if (ignore_initial_coordinates):
1482 if bounding_box
is not None:
1498 gcpf.get_close_pairs(
1500 otherparticleindexes,
1504 if (ignore_initial_coordinates):
1505 print (rb.get_name(),
IMP.core.XYZ(rb).get_coordinates())
1508 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1509 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
1510 if niter == niterations:
1511 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1515 print(
'shuffling', len(self.floppy_bodies),
'floppy bodies')
1516 for fb
in self.floppy_bodies:
1517 if (avoidcollision):
1522 otherparticleindexes = list(
1523 set(allparticleindexes) - set(fbindexes))
1524 if len(otherparticleindexes)
is None:
1542 while niter < niterations:
1543 if (ignore_initial_coordinates):
1549 if bounding_box
is not None:
1561 if (avoidcollision):
1564 gcpf.get_close_pairs(
1566 otherparticleindexes,
1570 if (ignore_initial_coordinates):
1571 print (fb.get_name(),
IMP.core.XYZ(fb).get_coordinates())
1574 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1575 print(
"shuffle_configuration: floppy body name: " + fb.get_name())
1576 if niter == niterations:
1577 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1581 def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1585 self.reference_structures[label] = (
1589 def get_all_rmsds(self):
1592 for label
in self.reference_structures:
1594 current_coordinates = [
IMP.core.XYZ(p).get_coordinates()
1595 for p
in self.reference_structures[label][1]]
1596 reference_coordinates = self.reference_structures[label][0]
1597 if len(reference_coordinates) != len(current_coordinates):
1598 print(
"calculate_all_rmsds: reference and actual coordinates are not the same")
1601 current_coordinates,
1602 reference_coordinates)
1604 reference_coordinates,
1605 current_coordinates)
1609 reference_coordinates,
1610 current_coordinates)
1611 rmsds[label +
"_GlobalRMSD"] = rmsd_global
1612 rmsds[label +
"_RelativeDRMS"] = rmsd_relative
1615 def setup_component_geometry(self, name, color=None, resolution=1.0):
1617 color = self.color_dict[name]
1621 self.hier_geometry_pairs[name] = []
1622 protein_h = self.hier_dict[name]
1624 pbr = IMP.pmi1.tools.sort_by_residues(pbr)
1626 for n
in range(len(pbr) - 1):
1627 self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1630 self, name, resolution=10, scale=1.0):
1632 Generate restraints between contiguous fragments in the hierarchy.
1633 The linkers are generated at resolution 10 by default.
1640 protein_h = self.hier_dict[name]
1642 frs = self.hier_db.get_preroot_fragments_by_resolution(
1648 start = fr.get_children()[0]
1653 end = fr.get_children()[-1]
1659 SortedSegments.append((start, end, startres))
1660 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1663 for x
in range(len(SortedSegments) - 1):
1664 last = SortedSegments[x][1]
1665 first = SortedSegments[x + 1][0]
1672 residuegap = firstresn - lastresn - 1
1673 if self.disorderedlength
and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1676 optdist = sqrt(5 / 3) * 1.93 * \
1677 (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1679 if self.upperharmonic:
1685 optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1686 if self.upperharmonic:
1692 pt0 = last.get_particle()
1693 pt1 = first.get_particle()
1696 print(
"Adding sequence connectivity restraint between", pt0.get_name(),
" and ", pt1.get_name(),
'of distance', optdist)
1697 sortedsegments_cr.add_restraint(r)
1698 self.linker_restraints_dict[
1699 "LinkerRestraint-" + pt0.get_name() +
"-" + pt1.get_name()] = r
1700 self.connected_intra_pairs.append((pt0, pt1))
1701 self.connected_intra_pairs.append((pt1, pt0))
1703 self.linker_restraints.add_restraint(sortedsegments_cr)
1704 self.linker_restraints.add_restraint(unmodeledregions_cr)
1706 self.sortedsegments_cr_dict[name] = sortedsegments_cr
1707 self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1709 def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1711 pts = IMP.pmi1.tools.ParticleToSampleList()
1712 for n, fb
in enumerate(self.floppy_bodies):
1713 pts.add_particle(fb,
"Floppy_Bodies", 1.0,
"Floppy_Body_" + str(n))
1714 if len(pts.get_particles_to_sample()) > 0:
1716 print(
"optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1719 print(
"optimize_floppy_bodies: no particle to optimize")
1722 nSymmetry=
None, skip_gaussian_in_clones=
False):
1724 The copies must not contain rigid bodies.
1725 The symmetry restraints are applied at each leaf.
1729 self.representation_is_modified =
True
1730 ncopies = len(copies) + 1
1733 for k
in range(len(copies)):
1734 if (nSymmetry
is None):
1735 rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1738 rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1740 rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1748 for n, p
in enumerate(main_hiers):
1749 if (skip_gaussian_in_clones):
1758 self.model.add_score_state(c)
1759 print(
"Completed setting " + str(maincopy) +
" as a reference for "
1760 + str(copies[k]) +
" by rotating it by "
1761 + str(rotation_angle / 2.0 / pi * 360)
1762 +
" degrees around the " + str(rotational_axis) +
" axis.")
1765 def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1768 self.representation_is_modified =
True
1770 mainparticles = particles_reference
1772 t=initial_transformation
1774 p.set_name(
"RigidBody_Symmetry")
1780 copyparticles = particles_copy
1784 for n, p
in enumerate(mainparticles):
1786 pc = copyparticles[n]
1788 mainpurged.append(p)
1794 copypurged.append(pc)
1801 for n, p
in enumerate(mainpurged):
1804 print(
"setting " + p.get_name() +
" as reference for " + pc.get_name())
1807 lc.add(pc.get_index())
1810 self.model.add_score_state(c)
1813 self.rigid_bodies.append(rb)
1814 self.rigid_body_symmetries.append(rb)
1815 rb.set_name(label+
".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1818 def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1819 longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1822 self.representation_is_modified =
True
1825 protein_h = self.hier_dict[maincopy]
1828 for ilc
in range(-longitudinal_copies, longitudinal_copies + 1):
1829 for iac
in range(axial_copies):
1830 copyname = maincopy +
"_a" + str(ilc) +
"_l" + str(iac)
1831 self.create_component(copyname, 0.0)
1832 for hier
in protein_h.get_children():
1834 copyhier = self.hier_dict[copyname]
1835 outhiers.append(copyhier)
1839 2 * pi / axial_copies * (float(iac)))
1840 translation_vector = tuple(
1841 [translation_value * float(ilc) * x
for x
in axis])
1842 print(translation_vector)
1847 for n, p
in enumerate(mainparts):
1854 lc.add(pc.get_index())
1856 self.model.add_score_state(c)
1862 Load coordinates in the current representation.
1863 This should be done only if the hierarchy self.prot is identical
1864 to the one as stored in the rmf i.e. all components were added.
1869 rh = RMF.open_rmf_file_read_only(rmfname)
1877 create the representation (i.e. hierarchies) from the rmf file.
1878 it will be stored in self.prot, which will be overwritten.
1879 load the coordinates from the rmf file at frameindex.
1881 rh = RMF.open_rmf_file_read_only(rmfname)
1887 for p
in self.prot.get_children():
1888 self.create_component(p.get_name())
1889 self.hier_dict[p.get_name()] = p
1891 still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1892 save floppy bodies in self.floppy_bodies
1893 get the connectivity restraints
1898 Construct a rigid body from hierarchies (and optionally particles).
1900 @param hiers list of hierarchies to make rigid
1901 @param particles list of particles to add to the rigid body
1904 if particles
is None:
1907 rigid_parts = set(particles)
1910 print(
"set_rigid_body_from_hierarchies> setting up a new rigid body")
1917 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1920 name += hier.get_name() +
"-"
1921 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1923 if len(list(rigid_parts)) != 0:
1925 rb.set_coordinates_are_optimized(
True)
1926 rb.set_name(name +
"rigid_body")
1927 self.rigid_bodies.append(rb)
1931 print(
"set_rigid_body_from_hierarchies> rigid body could not be setup")
1935 Construct a rigid body from a list of subunits.
1937 Each subunit is a tuple that identifies the residue ranges and the
1938 component name (as used in create_component()).
1940 subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1942 [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1944 example: ["prot1","prot2",("prot3",(1,10))]
1946 sometimes, we know about structure of an interaction
1947 and here we make such PPIs rigid
1952 if type(s) == type(tuple())
and len(s) == 2:
1956 residue_indexes=list(range(s[1][0],
1958 if len(sel.get_selected_particles()) == 0:
1959 print(
"set_rigid_bodies: selected particle does not exist")
1960 for p
in sel.get_selected_particles():
1964 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1968 elif type(s) == type(str()):
1970 if len(sel.get_selected_particles()) == 0:
1971 print(
"set_rigid_bodies: selected particle does not exist")
1972 for p
in sel.get_selected_particles():
1976 print(
"set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1981 rb.set_coordinates_are_optimized(
True)
1982 rb.set_name(
''.join(str(subunits)) +
"_rigid_body")
1983 self.rigid_bodies.append(rb)
1986 def set_super_rigid_body_from_hierarchies(
1992 super_rigid_xyzs = set()
1993 super_rigid_rbs = set()
1995 print(
"set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
2002 super_rigid_rbs.add(rb)
2004 super_rigid_xyzs.add(p)
2005 print(
"set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
2006 if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
2009 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2012 self.super_rigid_bodies.append(
2013 (super_rigid_xyzs, super_rigid_rbs, axis))
2015 def fix_rigid_bodies(self, rigid_bodies):
2016 self.fixed_rigid_bodies += rigid_bodies
2020 self, hiers, min_length=
None, max_length=
None, axis=
None):
2022 Make a chain of super rigid bodies from a list of hierarchies.
2024 Takes a linear list of hierarchies (they are supposed to be
2025 sequence-contiguous) and produces a chain of super rigid bodies
2026 with given length range, specified by min_length and max_length.
2029 hiers = IMP.pmi1.tools.flatten_list(hiers)
2033 self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2035 def set_super_rigid_bodies(self, subunits, coords=None):
2036 super_rigid_xyzs = set()
2037 super_rigid_rbs = set()
2040 if type(s) == type(tuple())
and len(s) == 3:
2044 residue_indexes=list(range(s[0],
2046 if len(sel.get_selected_particles()) == 0:
2047 print(
"set_rigid_bodies: selected particle does not exist")
2048 for p
in sel.get_selected_particles():
2051 super_rigid_rbs.add(rb)
2053 super_rigid_xyzs.add(p)
2054 elif type(s) == type(str()):
2056 if len(sel.get_selected_particles()) == 0:
2057 print(
"set_rigid_bodies: selected particle does not exist")
2058 for p
in sel.get_selected_particles():
2062 super_rigid_rbs.add(rb)
2064 super_rigid_xyzs.add(p)
2065 self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2069 Remove leaves of hierarchies from the floppy bodies list based
2070 on the component name
2073 ps=[h.get_particle()
for h
in hs]
2074 for p
in self.floppy_bodies:
2076 if p
in ps: self.floppy_bodies.remove(p)
2077 if p
in hs: self.floppy_bodies.remove(p)
2084 Remove leaves of hierarchies from the floppy bodies list.
2086 Given a list of hierarchies, get the leaves and remove the
2087 corresponding particles from the floppy bodies list. We need this
2088 function because sometimes
2089 we want to constrain the floppy bodies in a rigid body - for instance
2090 when you want to associate a bead with a density particle.
2092 for h
in hierarchies:
2095 if p
in self.floppy_bodies:
2096 print(
"remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2097 self.floppy_bodies.remove(p)
2100 def set_floppy_bodies(self):
2101 for p
in self.floppy_bodies:
2103 p.set_name(name +
"_floppy_body")
2105 print(
"I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2108 rb.set_is_rigid_member(p.get_index(),
False)
2111 rb.set_is_rigid_member(p.get_particle_index(),
False)
2112 p.set_name(p.get_name() +
"_rigid_body_member")
2114 def set_floppy_bodies_from_hierarchies(self, hiers):
2119 self.floppy_bodies.append(p)
2126 selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2127 @return the particles
2130 print(selection_tuples)
2131 for s
in selection_tuples:
2132 ps = IMP.pmi1.tools.select_by_tuple(
2133 representation=self, tupleselection=tuple(s),
2134 resolution=
None, name_is_ambiguous=
False)
2138 def get_connected_intra_pairs(self):
2139 return self.connected_intra_pairs
2141 def set_rigid_bodies_max_trans(self, maxtrans):
2142 self.maxtrans_rb = maxtrans
2144 def set_rigid_bodies_max_rot(self, maxrot):
2145 self.maxrot_rb = maxrot
2147 def set_super_rigid_bodies_max_trans(self, maxtrans):
2148 self.maxtrans_srb = maxtrans
2150 def set_super_rigid_bodies_max_rot(self, maxrot):
2151 self.maxrot_srb = maxrot
2153 def set_floppy_bodies_max_trans(self, maxtrans):
2154 self.maxtrans_fb = maxtrans
2158 Fix rigid bodies in their actual position.
2159 The get_particles_to_sample() function will return
2160 just the floppy bodies (if they are not fixed).
2162 self.rigidbodiesarefixed = rigidbodiesarefixed
2166 Fix floppy bodies in their actual position.
2167 The get_particles_to_sample() function will return
2168 just the rigid bodies (if they are not fixed).
2170 self.floppybodiesarefixed=floppybodiesarefixed
2172 def draw_hierarchy_graph(self):
2174 print(
"Drawing hierarchy graph for " + c.get_name())
2175 IMP.pmi1.output.get_graph_from_hierarchy(c)
2177 def get_geometries(self):
2180 for name
in self.hier_geometry_pairs:
2181 for pt
in self.hier_geometry_pairs[name]:
2195 def setup_bonds(self):
2198 for name
in self.hier_geometry_pairs:
2199 for pt
in self.hier_geometry_pairs[name]:
2214 def show_component_table(self, name):
2215 if name
in self.sequence_dict:
2216 lastresn = len(self.sequence_dict[name])
2219 residues = self.hier_db.get_residue_numbers(name)
2220 firstresn = min(residues)
2221 lastresn = max(residues)
2223 for nres
in range(firstresn, lastresn + 1):
2225 resolutions = self.hier_db.get_residue_resolutions(name, nres)
2227 for r
in resolutions:
2228 ps += self.hier_db.get_particles(name, nres, r)
2229 print(
"%20s %7s" % (name, nres),
" ".join([
"%20s %7s" % (str(p.get_name()),
2232 print(
"%20s %20s" % (name, nres),
"**** not represented ****")
2234 def draw_hierarchy_composition(self):
2236 ks = list(self.elements.keys())
2240 for k
in self.elements:
2241 for l
in self.elements[k]:
2246 self.draw_component_composition(k, max)
2248 def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2249 from matplotlib
import pyplot
2250 import matplotlib
as mpl
2252 tmplist = sorted(self.elements[k], key=itemgetter(0))
2254 endres = tmplist[-1][1]
2256 print(
"draw_component_composition: missing information for component %s" % name)
2258 fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2259 ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2264 norm = mpl.colors.Normalize(vmin=5, vmax=10)
2268 for n, l
in enumerate(tmplist):
2272 if bounds[-1] != l[0]:
2273 colors.append(
"white")
2276 colors.append(
"#99CCFF")
2278 colors.append(
"#FFFF99")
2280 colors.append(
"#33CCCC")
2282 bounds.append(l[1] + 1)
2285 colors.append(
"#99CCFF")
2287 colors.append(
"#FFFF99")
2289 colors.append(
"#33CCCC")
2291 bounds.append(l[1] + 1)
2293 if bounds[-1] - 1 == l[0]:
2297 colors.append(
"white")
2300 bounds.append(bounds[-1])
2301 colors.append(
"white")
2302 cmap = mpl.colors.ListedColormap(colors)
2303 cmap.set_over(
'0.25')
2304 cmap.set_under(
'0.75')
2306 norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2307 cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2313 spacing=
'proportional',
2314 orientation=
'horizontal')
2323 mid = 1.0 / endres * float(l[0])
2328 l[2], xy=(mid, 1), xycoords=
'axes fraction',
2329 xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords=
'axes fraction',
2330 arrowprops=dict(arrowstyle=
"->",
2331 connectionstyle=
"angle,angleA=0,angleB=90,rad=10"),
2333 extra_artists.append(t)
2336 title = ax.text(-0.005, 0.5, k, ha=
"right", va=
"center", rotation=90,
2339 extra_artists.append(title)
2341 labels = len(bounds) * [
" "]
2342 ax.set_xticklabels(labels)
2343 mid = 1.0 / endres * float(bounds[0])
2344 t = ax.annotate(bounds[0], xy=(mid, 0), xycoords=
'axes fraction',
2345 xytext=(mid - 0.01, -0.5), textcoords=
'axes fraction',)
2346 extra_artists.append(t)
2347 offsets = [0, -0.5, -1.0]
2349 for n
in range(1, len(bounds)):
2350 if bounds[n] == bounds[n - 1]:
2352 mid = 1.0 / endres * float(bounds[n])
2353 if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2355 offset = offsets[nclashes % 3]
2361 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2362 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction')
2365 bounds[n], xy=(mid, 0), xycoords=
'axes fraction',
2366 xytext=(mid, -0.5 + offset), textcoords=
'axes fraction', arrowprops=dict(arrowstyle=
"-"))
2367 extra_artists.append(t)
2369 cb2.add_lines(bounds, [
"black"] * len(bounds), [1] * len(bounds))
2373 k +
"structure.pdf",
2376 bbox_extra_artists=(extra_artists),
2377 bbox_inches=
'tight')
2380 def draw_coordinates_projection(self):
2381 import matplotlib.pyplot
as pp
2384 for name
in self.hier_geometry_pairs:
2385 for pt
in self.hier_geometry_pairs[name]:
2395 xpairs.append([x1, x2])
2396 ypairs.append([y1, y2])
2399 for xends, yends
in zip(xpairs, ypairs):
2404 pp.plot(xlist, ylist,
'b-', alpha=0.1)
2407 def get_prot_name_from_particle(self, particle):
2408 names = self.get_component_names()
2409 particle0 = particle
2411 while not name
in names:
2414 particle0 = h.get_particle()
2418 def get_particles_to_sample(self):
2428 if not self.rigidbodiesarefixed:
2429 for rb
in self.rigid_bodies:
2434 if rb
not in self.fixed_rigid_bodies:
2437 if not self.floppybodiesarefixed:
2438 for fb
in self.floppy_bodies:
2445 for srb
in self.super_rigid_bodies:
2448 rigid_bodies = list(srb[1])
2449 filtered_rigid_bodies = []
2450 for rb
in rigid_bodies:
2451 if rb
not in self.fixed_rigid_bodies:
2452 filtered_rigid_bodies.append(rb)
2453 srbtmp.append((srb[0], filtered_rigid_bodies))
2455 self.rigid_bodies = rbtmp
2456 self.floppy_bodies = fbtmp
2457 self.super_rigid_bodies = srbtmp
2459 ps[
"Rigid_Bodies_SimplifiedModel"] = (
2463 ps[
"Floppy_Bodies_SimplifiedModel"] = (
2466 ps[
"SR_Bodies_SimplifiedModel"] = (
2467 self.super_rigid_bodies,
2472 def set_output_level(self, level):
2473 self.output_level = level
2475 def _evaluate(self, deriv):
2476 """Evaluate the total score of all added restraints"""
2478 return r.evaluate(deriv)
2480 def get_output(self):
2484 output[
"SimplifiedModel_Total_Score_" +
2485 self.label] = str(self._evaluate(
False))
2486 output[
"SimplifiedModel_Linker_Score_" +
2487 self.label] = str(self.linker_restraints.unprotected_evaluate(
None))
2488 for name
in self.sortedsegments_cr_dict:
2489 partialscore = self.sortedsegments_cr_dict[name].evaluate(
False)
2490 score += partialscore
2492 "SimplifiedModel_Link_SortedSegments_" +
2497 partialscore = self.unmodeledregions_cr_dict[name].evaluate(
False)
2498 score += partialscore
2500 "SimplifiedModel_Link_UnmodeledRegions_" +
2505 for rb
in self.rigid_body_symmetries:
2507 output[name +
"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2508 for name
in self.linker_restraints_dict:
2513 self.linker_restraints_dict[
2514 name].unprotected_evaluate(
2517 if len(self.reference_structures.keys()) != 0:
2518 rmsds = self.get_all_rmsds()
2521 "SimplifiedModel_" +
2524 self.label] = rmsds[
2527 if self.output_level ==
"high":
2530 output[
"Coordinates_" +
2531 p.get_name() +
"_" + self.label] = str(d)
2533 output[
"_TotalScore"] = str(score)
2536 def get_test_output(self):
2538 output = self.get_output()
2539 for n, p
in enumerate(self.get_particles_to_sample()):
2540 output[
"Particle_to_sample_" + str(n)] = str(p)
2542 output[
"Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2543 output[
"Number_of_floppy_bodies"] = len(self.floppy_bodies)
2544 output[
"Number_of_rigid_bodies"] = len(self.rigid_bodies)
2545 output[
"Number_of_super_bodies"] = len(self.super_rigid_bodies)
2546 output[
"Selection_resolution_1"] = len(
2548 output[
"Selection_resolution_5"] = len(
2550 output[
"Selection_resolution_7"] = len(
2552 output[
"Selection_resolution_10"] = len(
2554 output[
"Selection_resolution_100"] = len(
2557 output[
"Selection_resolution=1"] = len(
2559 output[
"Selection_resolution=1,resid=10"] = len(
2561 for resolution
in self.hier_resolution:
2562 output[
"Hier_resolution_dictionary" +
2563 str(resolution)] = len(self.hier_resolution[resolution])
2564 for name
in self.hier_dict:
2566 "Selection_resolution=1,resid=10,name=" +
2574 "Selection_resolution=1,resid=10,name=" +
2576 ",ambiguous"] = len(
2581 name_is_ambiguous=
True,
2584 "Selection_resolution=1,resid=10,name=" +
2586 ",ambiguous"] = len(
2591 name_is_ambiguous=
True,
2594 "Selection_resolution=1,resrange=(10,20),name=" +
2603 "Selection_resolution=1,resrange=(10,20),name=" +
2605 ",ambiguous"] = len(
2610 name_is_ambiguous=
True,
2614 "Selection_resolution=10,resrange=(10,20),name=" +
2623 "Selection_resolution=10,resrange=(10,20),name=" +
2625 ",ambiguous"] = len(
2630 name_is_ambiguous=
True,
2634 "Selection_resolution=100,resrange=(10,20),name=" +
2643 "Selection_resolution=100,resrange=(10,20),name=" +
2645 ",ambiguous"] = len(
2650 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)
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.