IMP logo
IMP Reference Guide  2.10.1
The Integrative Modeling Platform
pmi/representation.py
1 #!/usr/bin/env python
2 
3 """@namespace IMP.pmi.representation
4  Representation of the system.
5 """
6 
7 from __future__ import print_function
8 import IMP
9 import IMP.core
10 import IMP.algebra
11 import IMP.atom
12 import IMP.display
13 import IMP.isd
14 import IMP.pmi
15 import IMP.pmi.tools
16 import IMP.pmi.output
17 import IMP.rmf
18 import IMP.pmi.topology
19 import RMF
20 from math import pi, sqrt
21 from operator import itemgetter
22 import os
23 import weakref
24 
25 class _Repo(object):
26  def __init__(self, doi, root):
27  self.doi = doi
28  self._root = root
29 
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)
33 
34 class _StateInfo(object):
35  """Score state-specific information about this representation."""
36  short_name = None
37  long_name = None
38 
40  "use IMP.pmi.topology.System instead")
41 class Representation(object):
42  # Authors: Peter Cimermancic, Riccardo Pellarin, Charles Greenberg
43 
44  '''
45  Set up the representation of all proteins and nucleic acid macromolecules.
46 
47  Create the molecular hierarchies, representation,
48  sequence connectivity for the various involved proteins and
49  nucleic acid macromolecules:
50 
51  Create a protein, DNA or RNA, represent it as a set of connected balls of appropriate
52  radii and number of residues, PDB at given resolution(s), or ideal helices.
53 
54  How to use the SimplifiedModel class (typical use):
55 
56  see test/test_hierarchy_contruction.py
57 
58  examples:
59 
60  1) Create a chain of helices and flexible parts
61 
62  c_1_119 =self.add_component_necklace("prot1",1,119,20)
63  c_120_131 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(120,131))
64  c_132_138 =self.add_component_beads("prot1",[(132,138)])
65  c_139_156 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(139,156))
66  c_157_174 =self.add_component_beads("prot1",[(157,174)])
67  c_175_182 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(175,182))
68  c_183_194 =self.add_component_beads("prot1",[(183,194)])
69  c_195_216 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(195,216))
70  c_217_250 =self.add_component_beads("prot1",[(217,250)])
71 
72 
73  self.set_rigid_body_from_hierarchies(c_120_131)
74  self.set_rigid_body_from_hierarchies(c_139_156)
75  self.set_rigid_body_from_hierarchies(c_175_182)
76  self.set_rigid_body_from_hierarchies(c_195_216)
77 
78  clist=[c_1_119,c_120_131,c_132_138,c_139_156,c_157_174,c_175_182,c_183_194,c_195_216,
79  c_217_250]
80 
81  self.set_chain_of_super_rigid_bodies(clist,2,3)
82 
83  self.set_super_rigid_bodies(["prot1"])
84 
85  '''
86 
87  def __init__(self, model, upperharmonic=True, disorderedlength=True):
88  """Constructor.
89  @param model the model
90  @param upperharmonic This flag uses either harmonic (False)
91  or upperharmonic (True) in the intra-pair
92  connectivity restraint.
93  @param disorderedlength This flag uses either disordered length
94  calculated for random coil peptides (True) or zero
95  surface-to-surface distance between beads (False)
96  as optimal distance for the sequence connectivity
97  restraint.
98  """
99 
100  self.state = _StateInfo()
101  self._metadata = []
102  self._file_dataset = {}
103  self._protocol_output = []
104  self._non_modeled_components = {}
105 
106  # this flag uses either harmonic (False) or upperharmonic (True)
107  # in the intra-pair connectivity restraint. Harmonic is used whe you want to
108  # remove the intra-ev term from energy calculations, e.g.:
109  # upperharmonic=False
110  # ip=self.get_connected_intra_pairs()
111  # ev.add_excluded_particle_pairs(ip)
112 
113  self.upperharmonic = upperharmonic
114  self.disorderedlength = disorderedlength
115  self.rigid_bodies = []
116  self.fixed_rigid_bodies = []
117  self.fixed_floppy_bodies = []
118  self.floppy_bodies = []
119  # self.super_rigid_bodies is a list of tuples.
120  # each tuple, corresponding to a particular super rigid body
121  # the tuple is (super_rigid_xyzs,super_rigid_rbs)
122  # where super_rigid_xyzs are the flexible xyz particles
123  # and super_rigid_rbs is the list of rigid bodies.
124  self.super_rigid_bodies = []
125  self.rigid_body_symmetries = []
126  self.output_level = "low"
127  self.label = "None"
128 
129  self.maxtrans_rb = 2.0
130  self.maxrot_rb = 0.04
131  self.maxtrans_srb = 2.0
132  self.maxrot_srb = 0.2
133  self.rigidbodiesarefixed = False
134  self.floppybodiesarefixed = False
135  self.maxtrans_fb = 3.0
136  self.resolution = 10.0
137  self.bblenght = 100.0
138  self.kappa = 100.0
139  self.model = model
140 
141  self.representation_is_modified = False
142  self.unmodeledregions_cr_dict = {}
143  self.sortedsegments_cr_dict = {}
144  self.prot = IMP.atom.Hierarchy.setup_particle(IMP.Particle(self.model))
145  self.connected_intra_pairs = []
146  self.hier_dict = {}
147  self.color_dict = {}
148  self.sequence_dict = {}
149  self.hier_geometry_pairs = {}
150  self.hier_db = IMP.pmi.tools.HierarchyDatabase()
151  # this dictionary stores the hierarchies by component name and representation type
152  # self.hier_representation[name][representation_type]
153  # where representation type is Res:X, Beads, Densities, Representation,
154  # etc...
155  self.hier_representation = {}
156  self.hier_resolution = {}
157  # reference structures is a dictionary that contains the coordinates of
158  # structures that are used to calculate the rmsd
159  self.reference_structures = {}
160  self.elements = {}
161  self.linker_restraints = IMP.RestraintSet(self.model, "linker_restraints")
162  self.linker_restraints.set_was_used(True)
163  self.linker_restraints_dict = {}
164  self.threetoone = {'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D',
165  'CYS': 'C', 'GLU': 'E', 'GLN': 'Q', 'GLY': 'G',
166  'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'LYS': 'K',
167  'MET': 'M', 'PHE': 'F', 'PRO': 'P', 'SER': 'S',
168  'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V', 'UNK': 'X'}
169 
170  self.onetothree = dict((v, k) for k, v in self.threetoone.items())
171 
172  self.residuenamekey = IMP.StringKey("ResidueName")
173 
174  @property
175  @IMP.deprecated_method("2.10", "Model should be accessed with `.model`.")
176  def m(self):
177  return self.model
178 
179  def add_metadata(self, m):
180  """Associate some metadata with this modeling.
181  @param m an instance of an ihm metadata class, such as
182  ihm.Software, ihm.Citation, or ihm.location.Repository.
183  """
184  self._metadata.append(m)
185 
186  def set_file_dataset(self, fname, dataset):
187  """Associate a dataset with a filename.
188  This can be used to identify how the file was produced (in many
189  cases IMP can determine this automatically from a file header or
190  other metadata, but not always). For example, a manually-produced
191  PDB file (not from the PDB database or Modeller) can be
192  identified this way.
193  @param fname filename
194  @dataset the ihm.dataset.Dataset object to associate.
195  """
196  self._file_dataset[os.path.abspath(fname)] = dataset
197 
198  def get_file_dataset(self, fname):
199  """Get the dataset associated with a filename, or None.
200  @param fname filename
201  @return an ihm.dataset.Dataset, or None.
202  """
203  return self._file_dataset.get(os.path.abspath(fname), None)
204 
205  def add_protocol_output(self, p):
206  """Capture details of the modeling protocol.
207  @param p an instance of IMP.pmi.output.ProtocolOutput or a subclass.
208  """
209  state = p._add_state(self)
210  self._protocol_output.append((p, state))
211  p._each_metadata.append(self._metadata)
212  p._file_datasets.append(self._file_dataset)
213  state.model = self.model
214  state.prot = self.prot
215  protocol_output = property(lambda self:
216  [x[0] for x in self._protocol_output])
217 
218  def set_label(self, label):
219  self.label = label
220 
221  def create_component(self, name, color=0.0):
222  # Note that we set up the component as *both* a Chain and a Molecule.
223  # This is because old PMI1 code expects the top-level particle to be
224  # a Molecule, but we also need it to be a Chain to set the sequence.
225  # This looks a little odd but is a valid IMP atom::Hierarchy.
226  protein_h = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), 'X')
228  protein_h.set_name(name)
229  self.hier_dict[name] = protein_h
230  self.hier_representation[name] = {}
231  self.hier_db.add_name(name)
232  self.prot.add_child(protein_h)
233  self.color_dict[name] = color
234  self.elements[name] = []
235  for p, state in self._protocol_output:
236  p.create_component(state, name, True)
237 
238  def create_transformed_component(self, name, original, transform):
239  """Create a transformed view of a component.
240  This does not copy the coordinates or representation of the
241  component, and the transformed component is not visible to IMP,
242  but when the model is written to a ProtocolOutput, the transformed
243  component is added. This can be used to easily add symmetry-related
244  'copies' of a component without copying the underlying
245  IMP objects."""
246  for p, state in self._protocol_output:
247  p.create_transformed_component(state, name, original, transform)
248 
250  """Create a component that isn't used in the modeling.
251  No coordinates or other structural data for this component will
252  be read or written, but a primary sequence can be assigned. This
253  is useful if the input experimental data is of a system larger
254  than that modeled. Any references to these non-modeled components
255  can then be correctly resolved later."""
256  self._non_modeled_components[name] = None
257  self.elements[name] = []
258  for p, state in self._protocol_output:
259  p.create_component(state, name, False)
260 
261  def get_component_names(self):
262  return list(self.hier_dict.keys())
263 
264  def add_component_sequence(self, name, filename, id=None, offs=None,
265  format="FASTA"):
266  '''
267  Add the primary sequence for a single component.
268 
269  @param name Human-readable name of the component
270  @param filename Name of the FASTA file
271  @param id Identifier of the sequence in the FASTA file header
272  (if not provided, use `name` instead)
273  '''
274  record_dict = IMP.pmi.topology.Sequences(filename)
275  if id is None:
276  id = name
277  if id not in record_dict:
278  raise KeyError("id %s not found in fasta file" % id)
279  length = len(record_dict[id])
280  self.sequence_dict[name] = str(record_dict[id])
281  # No Hierarchy for this component if it is non-modeled
282  if name not in self._non_modeled_components:
283  protein_h = self.hier_dict[name]
284  protein_h.set_sequence(self.sequence_dict[name])
285  if offs is not None:
286  offs_str="-"*offs
287  self.sequence_dict[name]=offs_str+self.sequence_dict[name]
288 
289  self.elements[name].append((length, length, " ", "end"))
290  for p, state in self._protocol_output:
291  p.add_component_sequence(state, name, self.sequence_dict[name])
292 
293  def autobuild_model(self, name, pdbname, chain,
294  resolutions=None, resrange=None,
295  missingbeadsize=20,
296  color=None, pdbresrange=None, offset=0,
297  show=False, isnucleicacid=False,
298  attachbeads=False):
299 
300  self.representation_is_modified = True
301  outhiers = []
302 
303  if color is None:
304  color = self.color_dict[name]
305  else:
306  self.color_dict[name] = color
307 
308  if resolutions is None:
309  resolutions = [1]
310  print("autobuild_model: constructing %s from pdb %s and chain %s" % (name, pdbname, str(chain)))
311 
312  # get the initial and end residues of the pdb
313  t = IMP.atom.read_pdb(pdbname, self.model,
315 
316  # find start and end indexes
317 
318  start = IMP.atom.Residue(
319  t.get_children()[0].get_children()[0]).get_index()
320  end = IMP.atom.Residue(
321  t.get_children()[0].get_children()[-1]).get_index()
322 
323  # check if resrange was defined, otherwise
324  # use the sequence, or the pdb resrange
325 
326  if resrange is None:
327  if name in self.sequence_dict:
328  resrange = (1, len(self.sequence_dict[name]))
329  else:
330  resrange = (start + offset, end + offset)
331  else:
332  if resrange[1] in (-1, 'END'):
333  resrange = (resrange[0],end)
334  start = resrange[0] - offset
335  end = resrange[1] - offset
336 
338  t,
339  resrange[0],
340  resrange[1])
341 
343  t,
344  resrange[0],
345  terminus="N")
347  t,
348  resrange[1],
349  terminus="C")
350  # Done with the PDB
351  IMP.atom.destroy(t)
352 
353  # construct pdb fragments and intervening beads
354  for n, g in enumerate(gaps):
355  first = g[0]
356  last = g[1]
357  if g[2] == "cont":
358  print("autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
359  outhiers += self.add_component_pdb(name, pdbname,
360  chain, resolutions=resolutions,
361  color=color, cacenters=True,
362  resrange=(first, last),
363  offset=offset, isnucleicacid=isnucleicacid)
364  elif g[2] == "gap" and n > 0:
365  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
366  parts = self.hier_db.get_particles_at_closest_resolution(name,
367  first + offset - 1,
368  1)
369  xyz = IMP.core.XYZ(parts[0]).get_coordinates()
370  outhiers += self.add_component_necklace(name,
371  first+offset, last+offset, missingbeadsize, incoord=xyz)
372 
373  elif g[2] == "gap" and n == 0:
374  # add pre-beads
375  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
376  outhiers += self.add_component_necklace(name,
377  first+offset, last+offset, missingbeadsize, incoord=xyznter)
378 
379  return outhiers
380 
381  def add_component_pdb(self, name, pdbname, chain, resolutions, color=None,
382  resrange=None, offset=0, cacenters=True, show=False,
383  isnucleicacid=False, readnonwateratoms=False,
384  read_ca_cb_only=False):
385  '''
386  Add a component that has an associated 3D structure in a PDB file.
387 
388  Reads the PDB, and constructs the fragments corresponding to contiguous
389  sequence stretches.
390 
391  @return a list of hierarchies.
392 
393  @param name (string) the name of the component
394  @param pdbname (string) the name of the PDB file
395  @param chain (string or integer) can be either a string (eg, "A")
396  or an integer (eg, 0 or 1) in case you want
397  to get the corresponding chain number in the PDB.
398  @param resolutions (integers) a list of integers that corresponds
399  to the resolutions that have to be generated
400  @param color (float from 0 to 1) the color applied to the
401  hierarchies generated
402  @param resrange (tuple of integers): the residue range to extract
403  from the PDB. It is a tuple (beg,end). If not specified,
404  all residues belonging to the specified chain are read.
405  @param offset (integer) specifies the residue index offset to be
406  applied when reading the PDB (the FASTA sequence is
407  assumed to start from residue 1, so use this parameter
408  if the PDB file does not start at residue 1)
409  @param cacenters (boolean) if True generates resolution=1 beads
410  centered on C-alpha atoms.
411  @param show (boolean) print out the molecular hierarchy at the end.
412  @param isnucleicacid (boolean) use True if you're reading a PDB
413  with nucleic acids.
414  @param readnonwateratoms (boolean) if True fixes some pathological PDB.
415  @param read_ca_cb_only (boolean) if True, only reads CA/CB
416  '''
417 
418  self.representation_is_modified = True
419  if color is None:
420  # if the color is not passed, then get the stored color
421  color = self.color_dict[name]
422  protein_h = self.hier_dict[name]
423  outhiers = []
424 
425  # determine selector
427  if read_ca_cb_only:
428  cacbsel = IMP.atom.OrPDBSelector(
431  sel = IMP.atom.AndPDBSelector(cacbsel, sel)
432  if type(chain) == str:
435  sel)
436  t = IMP.atom.read_pdb(pdbname, self.model, sel)
437 
438  # get the first and last residue
439  start = IMP.atom.Residue(
440  t.get_children()[0].get_children()[0]).get_index()
441  end = IMP.atom.Residue(
442  t.get_children()[0].get_children()[-1]).get_index()
443  c = IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
444  else:
445  t = IMP.atom.read_pdb(pdbname, self.model, sel)
446  c = IMP.atom.Chain(
447  IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
448 
449  # get the first and last residue
450  start = IMP.atom.Residue(c.get_children()[0]).get_index()
451  end = IMP.atom.Residue(c.get_children()[-1]).get_index()
452  chain = c.get_id()
453 
454  if not resrange is None:
455  if resrange[0] > start:
456  start = resrange[0]
457  if resrange[1] < end:
458  end = resrange[1]
459 
460  if not isnucleicacid:
461  # do what you have to do for proteins
462  sel = IMP.atom.Selection(
463  c,
464  residue_indexes=list(range(
465  start,
466  end + 1)),
467  atom_type=IMP.atom.AT_CA)
468 
469  else:
470  # do what you have to do for nucleic-acids
471  # to work, nucleic acids should not be indicated as HETATM in the pdb
472  sel = IMP.atom.Selection(
473  c,
474  residue_indexes=list(range(
475  start,
476  end + 1)),
477  atom_type=IMP.atom.AT_P)
478 
479  ps = sel.get_selected_particles()
480  if len(ps) == 0:
481  raise ValueError("%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
482  % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
483  c0 = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), "X")
484 
485  for p in ps:
486  par = IMP.atom.Atom(p).get_parent()
487  ri = IMP.atom.Residue(par).get_index()
488  # Move residue from original PDB hierarchy to new chain c0
489  IMP.atom.Residue(par).set_index(ri + offset)
490  if par.get_parent() != c0:
491  par.get_parent().remove_child(par)
492  c0.add_child(par)
493  start = start + offset
494  end = end + offset
495 
496  self.elements[name].append(
497  (start, end, pdbname.split("/")[-1] + ":" + chain, "pdb"))
498 
499  hiers = self.coarse_hierarchy(name, start, end,
500  resolutions, isnucleicacid, c0, protein_h, "pdb", color)
501  self._copy_pdb_provenance(t, hiers[0], offset)
502  outhiers += hiers
503  for p, state in self._protocol_output:
504  p.add_pdb_element(state, name, start, end, offset, pdbname, chain,
505  hiers[0])
506 
507  if show:
509 
510  # We cannot simply destroy(c0) since it might not be a well-behaved
511  # hierarchy; in some cases it could contain a given residue more than
512  # once (this is surely a bug but we need to keep this behavior for
513  # backwards compatibility).
514  residues = {}
515  for p in ps:
516  par = IMP.atom.Atom(p).get_parent()
517  residues[par] = None
518  for r in residues.keys():
519  IMP.atom.destroy(r)
520  self.model.remove_particle(c0)
521 
522  IMP.atom.destroy(t)
523 
524  return outhiers
525 
526  def _copy_pdb_provenance(self, pdb, h, offset):
527  """Copy the provenance information from the PDB to our hierarchy"""
528  c = IMP.atom.Chain(IMP.atom.get_by_type(pdb, IMP.atom.CHAIN_TYPE)[0])
529  prov = IMP.core.Provenanced(c).get_provenance()
530  newprov = IMP.core.create_clone(prov)
533  IMP.core.StructureProvenance(newprov).set_residue_offset(offset)
534 
535  def add_component_ideal_helix(
536  self,
537  name,
538  resolutions,
539  resrange,
540  color=None,
541  show=False):
542 
543  self.representation_is_modified = True
544  from math import pi, cos, sin
545 
546  protein_h = self.hier_dict[name]
547  outhiers = []
548  if color is None:
549  color = self.color_dict[name]
550 
551  start = resrange[0]
552  end = resrange[1]
553  self.elements[name].append((start, end, " ", "helix"))
554  c0 = IMP.atom.Chain.setup_particle(IMP.Particle(self.model), "X")
555  for n, res in enumerate(range(start, end + 1)):
556  if name in self.sequence_dict:
557  try:
558  rtstr = self.onetothree[
559  self.sequence_dict[name][res-1]]
560  except:
561  rtstr = "UNK"
562  rt = IMP.atom.ResidueType(rtstr)
563  else:
564  rt = IMP.atom.ResidueType("ALA")
565 
566  # get the residue volume
567  try:
569  # mass=IMP.atom.get_mass_from_residue_type(rt)
570  except IMP.ValueException:
572  IMP.atom.ResidueType("ALA"))
573  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
575 
576  r = IMP.atom.Residue.setup_particle(IMP.Particle(self.model), rt, res)
577  p = IMP.Particle(self.model)
579  x = 2.3 * cos(n * 2 * pi / 3.6)
580  y = 2.3 * sin(n * 2 * pi / 3.6)
581  z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
582  d.set_coordinates(IMP.algebra.Vector3D(x, y, z))
583  d.set_radius(radius)
584  # print d
585  a = IMP.atom.Atom.setup_particle(p, IMP.atom.AT_CA)
586  r.add_child(a)
587  c0.add_child(r)
588 
589  outhiers += self.coarse_hierarchy(name, start, end,
590  resolutions, False, c0, protein_h, "helix", color)
591 
592  if show:
594 
595  IMP.atom.destroy(c0)
596  return outhiers
597 
598  def add_component_beads(self, name, ds, colors=None, incoord=None):
599  """ add beads to the representation
600  @param name the component name
601  @param ds a list of tuples corresponding to the residue ranges
602  of the beads
603  @param colors a list of colors associated to the beads
604  @param incoord the coordinate tuple corresponding to the position
605  of the beads
606  """
607 
608  from math import pi
609  self.representation_is_modified = True
610 
611  protein_h = self.hier_dict[name]
612  outhiers = []
613  if colors is None:
614  colors = [self.color_dict[name]]
615 
616 
617  for n, dss in enumerate(ds):
618  ds_frag = (dss[0], dss[1])
619  self.elements[name].append((dss[0], dss[1], " ", "bead"))
620  prt = IMP.Particle(self.model)
621  if ds_frag[0] == ds_frag[1]:
622  # if the bead represent a single residue
623  if name in self.sequence_dict:
624  try:
625  rtstr = self.onetothree[
626  self.sequence_dict[name][ds_frag[0]-1]]
627  except:
628  rtstr = "UNK"
629  rt = IMP.atom.ResidueType(rtstr)
630  else:
631  rt = IMP.atom.ResidueType("ALA")
632  h = IMP.atom.Residue.setup_particle(prt, rt, ds_frag[0])
633  h.set_name(name + '_%i_bead' % (ds_frag[0]))
634  prt.set_name(name + '_%i_bead' % (ds_frag[0]))
635  resolution = 1
636  else:
638  h.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
639  prt.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
640  h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
641  resolution = len(h.get_residue_indexes())
642  if "Beads" not in self.hier_representation[name]:
644  root.set_name("Beads")
645  self.hier_representation[name]["Beads"] = root
646  protein_h.add_child(root)
647  self.hier_representation[name]["Beads"].add_child(h)
648 
649  for kk in range(ds_frag[0], ds_frag[1] + 1):
650  self.hier_db.add_particles(name, kk, resolution, [h])
651 
652  try:
653  clr = IMP.display.get_rgb_color(colors[n])
654  except:
655  clr = IMP.display.get_rgb_color(colors[0])
656 
658 
659  # decorate particles according to their resolution
660  IMP.pmi.Resolution.setup_particle(prt, resolution)
661 
663  ptem = IMP.core.XYZR(prt)
665  if resolution == 1:
666  try:
668  except IMP.ValueException:
670  IMP.atom.ResidueType("ALA"))
673  ptem.set_radius(radius)
674  else:
675  volume = IMP.atom.get_volume_from_mass(mass)
676  radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
678  ptem.set_radius(radius)
679  try:
680  if not tuple(incoord) is None:
681  ptem.set_coordinates(incoord)
682  except TypeError:
683  pass
686  self.floppy_bodies.append(prt)
687  IMP.core.XYZ(prt).set_coordinates_are_optimized(True)
688  outhiers += [h]
689 
690  for p, state in self._protocol_output:
691  p.add_bead_element(state, name, ds[0][0], ds[-1][1], len(ds),
692  outhiers[0])
693 
694  return outhiers
695 
696  def add_component_necklace(self, name, begin, end, length, color=None,
697  incoord=None):
698  '''
699  Generates a string of beads with given length.
700  '''
701  self.representation_is_modified = True
702  outhiers = []
703  if color is None:
704  colors=None
705  else:
706  colors=[color]
707  for chunk in IMP.pmi.tools.list_chunks_iterator(range(begin, end + 1), length):
708  outhiers += self.add_component_beads(name,
709  [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
710  return outhiers
711 
713  self, name, hierarchies=None, selection_tuples=None,
714  particles=None,
715  resolution=0.0, num_components=10,
716  inputfile=None, outputfile=None,
717  outputmap=None,
718  kernel_type=None,
719  covariance_type='full', voxel_size=1.0,
720  out_hier_name='',
721  sampled_points=1000000, num_iter=100,
722  simulation_res=1.0,
723  multiply_by_total_mass=True,
724  transform=None,
725  intermediate_map_fn=None,
726  density_ps_to_copy=None,
727  use_precomputed_gaussians=False):
728  '''
729  Sets up a Gaussian Mixture Model for this component.
730  Can specify input GMM file or it will be computed.
731  @param name component name
732  @param hierarchies set up GMM for some hierarchies
733  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
734  @param particles set up GMM for particles directly
735  @param resolution usual PMI resolution for selecting particles from the hierarchies
736  @param inputfile read the GMM from this file
737  @param outputfile fit and write the GMM to this file (text)
738  @param outputmap after fitting, create GMM density file (mrc)
739  @param kernel_type for creating the intermediate density (points are sampled to make GMM). Options are IMP.em.GAUSSIAN, IMP.em.SPHERE, and IMP.em.BINARIZED_SPHERE
740  @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
741  @param voxel_size for creating the intermediate density map and output map.
742  lower number increases accuracy but also rasterizing time grows
743  @param out_hier_name name of the output density hierarchy
744  @param sampled_points number of points to sample. more will increase accuracy and fitting time
745  @param num_iter num GMM iterations. more will increase accuracy and fitting time
746  @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
747  @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
748  @param intermediate_map_fn for debugging, this will write the intermediate (simulated) map
749  @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
750  @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
751  '''
752  import numpy as np
753  import sys
754  import IMP.em
755  import IMP.isd.gmm_tools
756 
757  # prepare output
758  self.representation_is_modified = True
759  out_hier = []
760  protein_h = self.hier_dict[name]
761  if "Densities" not in self.hier_representation[name]:
763  root.set_name("Densities")
764  self.hier_representation[name]["Densities"] = root
765  protein_h.add_child(root)
766 
767  # gather passed particles
768  fragment_particles = []
769  if not particles is None:
770  fragment_particles += particles
771  if not hierarchies is None:
772  fragment_particles += IMP.pmi.tools.select(
773  self, resolution=resolution,
774  hierarchies=hierarchies)
775  if not selection_tuples is None:
776  for st in selection_tuples:
777  fragment_particles += IMP.pmi.tools.select_by_tuple(
778  self, tupleselection=st,
779  resolution=resolution,
780  name_is_ambiguous=False)
781 
782  # compute or read gaussians
783  density_particles = []
784  if inputfile:
786  inputfile, density_particles,
787  self.model, transform)
788  elif density_ps_to_copy:
789  for ip in density_ps_to_copy:
790  p = IMP.Particle(self.model)
791  shape = IMP.core.Gaussian(ip).get_gaussian()
792  mass = IMP.atom.Mass(ip).get_mass()
795  density_particles.append(p)
796  elif use_precomputed_gaussians:
797  if len(fragment_particles) == 0:
798  print("add_component_density: no particle was selected")
799  return out_hier
800  for p in fragment_particles:
801  if not (IMP.atom.Fragment.get_is_setup(self.model,p.get_particle_index()) and
802  IMP.core.XYZ.get_is_setup(self.model,p.get_particle_index())):
803  raise Exception("The particles you selected must be Fragments and XYZs")
804  nres=len(IMP.atom.Fragment(self.model,p.get_particle_index()).get_residue_indexes())
805  pos=IMP.core.XYZ(self.model,p.get_particle_index()).get_coordinates()
806  density_particles=[]
807  try:
808  IMP.isd.get_data_path("beads/bead_%i.txt"%nres)
809  except:
810  raise Exception("We haven't created a bead file for",nres,"residues")
811  transform = IMP.algebra.Transformation3D(pos)
813  IMP.isd.get_data_path("beads/bead_%i.txt"%nres), density_particles,
814  self.model, transform)
815  else:
816  #compute the gaussians here
817  if len(fragment_particles) == 0:
818  print("add_component_density: no particle was selected")
819  return out_hier
820 
821  density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
822  self.model,
823  fragment_particles,
824  num_components,
825  sampled_points,
826  simulation_res,
827  voxel_size,
828  num_iter,
829  covariance_type,
830  multiply_by_total_mass,
831  outputmap,
832  outputfile)
833 
834  # prepare output hierarchy
836  s0.set_name(out_hier_name)
837  self.hier_representation[name]["Densities"].add_child(s0)
838  out_hier.append(s0)
839  for nps, p in enumerate(density_particles):
840  s0.add_child(p)
841  p.set_name(s0.get_name() + '_gaussian_%i' % nps)
842  return out_hier
843 
844  def get_component_density(self, name):
845  return self.hier_representation[name]["Densities"]
846 
847  def add_all_atom_densities(self, name, hierarchies=None,
848  selection_tuples=None,
849  particles=None,
850  resolution=0,
851  output_map=None,
852  voxel_size=1.0):
853  '''Decorates all specified particles as Gaussians directly.
854  @param name component name
855  @param hierarchies set up GMM for some hierarchies
856  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
857  @param particles set up GMM for particles directly
858  @param resolution usual PMI resolution for selecting particles from the hierarchies
859  '''
860 
861  import IMP.em
862  import numpy as np
863  import sys
864  from math import sqrt
865  self.representation_is_modified = True
866 
867  if particles is None:
868  fragment_particles = []
869  else:
870  fragment_particles = particles
871 
872  if not hierarchies is None:
873  fragment_particles += IMP.pmi.tools.select(
874  self, resolution=resolution,
875  hierarchies=hierarchies)
876 
877  if not selection_tuples is None:
878  for st in selection_tuples:
879  fragment_particles += IMP.pmi.tools.select_by_tuple(
880  self, tupleselection=st,
881  resolution=resolution,
882  name_is_ambiguous=False)
883 
884  if len(fragment_particles) == 0:
885  print("add all atom densities: no particle was selected")
886  return
887 
888  # create a spherical gaussian for each particle based on atom type
889  print('setting up all atom gaussians num_particles',len(fragment_particles))
890  for n,p in enumerate(fragment_particles):
891  if IMP.core.Gaussian.get_is_setup(p): continue
892  center=IMP.core.XYZ(p).get_coordinates()
893  rad=IMP.core.XYZR(p).get_radius()
894  mass=IMP.atom.Mass(p).get_mass()
898  print('setting up particle',p.get_name(), " as individual gaussian particle")
899 
900  if not output_map is None:
901  print('writing map to', output_map)
903  fragment_particles,
904  output_map,
905  voxel_size)
906 
907  def add_component_hierarchy_clone(self, name, hierarchy):
908  '''
909  Make a copy of a hierarchy and append it to a component.
910  '''
911  outhier = []
912  self.representation_is_modified = True
913  protein_h = self.hier_dict[name]
914  hierclone = IMP.atom.create_clone(hierarchy)
915  hierclone.set_name(hierclone.get_name() + "_clone")
916  protein_h.add_child(hierclone)
917  outhier.append(hierclone)
918 
919  psmain = IMP.atom.get_leaves(hierarchy)
920  psclone = IMP.atom.get_leaves(hierclone)
921 
922  # copying attributes
923  for n, pmain in enumerate(psmain):
924  pclone = psclone[n]
926  resolution = IMP.pmi.Resolution(pmain).get_resolution()
927  IMP.pmi.Resolution.setup_particle(pclone, resolution)
928  for kk in IMP.pmi.tools.get_residue_indexes(pclone):
929  self.hier_db.add_particles(
930  name,
931  kk,
933  [pclone])
934 
936  uncertainty = IMP.pmi.Uncertainty(pmain).get_uncertainty()
937  IMP.pmi.Uncertainty.setup_particle(pclone, uncertainty)
938 
940  symmetric = IMP.pmi.Symmetric(pmain).get_symmetric()
941  IMP.pmi.Symmetric.setup_particle(pclone, symmetric)
942 
943  return outhier
944 
945  def dump_particle_descriptors(self):
946  import numpy
947  import pickle
948  import IMP.isd
949  import IMP.isd.gmm_tools
950 
951  particles_attributes={}
952  floppy_body_attributes={}
953  gaussians=[]
954  for h in IMP.atom.get_leaves(self.prot):
955  leaf=h
956  name=h.get_name()
957  hroot=self.prot
958  hparent=h.get_parent()
959  while hparent != hroot:
960  hparent=h.get_parent()
961  name+="|"+hparent.get_name()
962  h=hparent
963  particles_attributes[name]={"COORDINATES":numpy.array(IMP.core.XYZR(leaf.get_particle()).get_coordinates()),
964  "RADIUS":IMP.core.XYZR(leaf.get_particle()).get_radius(),
965  "MASS":IMP.atom.Mass(leaf.get_particle()).get_mass()}
966  if IMP.core.Gaussian.get_is_setup(leaf.get_particle()):
967  gaussians.append(IMP.core.Gaussian(leaf.get_particle()))
968 
969  rigid_body_attributes={}
970  for rb in self.rigid_bodies:
971  name=rb.get_name()
972  rf=rb.get_reference_frame()
973  t=rf.get_transformation_to()
974  trans=t.get_translation()
975  rot=t.get_rotation()
976  rigid_body_attributes[name]={"TRANSLATION":numpy.array(trans),
977  "ROTATION":numpy.array(rot.get_quaternion()),
978  "COORDINATES_NONRIGID_MEMBER":{},
979  "COORDINATES_RIGID_MEMBER":{}}
980  for mi in rb.get_member_indexes():
981  rm=self.model.get_particle(mi)
983  name_part=rm.get_name()
984  xyz=[self.model.get_attribute(fk, rm) for fk in [IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]]
985  rigid_body_attributes[name]["COORDINATES_NONRIGID_MEMBER"][name_part]=numpy.array(xyz)
986  else:
987  name_part=rm.get_name()
988  xyz=IMP.core.XYZ(rm).get_coordinates()
989  rigid_body_attributes[name]["COORDINATES_RIGID_MEMBER"][name_part]=numpy.array(xyz)
990 
991 
992  IMP.isd.gmm_tools.write_gmm_to_text(gaussians,"model_gmm.txt")
993  pickle.dump(particles_attributes,
994  open("particles_attributes.pkl", "wb"))
995  pickle.dump(rigid_body_attributes,
996  open("rigid_body_attributes.pkl", "wb"))
997 
998 
999 
1000  def load_particle_descriptors(self):
1001  import numpy
1002  import pickle
1003  import IMP.isd
1004  import IMP.isd.gmm_tools
1005 
1006  particles_attributes = pickle.load(open("particles_attributes.pkl",
1007  "rb"))
1008  rigid_body_attributes = pickle.load(open("rigid_body_attributes.pkl",
1009  "rb"))
1010 
1011  particles=[]
1012  hierarchies=[]
1013  gaussians=[]
1014  for h in IMP.atom.get_leaves(self.prot):
1015  leaf=h
1016  name=h.get_name()
1017  hroot=self.prot
1018  hparent=h.get_parent()
1019  while hparent != hroot:
1020  hparent=h.get_parent()
1021  name+="|"+hparent.get_name()
1022  h=hparent
1023 
1024  xyzr=IMP.core.XYZR(leaf.get_particle())
1025  xyzr.set_coordinates(particles_attributes[name]["COORDINATES"])
1026  #xyzr.set_radius(particles_attributes[name]["RADIUS"])
1027  #IMP.atom.Mass(leaf.get_particle()).set_mass(particles_attributes[name]["MASS"])
1028  if IMP.core.Gaussian.get_is_setup(leaf.get_particle()):
1029  gaussians.append(IMP.core.Gaussian(leaf.get_particle()))
1030 
1031  for rb in self.rigid_bodies:
1032  name=rb.get_name()
1033  trans=rigid_body_attributes[name]["TRANSLATION"]
1034  rot=rigid_body_attributes[name]["ROTATION"]
1037  rb.set_reference_frame(rf)
1038  coor_nrm_ref=rigid_body_attributes[name]["COORDINATES_NONRIGID_MEMBER"]
1039  coor_rm_ref_dict=rigid_body_attributes[name]["COORDINATES_RIGID_MEMBER"]
1040  coor_rm_model=[]
1041  coor_rm_ref=[]
1042  for mi in rb.get_member_indexes():
1043  rm=self.model.get_particle(mi)
1045  name_part=rm.get_name()
1046  xyz=coor_nrm_ref[name_part]
1047  for n,fk in enumerate([IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]):
1048  self.model.set_attribute(fk, rm,xyz[n])
1049  else:
1050  name_part=rm.get_name()
1051  coor_rm_ref.append(IMP.algebra.Vector3D(coor_rm_ref_dict[name_part]))
1052  coor_rm_model.append(IMP.core.XYZ(rm).get_coordinates())
1053  if len(coor_rm_model)==0: continue
1054  t=IMP.algebra.get_transformation_aligning_first_to_second(coor_rm_model,coor_rm_ref)
1055  IMP.core.transform(rb,t)
1056 
1057  IMP.isd.gmm_tools.decorate_gmm_from_text("model_gmm.txt",gaussians,self.model)
1058 
1059  def _compare_rmf_repr_names(self, rmfname, reprname, component_name):
1060  """Print a warning if particle names in RMF and model don't match"""
1061  def match_any_suffix():
1062  # Handle common mismatches like 743 != Nup85_743_pdb
1063  suffixes = ["pdb", "bead_floppy_body_rigid_body_member_floppy_body",
1064  "bead_floppy_body_rigid_body_member",
1065  "bead_floppy_body"]
1066  for s in suffixes:
1067  if "%s_%s_%s" % (component_name, rmfname, s) == reprname:
1068  return True
1069  if rmfname != reprname and not match_any_suffix():
1070  print("set_coordinates_from_rmf: WARNING rmf particle and "
1071  "representation particle names don't match %s %s"
1072  % (rmfname, reprname))
1073 
1074  def set_coordinates_from_rmf(self, component_name, rmf_file_name,
1075  rmf_frame_number,
1076  rmf_component_name=None,
1077  check_number_particles=True,
1078  representation_name_to_rmf_name_map=None,
1079  state_number=0,
1080  skip_gaussian_in_rmf=False,
1081  skip_gaussian_in_representation=False,
1082  save_file=False,
1083  force_rigid_update=False):
1084  '''Read and replace coordinates from an RMF file.
1085  Replace the coordinates of particles with the same name.
1086  It assumes that the RMF and the representation have the particles
1087  in the same order.
1088  @param component_name Component name
1089  @param rmf_component_name Name of the component in the RMF file
1090  (if not specified, use `component_name`)
1091  @param representation_name_to_rmf_name_map a dictionary that map
1092  the original rmf particle name to the recipient particle component name
1093  @param save_file: save a file with the names of particles of the component
1094  @param force_rigid_update: update the coordinates of rigid bodies
1095  (normally this should be called before rigid bodies are set up)
1096  '''
1097  import IMP.pmi.analysis
1098 
1099  tempm = IMP.Model()
1100  prots = IMP.pmi.analysis.get_hiers_from_rmf(
1101  tempm,
1102  rmf_frame_number,
1103  rmf_file_name)
1104 
1105  if not prots:
1106  raise ValueError("cannot read hierarchy from rmf")
1107 
1108  prot=prots[0]
1109  # Make sure coordinates of rigid body members in the RMF are correct
1110  if force_rigid_update:
1111  tempm.update()
1112 
1113  # if len(self.rigid_bodies)!=0:
1114  # print "set_coordinates_from_rmf: cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies"
1115  # exit()
1116 
1117  allpsrmf = IMP.atom.get_leaves(prot)
1118  psrmf = []
1119  for p in allpsrmf:
1120  (protname, is_a_bead) = IMP.pmi.tools.get_prot_name_from_particle(
1121  p, self.hier_dict.keys())
1122  if (protname is None) and (rmf_component_name is not None):
1123  (protname, is_a_bead) = IMP.pmi.tools.get_prot_name_from_particle(
1124  p, rmf_component_name)
1125  if (skip_gaussian_in_rmf):
1127  continue
1128  if (rmf_component_name is not None) and (protname == rmf_component_name):
1129  psrmf.append(p)
1130  elif (rmf_component_name is None) and (protname == component_name):
1131  psrmf.append(p)
1132 
1133  psrepr = IMP.atom.get_leaves(self.hier_dict[component_name])
1134  if (skip_gaussian_in_representation):
1135  allpsrepr = psrepr
1136  psrepr = []
1137  for p in allpsrepr:
1138  #(protname, is_a_bead) = IMP.pmi.tools.get_prot_name_from_particle(
1139  # p, self.hier_dict.keys())
1141  continue
1142  psrepr.append(p)
1143 
1144  import itertools
1145  reprnames=[p.get_name() for p in psrepr]
1146  rmfnames=[p.get_name() for p in psrmf]
1147 
1148  if save_file:
1149  fl=open(component_name+".txt","w")
1150  for i in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+","+str(i[1])+"\n")
1151 
1152 
1153  if check_number_particles and not representation_name_to_rmf_name_map:
1154  if len(psrmf) != len(psrepr):
1155  fl=open(component_name+".txt","w")
1156  for i in itertools.izip_longest(reprnames,rmfnames): fl.write(str(i[0])+","+str(i[1])+"\n")
1157  raise ValueError("%s cannot proceed the rmf and the representation don't have the same number of particles; "
1158  "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
1159 
1160 
1161  if not representation_name_to_rmf_name_map:
1162  for n, prmf in enumerate(psrmf):
1163 
1164  prmfname = prmf.get_name()
1165  preprname = psrepr[n].get_name()
1166  if force_rigid_update:
1167  if IMP.core.RigidBody.get_is_setup(psrepr[n]) \
1168  and not IMP.core.RigidBodyMember.get_is_setup(psrepr[n]):
1169  continue
1170  else:
1172  raise ValueError("component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
1173 
1174  self._compare_rmf_repr_names(prmfname, preprname,
1175  component_name)
1176  if IMP.core.XYZ.get_is_setup(prmf) and IMP.core.XYZ.get_is_setup(psrepr[n]):
1177  xyz = IMP.core.XYZ(prmf).get_coordinates()
1178  IMP.core.XYZ(psrepr[n]).set_coordinates(xyz)
1180  # Set rigid body so that coordinates are preserved
1181  # on future model updates
1182  rbm = IMP.core.RigidBodyMember(psrepr[n])
1183  rbm.set_internal_coordinates(xyz)
1185  rb = rbm.get_rigid_body()
1187  raise ValueError("Cannot handle nested "
1188  "rigid bodies yet")
1189  rb.set_reference_frame_lazy(tr)
1190  else:
1191  print("set_coordinates_from_rmf: WARNING particles are not XYZ decorated %s %s " % (str(IMP.core.XYZ.get_is_setup(prmf)), str(IMP.core.XYZ.get_is_setup(psrepr[n]))))
1192 
1194  gprmf=IMP.core.Gaussian(prmf)
1195  grepr=IMP.core.Gaussian(psrepr[n])
1196  g=gprmf.get_gaussian()
1197  grepr.set_gaussian(g)
1198 
1199  else:
1200  repr_name_particle_map={}
1201  rmf_name_particle_map={}
1202  for p in psrmf:
1203  rmf_name_particle_map[p.get_name()]=p
1204  #for p in psrepr:
1205  # repr_name_particle_map[p.get_name()]=p
1206 
1207  for prepr in psrepr:
1208  try:
1209  prmf=rmf_name_particle_map[representation_name_to_rmf_name_map[prepr.get_name()]]
1210  except KeyError:
1211  print("set_coordinates_from_rmf: WARNING missing particle name in representation_name_to_rmf_name_map, skipping...")
1212  continue
1213  xyz = IMP.core.XYZ(prmf).get_coordinates()
1214  IMP.core.XYZ(prepr).set_coordinates(xyz)
1215 
1216 
1217 
1218  def check_root(self, name, protein_h, resolution):
1219  '''
1220  If the root hierarchy does not exist, construct it.
1221  '''
1222  if "Res:" + str(int(resolution)) not in self.hier_representation[name]:
1224  root.set_name(name + "_Res:" + str(int(resolution)))
1225  self.hier_representation[name][
1226  "Res:" + str(int(resolution))] = root
1227  protein_h.add_child(root)
1228 
1229  def coarse_hierarchy(self, name, start, end, resolutions, isnucleicacid,
1230  input_hierarchy, protein_h, type, color):
1231  '''
1232  Generate all needed coarse grained layers.
1233 
1234  @param name name of the protein
1235  @param resolutions list of resolutions
1236  @param protein_h root hierarchy
1237  @param input_hierarchy hierarchy to coarse grain
1238  @param type a string, typically "pdb" or "helix"
1239  '''
1240  outhiers = []
1241 
1242  if (1 in resolutions) or (0 in resolutions):
1243  # in that case create residues and append atoms
1244 
1245  if 1 in resolutions:
1246  self.check_root(name, protein_h, 1)
1248  s1.set_name('%s_%i-%i_%s' % (name, start, end, type))
1249  # s1.set_residue_indexes(range(start,end+1))
1250  self.hier_representation[name]["Res:1"].add_child(s1)
1251  outhiers += [s1]
1252  if 0 in resolutions:
1253  self.check_root(name, protein_h, 0)
1255  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1256  # s0.set_residue_indexes(range(start,end+1))
1257  self.hier_representation[name]["Res:0"].add_child(s0)
1258  outhiers += [s0]
1259 
1260  if not isnucleicacid:
1261  sel = IMP.atom.Selection(
1262  input_hierarchy,
1263  atom_type=IMP.atom.AT_CA)
1264  else:
1265  sel = IMP.atom.Selection(
1266  input_hierarchy,
1267  atom_type=IMP.atom.AT_P)
1268 
1269  for p in sel.get_selected_particles():
1270  resobject = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
1271  if 0 in resolutions:
1272  # if you ask for atoms
1273  resclone0 = IMP.atom.create_clone(resobject)
1274  resindex = IMP.atom.Residue(resclone0).get_index()
1275  s0.add_child(resclone0)
1276  self.hier_db.add_particles(
1277  name,
1278  resindex,
1279  0,
1280  resclone0.get_children())
1281 
1282  chil = resclone0.get_children()
1283  for ch in chil:
1285  try:
1286  clr = IMP.display.get_rgb_color(color)
1287  except:
1288  clr = IMP.display.get_rgb_color(1.0)
1290 
1291  if 1 in resolutions:
1292  # else clone the residue
1293  resclone1 = IMP.atom.create_clone_one(resobject)
1294  resindex = IMP.atom.Residue(resclone1).get_index()
1295  s1.add_child(resclone1)
1296  self.hier_db.add_particles(name, resindex, 1, [resclone1])
1297 
1298  rt = IMP.atom.Residue(resclone1).get_residue_type()
1299  xyz = IMP.core.XYZ(p).get_coordinates()
1300  prt = resclone1.get_particle()
1301  prt.set_name('%s_%i_%s' % (name, resindex, type))
1302  IMP.core.XYZ.setup_particle(prt).set_coordinates(xyz)
1303 
1304  try:
1306  # mass=IMP.atom.get_mass_from_residue_type(rt)
1307  except IMP.ValueException:
1309  IMP.atom.ResidueType("ALA"))
1310  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
1312  IMP.core.XYZR.setup_particle(prt).set_radius(radius)
1317  try:
1318  clr = IMP.display.get_rgb_color(color)
1319  except:
1320  clr = IMP.display.get_rgb_color(1.0)
1322 
1323  for r in resolutions:
1324  if r != 0 and r != 1:
1325  self.check_root(name, protein_h, r)
1327  input_hierarchy,
1328  r)
1329 
1330  chil = s.get_children()
1332  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1333  # Move all children from s to s0
1334  for ch in chil:
1335  s.remove_child(ch)
1336  s0.add_child(ch)
1337  self.hier_representation[name][
1338  "Res:" + str(int(r))].add_child(s0)
1339  outhiers += [s0]
1340  IMP.atom.destroy(s)
1341 
1342  for prt in IMP.atom.get_leaves(s0):
1343  ri = IMP.atom.Fragment(prt).get_residue_indexes()
1344  first = ri[0]
1345  last = ri[-1]
1346  if first == last:
1347  prt.set_name('%s_%i_%s' % (name, first, type))
1348  else:
1349  prt.set_name('%s_%i_%i_%s' % (name, first, last, type))
1350  for kk in ri:
1351  self.hier_db.add_particles(name, kk, r, [prt])
1352 
1353  radius = IMP.core.XYZR(prt).get_radius()
1357 
1358  # setting up color for each particle in the
1359  # hierarchy, if colors missing in the colors list set it to
1360  # red
1361  try:
1362  clr = IMP.display.get_rgb_color(color)
1363  except:
1364  colors.append(1.0)
1365  clr = IMP.display.get_rgb_color(colors[pdb_part_count])
1367 
1368  return outhiers
1369 
1371  '''
1372  Get the hierarchies at the given resolution.
1373 
1374  The map between resolution and hierarchies is cached to accelerate
1375  the selection algorithm. The cache is invalidated when the
1376  representation was changed by any add_component_xxx.
1377  '''
1378 
1379  if self.representation_is_modified:
1380  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1381  resolution)
1382  self.hier_resolution[resolution] = rhbr
1383  self.representation_is_modified = False
1384  return rhbr
1385  else:
1386  if resolution in self.hier_resolution:
1387  return self.hier_resolution[resolution]
1388  else:
1389  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1390  resolution)
1391  self.hier_resolution[resolution] = rhbr
1392  return rhbr
1393 
1395  self, max_translation=300., max_rotation=2.0 * pi,
1396  avoidcollision=True, cutoff=10.0, niterations=100,
1397  bounding_box=None,
1398  excluded_rigid_bodies=None,
1399  ignore_initial_coordinates=False,
1400  hierarchies_excluded_from_collision=None):
1401  '''
1402  Shuffle configuration; used to restart the optimization.
1403 
1404  The configuration of the system is initialized by placing each
1405  rigid body and each bead randomly in a box with a side of
1406  `max_translation` angstroms, and far enough from each other to
1407  prevent any steric clashes. The rigid bodies are also randomly rotated.
1408 
1409  @param avoidcollision check if the particle/rigid body was
1410  placed close to another particle; uses the optional
1411  arguments cutoff and niterations
1412  @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1413  '''
1414 
1415  if excluded_rigid_bodies is None:
1416  excluded_rigid_bodies = []
1417  if hierarchies_excluded_from_collision is None:
1418  hierarchies_excluded_from_collision = []
1419 
1420  if len(self.rigid_bodies) == 0:
1421  print("shuffle_configuration: rigid bodies were not intialized")
1422 
1424  gcpf.set_distance(cutoff)
1425  ps = []
1426  hierarchies_excluded_from_collision_indexes = []
1427  for p in IMP.atom.get_leaves(self.prot):
1429  ps.append(p)
1431  # remove the densities particles out of the calculation
1432  hierarchies_excluded_from_collision_indexes += IMP.get_indexes([p])
1433  allparticleindexes = IMP.get_indexes(ps)
1434 
1435  if bounding_box is not None:
1436  ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1437  lb = IMP.algebra.Vector3D(x1, y1, z1)
1438  ub = IMP.algebra.Vector3D(x2, y2, z2)
1439  bb = IMP.algebra.BoundingBox3D(lb, ub)
1440 
1441  for h in hierarchies_excluded_from_collision:
1442  hierarchies_excluded_from_collision_indexes += IMP.get_indexes(IMP.atom.get_leaves(h))
1443 
1444 
1445  allparticleindexes = list(
1446  set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1447 
1448  print(hierarchies_excluded_from_collision)
1449  print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1450 
1451  print('shuffling', len(self.rigid_bodies) - len(excluded_rigid_bodies), 'rigid bodies')
1452  for rb in self.rigid_bodies:
1453  if rb not in excluded_rigid_bodies:
1454  if avoidcollision:
1455  rbindexes = rb.get_member_particle_indexes()
1456  rbindexes = list(
1457  set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1458  otherparticleindexes = list(
1459  set(allparticleindexes) - set(rbindexes))
1460 
1461  if len(otherparticleindexes) is None:
1462  continue
1463 
1464  niter = 0
1465  while niter < niterations:
1466  if (ignore_initial_coordinates):
1467  # Move the particle to the origin
1468  transformation = IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(), -IMP.core.XYZ(rb).get_coordinates())
1469  IMP.core.transform(rb, transformation)
1470  rbxyz = IMP.core.XYZ(rb).get_coordinates()
1471 
1472  if bounding_box is not None:
1473  # overrides the perturbation
1474  translation = IMP.algebra.get_random_vector_in(bb)
1476  transformation = IMP.algebra.Transformation3D(rotation, translation-rbxyz)
1477  else:
1479  rbxyz,
1480  max_translation,
1481  max_rotation)
1482 
1483  IMP.core.transform(rb, transformation)
1484 
1485  if avoidcollision:
1486  self.model.update()
1487  npairs = len(
1488  gcpf.get_close_pairs(
1489  self.model,
1490  otherparticleindexes,
1491  rbindexes))
1492  if npairs == 0:
1493  niter = niterations
1494  if (ignore_initial_coordinates):
1495  print (rb.get_name(), IMP.core.XYZ(rb).get_coordinates())
1496  else:
1497  niter += 1
1498  print("shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1499  print("shuffle_configuration: rigid body name: " + rb.get_name())
1500  if niter == niterations:
1501  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1502  else:
1503  break
1504 
1505  print('shuffling', len(self.floppy_bodies), 'floppy bodies')
1506  for fb in self.floppy_bodies:
1507  if (avoidcollision):
1510  if rm and nrm:
1511  fbindexes = IMP.get_indexes([fb])
1512  otherparticleindexes = list(
1513  set(allparticleindexes) - set(fbindexes))
1514  if len(otherparticleindexes) is None:
1515  continue
1516  else:
1517  continue
1518  else:
1521  if (rm or nrm):
1522  continue
1523 
1525  d=IMP.core.RigidBodyMember(fb).get_rigid_body()
1527  d=IMP.core.RigidBody(fb)
1528  elif IMP.core.XYZ.get_is_setup(fb):
1529  d=IMP.core.XYZ(fb)
1530 
1531  niter = 0
1532  while niter < niterations:
1533  if (ignore_initial_coordinates):
1534  # Move the particle to the origin
1535  transformation = IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(), -IMP.core.XYZ(fb).get_coordinates())
1536  IMP.core.transform(d, transformation)
1537  fbxyz = IMP.core.XYZ(fb).get_coordinates()
1538 
1539  if bounding_box is not None:
1540  # overrides the perturbation
1541  translation = IMP.algebra.get_random_vector_in(bb)
1542  transformation = IMP.algebra.Transformation3D(translation-fbxyz)
1543  else:
1545  fbxyz,
1546  max_translation,
1547  max_rotation)
1548 
1549  IMP.core.transform(d, transformation)
1550 
1551  if (avoidcollision):
1552  self.model.update()
1553  npairs = len(
1554  gcpf.get_close_pairs(
1555  self.model,
1556  otherparticleindexes,
1557  fbindexes))
1558  if npairs == 0:
1559  niter = niterations
1560  if (ignore_initial_coordinates):
1561  print (fb.get_name(), IMP.core.XYZ(fb).get_coordinates())
1562  else:
1563  niter += 1
1564  print("shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1565  print("shuffle_configuration: floppy body name: " + fb.get_name())
1566  if niter == niterations:
1567  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1568  else:
1569  break
1570 
1571  def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1572  # getting only coordinates from pdb
1573  ps = IMP.pmi.tools.select(self, resolution=1.0)
1574  # storing the reference coordinates and the particles
1575  self.reference_structures[label] = (
1576  [IMP.core.XYZ(p).get_coordinates() for p in ps],
1577  ps)
1578 
1579  def get_all_rmsds(self):
1580  rmsds = {}
1581 
1582  for label in self.reference_structures:
1583 
1584  current_coordinates = [IMP.core.XYZ(p).get_coordinates()
1585  for p in self.reference_structures[label][1]]
1586  reference_coordinates = self.reference_structures[label][0]
1587  if len(reference_coordinates) != len(current_coordinates):
1588  print("calculate_all_rmsds: reference and actual coordinates are not the same")
1589  continue
1591  current_coordinates,
1592  reference_coordinates)
1593  rmsd_global = IMP.algebra.get_rmsd(
1594  reference_coordinates,
1595  current_coordinates)
1596  # warning: temporary we are calculating the drms, and not the rmsd,
1597  # for the relative distance
1598  rmsd_relative = IMP.atom.get_drms(
1599  reference_coordinates,
1600  current_coordinates)
1601  rmsds[label + "_GlobalRMSD"] = rmsd_global
1602  rmsds[label + "_RelativeDRMS"] = rmsd_relative
1603  return rmsds
1604 
1605  def setup_component_geometry(self, name, color=None, resolution=1.0):
1606  if color is None:
1607  color = self.color_dict[name]
1608  # this function stores all particle pairs
1609  # ordered by residue number, to be used
1610  # to construct backbone traces
1611  self.hier_geometry_pairs[name] = []
1612  protein_h = self.hier_dict[name]
1613  pbr = IMP.pmi.tools.select(self, name=name, resolution=resolution)
1614  pbr = IMP.pmi.tools.sort_by_residues(pbr)
1615 
1616  for n in range(len(pbr) - 1):
1617  self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1618 
1620  self, name, resolution=10, scale=1.0):
1621  '''
1622  Generate restraints between contiguous fragments in the hierarchy.
1623  The linkers are generated at resolution 10 by default.
1624 
1625  '''
1626 
1627  unmodeledregions_cr = IMP.RestraintSet(self.model, "unmodeledregions")
1628  sortedsegments_cr = IMP.RestraintSet(self.model, "sortedsegments")
1629 
1630  protein_h = self.hier_dict[name]
1631  SortedSegments = []
1632  frs = self.hier_db.get_preroot_fragments_by_resolution(
1633  name,
1634  resolution)
1635 
1636  for fr in frs:
1637  try:
1638  start = fr.get_children()[0]
1639  except:
1640  start = fr
1641 
1642  try:
1643  end = fr.get_children()[-1]
1644  except:
1645  end = fr
1646 
1647  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
1648  endres = IMP.pmi.tools.get_residue_indexes(end)[-1]
1649  SortedSegments.append((start, end, startres))
1650  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1651 
1652  # connect the particles
1653  for x in range(len(SortedSegments) - 1):
1654  last = SortedSegments[x][1]
1655  first = SortedSegments[x + 1][0]
1656 
1657  nreslast = len(IMP.pmi.tools.get_residue_indexes(last))
1658  lastresn = IMP.pmi.tools.get_residue_indexes(last)[-1]
1659  nresfirst = len(IMP.pmi.tools.get_residue_indexes(first))
1660  firstresn = IMP.pmi.tools.get_residue_indexes(first)[0]
1661 
1662  residuegap = firstresn - lastresn - 1
1663  if self.disorderedlength and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1664  # calculate the distance between the sphere centers using Kohn
1665  # PNAS 2004
1666  optdist = sqrt(5 / 3) * 1.93 * \
1667  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1668  # optdist2=sqrt(5/3)*1.93*((nreslast)**0.6+(nresfirst)**0.6)/2
1669  if self.upperharmonic:
1670  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1671  else:
1672  hu = IMP.core.Harmonic(optdist, self.kappa)
1673  dps = IMP.core.DistancePairScore(hu)
1674  else: # default
1675  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1676  if self.upperharmonic: # default
1677  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1678  else:
1679  hu = IMP.core.Harmonic(optdist, self.kappa)
1681 
1682  pt0 = last.get_particle()
1683  pt1 = first.get_particle()
1684  r = IMP.core.PairRestraint(self.model, dps, (pt0.get_index(), pt1.get_index()))
1685 
1686  print("Adding sequence connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
1687  sortedsegments_cr.add_restraint(r)
1688  self.linker_restraints_dict[
1689  "LinkerRestraint-" + pt0.get_name() + "-" + pt1.get_name()] = r
1690  self.connected_intra_pairs.append((pt0, pt1))
1691  self.connected_intra_pairs.append((pt1, pt0))
1692 
1693  self.linker_restraints.add_restraint(sortedsegments_cr)
1694  self.linker_restraints.add_restraint(unmodeledregions_cr)
1695  IMP.pmi.tools.add_restraint_to_model(self.model, self.linker_restraints)
1696  self.sortedsegments_cr_dict[name] = sortedsegments_cr
1697  self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1698 
1699  def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1700  import IMP.pmi.samplers
1701  pts = IMP.pmi.tools.ParticleToSampleList()
1702  for n, fb in enumerate(self.floppy_bodies):
1703  pts.add_particle(fb, "Floppy_Bodies", 1.0, "Floppy_Body_" + str(n))
1704  if len(pts.get_particles_to_sample()) > 0:
1705  mc = IMP.pmi.samplers.MonteCarlo(self.model, [pts], temperature)
1706  print("optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1707  mc.optimize(nsteps)
1708  else:
1709  print("optimize_floppy_bodies: no particle to optimize")
1710 
1711  def create_rotational_symmetry(self, maincopy, copies, rotational_axis=IMP.algebra.Vector3D(0, 0, 1.0),
1712  nSymmetry=None, skip_gaussian_in_clones=False):
1713  '''
1714  The copies must not contain rigid bodies.
1715  The symmetry restraints are applied at each leaf.
1716  '''
1717 
1718  from math import pi
1719  self.representation_is_modified = True
1720  ncopies = len(copies) + 1
1721  main_hiers = IMP.atom.get_leaves(self.hier_dict[maincopy])
1722 
1723  for k in range(len(copies)):
1724  if (nSymmetry is None):
1725  rotation_angle = 2.0 * pi / float(ncopies) * float(k + 1)
1726  else:
1727  if ( k % 2 == 0 ):
1728  rotation_angle = 2.0 * pi / float(nSymmetry) * float((k + 2) / 2)
1729  else:
1730  rotation_angle = -2.0 * pi / float(nSymmetry) * float((k + 1) / 2)
1731 
1732  rotation3D = IMP.algebra.get_rotation_about_axis(rotational_axis, rotation_angle)
1733 
1734  sm = IMP.core.TransformationSymmetry(rotation3D)
1735  clone_hiers = IMP.atom.get_leaves(self.hier_dict[copies[k]])
1736 
1737  lc = IMP.container.ListSingletonContainer(self.model)
1738  for n, p in enumerate(main_hiers):
1739  if (skip_gaussian_in_clones):
1741  continue
1742  pc = clone_hiers[n]
1743  #print("setting " + p.get_name() + " as reference for " + pc.get_name())
1744  IMP.core.Reference.setup_particle(pc.get_particle(), p.get_particle())
1745  lc.add(pc.get_particle().get_index())
1746 
1747  c = IMP.container.SingletonsConstraint(sm, None, lc)
1748  self.model.add_score_state(c)
1749  print("Completed setting " + str(maincopy) + " as a reference for "
1750  + str(copies[k]) + " by rotating it by "
1751  + str(rotation_angle / 2.0 / pi * 360)
1752  + " degrees around the " + str(rotational_axis) + " axis.")
1753  self.model.update()
1754 
1755  def create_rigid_body_symmetry(self, particles_reference, particles_copy,label="None",
1756  initial_transformation=IMP.algebra.get_identity_transformation_3d()):
1757  from math import pi
1758  self.representation_is_modified = True
1759 
1760  mainparticles = particles_reference
1761 
1762  t=initial_transformation
1763  p=IMP.Particle(self.model)
1764  p.set_name("RigidBody_Symmetry")
1766 
1768 
1769 
1770  copyparticles = particles_copy
1771 
1772  mainpurged = []
1773  copypurged = []
1774  for n, p in enumerate(mainparticles):
1775  print(p.get_name())
1776  pc = copyparticles[n]
1777 
1778  mainpurged.append(p)
1781  else:
1782  IMP.pmi.Symmetric(p).set_symmetric(0)
1783 
1784  copypurged.append(pc)
1785  if not IMP.pmi.Symmetric.get_is_setup(pc):
1787  else:
1788  IMP.pmi.Symmetric(pc).set_symmetric(1)
1789 
1790  lc = IMP.container.ListSingletonContainer(self.model)
1791  for n, p in enumerate(mainpurged):
1792 
1793  pc = copypurged[n]
1794  print("setting " + p.get_name() + " as reference for " + pc.get_name())
1795 
1797  lc.add(pc.get_index())
1798 
1799  c = IMP.container.SingletonsConstraint(sm, None, lc)
1800  self.model.add_score_state(c)
1801 
1802  self.model.update()
1803  self.rigid_bodies.append(rb)
1804  self.rigid_body_symmetries.append(rb)
1805  rb.set_name(label+".rigid_body_symmetry."+str(len(self.rigid_body_symmetries)))
1806 
1807 
1808  def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1809  longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1810 
1811  from math import pi
1812  self.representation_is_modified = True
1813 
1814  outhiers = []
1815  protein_h = self.hier_dict[maincopy]
1816  mainparts = IMP.atom.get_leaves(protein_h)
1817 
1818  for ilc in range(-longitudinal_copies, longitudinal_copies + 1):
1819  for iac in range(axial_copies):
1820  copyname = maincopy + "_a" + str(ilc) + "_l" + str(iac)
1821  self.create_component(copyname, 0.0)
1822  for hier in protein_h.get_children():
1823  self.add_component_hierarchy_clone(copyname, hier)
1824  copyhier = self.hier_dict[copyname]
1825  outhiers.append(copyhier)
1826  copyparts = IMP.atom.get_leaves(copyhier)
1828  IMP.algebra.Vector3D(axis),
1829  2 * pi / axial_copies * (float(iac)))
1830  translation_vector = tuple(
1831  [translation_value * float(ilc) * x for x in axis])
1832  print(translation_vector)
1833  translation = IMP.algebra.Vector3D(translation_vector)
1835  IMP.algebra.Transformation3D(rotation3D, translation))
1836  lc = IMP.container.ListSingletonContainer(self.model)
1837  for n, p in enumerate(mainparts):
1838  pc = copyparts[n]
1841  if not IMP.pmi.Symmetric.get_is_setup(pc):
1844  lc.add(pc.get_index())
1845  c = IMP.container.SingletonsConstraint(sm, None, lc)
1846  self.model.add_score_state(c)
1847  self.model.update()
1848  return outhiers
1849 
1850  def link_components_to_rmf(self, rmfname, frameindex):
1851  '''
1852  Load coordinates in the current representation.
1853  This should be done only if the hierarchy self.prot is identical
1854  to the one as stored in the rmf i.e. all components were added.
1855  '''
1856  import IMP.rmf
1857  import RMF
1858 
1859  rh = RMF.open_rmf_file_read_only(rmfname)
1860  IMP.rmf.link_hierarchies(rh, [self.prot])
1861  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1862  del rh
1863 
1864  def create_components_from_rmf(self, rmfname, frameindex):
1865  '''
1866  still not working.
1867  create the representation (i.e. hierarchies) from the rmf file.
1868  it will be stored in self.prot, which will be overwritten.
1869  load the coordinates from the rmf file at frameindex.
1870  '''
1871  rh = RMF.open_rmf_file_read_only(rmfname)
1872  self.prot = IMP.rmf.create_hierarchies(rh, self.model)[0]
1874  IMP.rmf.link_hierarchies(rh, [self.prot])
1875  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1876  del rh
1877  for p in self.prot.get_children():
1878  self.create_component(p.get_name())
1879  self.hier_dict[p.get_name()] = p
1880  '''
1881  still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1882  save floppy bodies in self.floppy_bodies
1883  get the connectivity restraints
1884  '''
1885 
1886  def set_rigid_body_from_hierarchies(self, hiers, particles=None):
1887  '''
1888  Construct a rigid body from hierarchies (and optionally particles).
1889 
1890  @param hiers list of hierarchies to make rigid
1891  @param particles list of particles to add to the rigid body
1892  '''
1893 
1894  if particles is None:
1895  rigid_parts = set()
1896  else:
1897  rigid_parts = set(particles)
1898 
1899  name = ""
1900  print("set_rigid_body_from_hierarchies> setting up a new rigid body")
1901  for hier in hiers:
1902  ps = IMP.atom.get_leaves(hier)
1903 
1904  for p in ps:
1906  rb = IMP.core.RigidMember(p).get_rigid_body()
1907  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1908  else:
1909  rigid_parts.add(p)
1910  name += hier.get_name() + "-"
1911  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1912 
1913  if len(list(rigid_parts)) != 0:
1914  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1915  rb.set_coordinates_are_optimized(True)
1916  rb.set_name(name + "rigid_body")
1917  self.rigid_bodies.append(rb)
1918  return rb
1919 
1920  else:
1921  print("set_rigid_body_from_hierarchies> rigid body could not be setup")
1922 
1923  def set_rigid_bodies(self, subunits):
1924  '''
1925  Construct a rigid body from a list of subunits.
1926 
1927  Each subunit is a tuple that identifies the residue ranges and the
1928  component name (as used in create_component()).
1929 
1930  subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1931  or
1932  [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1933 
1934  example: ["prot1","prot2",("prot3",(1,10))]
1935 
1936  sometimes, we know about structure of an interaction
1937  and here we make such PPIs rigid
1938  '''
1939 
1940  rigid_parts = set()
1941  for s in subunits:
1942  if type(s) == type(tuple()) and len(s) == 2:
1943  sel = IMP.atom.Selection(
1944  self.prot,
1945  molecule=s[0],
1946  residue_indexes=list(range(s[1][0],
1947  s[1][1] + 1)))
1948  if len(sel.get_selected_particles()) == 0:
1949  print("set_rigid_bodies: selected particle does not exist")
1950  for p in sel.get_selected_particles():
1951  # if not p in self.floppy_bodies:
1953  rb = IMP.core.RigidMember(p).get_rigid_body()
1954  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1955  else:
1956  rigid_parts.add(p)
1957 
1958  elif type(s) == type(str()):
1959  sel = IMP.atom.Selection(self.prot, molecule=s)
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():
1963  # if not p in self.floppy_bodies:
1965  rb = IMP.core.RigidMember(p).get_rigid_body()
1966  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1967  else:
1968  rigid_parts.add(p)
1969 
1970  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1971  rb.set_coordinates_are_optimized(True)
1972  rb.set_name(''.join(str(subunits)) + "_rigid_body")
1973  self.rigid_bodies.append(rb)
1974  return rb
1975 
1976  def set_super_rigid_body_from_hierarchies(
1977  self,
1978  hiers,
1979  axis=None,
1980  min_size=1):
1981  # axis is the rotation axis for 2D rotation
1982  super_rigid_xyzs = set()
1983  super_rigid_rbs = set()
1984  name = ""
1985  print("set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
1986 
1987  for hier in hiers:
1988  ps = IMP.atom.get_leaves(hier)
1989  for p in ps:
1991  rb = IMP.core.RigidMember(p).get_rigid_body()
1992  super_rigid_rbs.add(rb)
1993  else:
1994  super_rigid_xyzs.add(p)
1995  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1996  if len(super_rigid_rbs|super_rigid_xyzs) < min_size:
1997  return
1998  if axis is None:
1999  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2000  else:
2001  # these will be 2D rotation SRB or a bond rotamer (axis can be a IMP.algebra.Vector3D or particle Pair)
2002  self.super_rigid_bodies.append(
2003  (super_rigid_xyzs, super_rigid_rbs, axis))
2004 
2005  def fix_rigid_bodies(self, rigid_bodies):
2006  self.fixed_rigid_bodies += rigid_bodies
2007 
2008 
2010  self, hiers, min_length=None, max_length=None, axis=None):
2011  '''
2012  Make a chain of super rigid bodies from a list of hierarchies.
2013 
2014  Takes a linear list of hierarchies (they are supposed to be
2015  sequence-contiguous) and produces a chain of super rigid bodies
2016  with given length range, specified by min_length and max_length.
2017  '''
2018  try:
2019  hiers = IMP.pmi.tools.flatten_list(hiers)
2020  except:
2021  pass
2022  for hs in IMP.pmi.tools.sublist_iterator(hiers, min_length, max_length):
2023  self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
2024 
2025  def set_super_rigid_bodies(self, subunits, coords=None):
2026  super_rigid_xyzs = set()
2027  super_rigid_rbs = set()
2028 
2029  for s in subunits:
2030  if type(s) == type(tuple()) and len(s) == 3:
2031  sel = IMP.atom.Selection(
2032  self.prot,
2033  molecule=s[2],
2034  residue_indexes=list(range(s[0],
2035  s[1] + 1)))
2036  if len(sel.get_selected_particles()) == 0:
2037  print("set_rigid_bodies: selected particle does not exist")
2038  for p in sel.get_selected_particles():
2040  rb = IMP.core.RigidMember(p).get_rigid_body()
2041  super_rigid_rbs.add(rb)
2042  else:
2043  super_rigid_xyzs.add(p)
2044  elif type(s) == type(str()):
2045  sel = IMP.atom.Selection(self.prot, molecule=s)
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():
2049  # if not p in self.floppy_bodies:
2051  rb = IMP.core.RigidMember(p).get_rigid_body()
2052  super_rigid_rbs.add(rb)
2053  else:
2054  super_rigid_xyzs.add(p)
2055  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
2056 
2057  def remove_floppy_bodies_from_component(self, component_name):
2058  '''
2059  Remove leaves of hierarchies from the floppy bodies list based
2060  on the component name
2061  '''
2062  hs=IMP.atom.get_leaves(self.hier_dict[component_name])
2063  ps=[h.get_particle() for h in hs]
2064  for p in self.floppy_bodies:
2065  try:
2066  if p in ps: self.floppy_bodies.remove(p)
2067  if p in hs: self.floppy_bodies.remove(p)
2068  except:
2069  continue
2070 
2071 
2072  def remove_floppy_bodies(self, hierarchies):
2073  '''
2074  Remove leaves of hierarchies from the floppy bodies list.
2075 
2076  Given a list of hierarchies, get the leaves and remove the
2077  corresponding particles from the floppy bodies list. We need this
2078  function because sometimes
2079  we want to constrain the floppy bodies in a rigid body - for instance
2080  when you want to associate a bead with a density particle.
2081  '''
2082  for h in hierarchies:
2083  ps = IMP.atom.get_leaves(h)
2084  for p in ps:
2085  if p in self.floppy_bodies:
2086  print("remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
2087  self.floppy_bodies.remove(p)
2088 
2089 
2090  def set_floppy_bodies(self):
2091  for p in self.floppy_bodies:
2092  name = p.get_name()
2093  p.set_name(name + "_floppy_body")
2095  print("I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
2096  rb = IMP.core.RigidMember(p).get_rigid_body()
2097  try:
2098  rb.set_is_rigid_member(p.get_index(), False)
2099  except:
2100  # some IMP versions still work with that
2101  rb.set_is_rigid_member(p.get_particle_index(), False)
2102  p.set_name(p.get_name() + "_rigid_body_member")
2103 
2104  def set_floppy_bodies_from_hierarchies(self, hiers):
2105  for hier in hiers:
2106  ps = IMP.atom.get_leaves(hier)
2107  for p in ps:
2108  IMP.core.XYZ(p).set_coordinates_are_optimized(True)
2109  self.floppy_bodies.append(p)
2110 
2112  self,
2113  selection_tuples,
2114  resolution=None):
2115  '''
2116  selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
2117  @return the particles
2118  '''
2119  particles = []
2120  print(selection_tuples)
2121  for s in selection_tuples:
2122  ps = IMP.pmi.tools.select_by_tuple(
2123  representation=self, tupleselection=tuple(s),
2124  resolution=None, name_is_ambiguous=False)
2125  particles += ps
2126  return particles
2127 
2128  def get_connected_intra_pairs(self):
2129  return self.connected_intra_pairs
2130 
2131  def set_rigid_bodies_max_trans(self, maxtrans):
2132  self.maxtrans_rb = maxtrans
2133 
2134  def set_rigid_bodies_max_rot(self, maxrot):
2135  self.maxrot_rb = maxrot
2136 
2137  def set_super_rigid_bodies_max_trans(self, maxtrans):
2138  self.maxtrans_srb = maxtrans
2139 
2140  def set_super_rigid_bodies_max_rot(self, maxrot):
2141  self.maxrot_srb = maxrot
2142 
2143  def set_floppy_bodies_max_trans(self, maxtrans):
2144  self.maxtrans_fb = maxtrans
2145 
2146  def set_rigid_bodies_as_fixed(self, rigidbodiesarefixed=True):
2147  '''
2148  Fix rigid bodies in their actual position.
2149  The get_particles_to_sample() function will return
2150  just the floppy bodies (if they are not fixed).
2151  '''
2152  self.rigidbodiesarefixed = rigidbodiesarefixed
2153 
2154  def set_floppy_bodies_as_fixed(self, floppybodiesarefixed=True):
2155  '''
2156  Fix floppy bodies in their actual position.
2157  The get_particles_to_sample() function will return
2158  just the rigid bodies (if they are not fixed).
2159  '''
2160  self.floppybodiesarefixed=floppybodiesarefixed
2161 
2162  def draw_hierarchy_graph(self):
2163  for c in IMP.atom.Hierarchy(self.prot).get_children():
2164  print("Drawing hierarchy graph for " + c.get_name())
2165  IMP.pmi.output.get_graph_from_hierarchy(c)
2166 
2167  def get_geometries(self):
2168  # create segments at the lowest resolution
2169  seggeos = []
2170  for name in self.hier_geometry_pairs:
2171  for pt in self.hier_geometry_pairs[name]:
2172  p1 = pt[0]
2173  p2 = pt[1]
2174  color = pt[2]
2175  try:
2176  clr = IMP.display.get_rgb_color(color)
2177  except:
2178  clr = IMP.display.get_rgb_color(1.0)
2179  coor1 = IMP.core.XYZ(p1).get_coordinates()
2180  coor2 = IMP.core.XYZ(p2).get_coordinates()
2181  seg = IMP.algebra.Segment3D(coor1, coor2)
2182  seggeos.append(IMP.display.SegmentGeometry(seg, clr))
2183  return seggeos
2184 
2185  def setup_bonds(self):
2186  # create segments at the lowest resolution
2187  seggeos = []
2188  for name in self.hier_geometry_pairs:
2189  for pt in self.hier_geometry_pairs[name]:
2190 
2191  p1 = pt[0]
2192  p2 = pt[1]
2193  if not IMP.atom.Bonded.get_is_setup(p1):
2195  if not IMP.atom.Bonded.get_is_setup(p2):
2197 
2200  IMP.atom.Bonded(p1),
2201  IMP.atom.Bonded(p2),
2202  1)
2203 
2204  def show_component_table(self, name):
2205  if name in self.sequence_dict:
2206  lastresn = len(self.sequence_dict[name])
2207  firstresn = 1
2208  else:
2209  residues = self.hier_db.get_residue_numbers(name)
2210  firstresn = min(residues)
2211  lastresn = max(residues)
2212 
2213  for nres in range(firstresn, lastresn + 1):
2214  try:
2215  resolutions = self.hier_db.get_residue_resolutions(name, nres)
2216  ps = []
2217  for r in resolutions:
2218  ps += self.hier_db.get_particles(name, nres, r)
2219  print("%20s %7s" % (name, nres), " ".join(["%20s %7s" % (str(p.get_name()),
2220  str(IMP.pmi.Resolution(p).get_resolution())) for p in ps]))
2221  except:
2222  print("%20s %20s" % (name, nres), "**** not represented ****")
2223 
2224  def draw_hierarchy_composition(self):
2225 
2226  ks = list(self.elements.keys())
2227  ks.sort()
2228 
2229  max = 0
2230  for k in self.elements:
2231  for l in self.elements[k]:
2232  if l[1] > max:
2233  max = l[1]
2234 
2235  for k in ks:
2236  self.draw_component_composition(k, max)
2237 
2238  def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
2239  from matplotlib import pyplot
2240  import matplotlib as mpl
2241  k = name
2242  tmplist = sorted(self.elements[k], key=itemgetter(0))
2243  try:
2244  endres = tmplist[-1][1]
2245  except:
2246  print("draw_component_composition: missing information for component %s" % name)
2247  return
2248  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
2249  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
2250 
2251  # Set the colormap and norm to correspond to the data for which
2252  # the colorbar will be used.
2253  cmap = mpl.cm.cool
2254  norm = mpl.colors.Normalize(vmin=5, vmax=10)
2255  bounds = [1]
2256  colors = []
2257 
2258  for n, l in enumerate(tmplist):
2259  firstres = l[0]
2260  lastres = l[1]
2261  if l[3] != "end":
2262  if bounds[-1] != l[0]:
2263  colors.append("white")
2264  bounds.append(l[0])
2265  if l[3] == "pdb":
2266  colors.append("#99CCFF")
2267  if l[3] == "bead":
2268  colors.append("#FFFF99")
2269  if l[3] == "helix":
2270  colors.append("#33CCCC")
2271  if l[3] != "end":
2272  bounds.append(l[1] + 1)
2273  else:
2274  if l[3] == "pdb":
2275  colors.append("#99CCFF")
2276  if l[3] == "bead":
2277  colors.append("#FFFF99")
2278  if l[3] == "helix":
2279  colors.append("#33CCCC")
2280  if l[3] != "end":
2281  bounds.append(l[1] + 1)
2282  else:
2283  if bounds[-1] - 1 == l[0]:
2284  bounds.pop()
2285  bounds.append(l[0])
2286  else:
2287  colors.append("white")
2288  bounds.append(l[0])
2289 
2290  bounds.append(bounds[-1])
2291  colors.append("white")
2292  cmap = mpl.colors.ListedColormap(colors)
2293  cmap.set_over('0.25')
2294  cmap.set_under('0.75')
2295 
2296  norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
2297  cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
2298  norm=norm,
2299  # to use 'extend', you must
2300  # specify two extra boundaries:
2301  boundaries=bounds,
2302  ticks=bounds, # optional
2303  spacing='proportional',
2304  orientation='horizontal')
2305 
2306  extra_artists = []
2307  npdb = 0
2308 
2309  if draw_pdb_names:
2310  for l in tmplist:
2311  if l[3] == "pdb":
2312  npdb += 1
2313  mid = 1.0 / endres * float(l[0])
2314  # t =ax.text(mid, float(npdb-1)/2.0+1.5, l[2], ha="left", va="center", rotation=0,
2315  # size=10)
2316  # t=ax.annotate(l[0],2)
2317  t = ax.annotate(
2318  l[2], xy=(mid, 1), xycoords='axes fraction',
2319  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords='axes fraction',
2320  arrowprops=dict(arrowstyle="->",
2321  connectionstyle="angle,angleA=0,angleB=90,rad=10"),
2322  )
2323  extra_artists.append(t)
2324 
2325  # set the title of the bar
2326  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
2327  size=20)
2328 
2329  extra_artists.append(title)
2330  # changing the xticks labels
2331  labels = len(bounds) * [" "]
2332  ax.set_xticklabels(labels)
2333  mid = 1.0 / endres * float(bounds[0])
2334  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
2335  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
2336  extra_artists.append(t)
2337  offsets = [0, -0.5, -1.0]
2338  nclashes = 0
2339  for n in range(1, len(bounds)):
2340  if bounds[n] == bounds[n - 1]:
2341  continue
2342  mid = 1.0 / endres * float(bounds[n])
2343  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
2344  nclashes += 1
2345  offset = offsets[nclashes % 3]
2346  else:
2347  nclashes = 0
2348  offset = offsets[0]
2349  if offset > -0.75:
2350  t = ax.annotate(
2351  bounds[n], xy=(mid, 0), xycoords='axes fraction',
2352  xytext=(mid, -0.5 + offset), textcoords='axes fraction')
2353  else:
2354  t = ax.annotate(
2355  bounds[n], xy=(mid, 0), xycoords='axes fraction',
2356  xytext=(mid, -0.5 + offset), textcoords='axes fraction', arrowprops=dict(arrowstyle="-"))
2357  extra_artists.append(t)
2358 
2359  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
2360  # cb2.set_label(k)
2361 
2362  pyplot.savefig(
2363  k + "structure.pdf",
2364  dpi=150,
2365  transparent="True",
2366  bbox_extra_artists=(extra_artists),
2367  bbox_inches='tight')
2368  pyplot.show()
2369 
2370  def draw_coordinates_projection(self):
2371  import matplotlib.pyplot as pp
2372  xpairs = []
2373  ypairs = []
2374  for name in self.hier_geometry_pairs:
2375  for pt in self.hier_geometry_pairs[name]:
2376  p1 = pt[0]
2377  p2 = pt[1]
2378  color = pt[2]
2379  coor1 = IMP.core.XYZ(p1).get_coordinates()
2380  coor2 = IMP.core.XYZ(p2).get_coordinates()
2381  x1 = coor1[0]
2382  x2 = coor2[0]
2383  y1 = coor1[1]
2384  y2 = coor2[1]
2385  xpairs.append([x1, x2])
2386  ypairs.append([y1, y2])
2387  xlist = []
2388  ylist = []
2389  for xends, yends in zip(xpairs, ypairs):
2390  xlist.extend(xends)
2391  xlist.append(None)
2392  ylist.extend(yends)
2393  ylist.append(None)
2394  pp.plot(xlist, ylist, 'b-', alpha=0.1)
2395  pp.show()
2396 
2397  def get_prot_name_from_particle(self, particle):
2398  names = self.get_component_names()
2399  particle0 = particle
2400  name = None
2401  while not name in names:
2402  h = IMP.atom.Hierarchy(particle0).get_parent()
2403  name = h.get_name()
2404  particle0 = h.get_particle()
2405  return name
2406 
2407 
2408  def get_particles_to_sample(self):
2409  # get the list of samplable particles with their type
2410  # and the mover displacement. Everything wrapped in a dictionary,
2411  # to be used by samplers modules
2412  ps = {}
2413 
2414  # remove symmetric particles: they are not sampled
2415  rbtmp = []
2416  fbtmp = []
2417  srbtmp = []
2418  if not self.rigidbodiesarefixed:
2419  for rb in self.rigid_bodies:
2421  if IMP.pmi.Symmetric(rb).get_symmetric() != 1:
2422  rbtmp.append(rb)
2423  else:
2424  if rb not in self.fixed_rigid_bodies:
2425  rbtmp.append(rb)
2426 
2427  if not self.floppybodiesarefixed:
2428  for fb in self.floppy_bodies:
2430  if IMP.pmi.Symmetric(fb).get_symmetric() != 1:
2431  fbtmp.append(fb)
2432  else:
2433  fbtmp.append(fb)
2434 
2435  for srb in self.super_rigid_bodies:
2436  # going to prune the fixed rigid bodies out
2437  # of the super rigid body list
2438  rigid_bodies = list(srb[1])
2439  filtered_rigid_bodies = []
2440  for rb in rigid_bodies:
2441  if rb not in self.fixed_rigid_bodies:
2442  filtered_rigid_bodies.append(rb)
2443  srbtmp.append((srb[0], filtered_rigid_bodies))
2444 
2445  self.rigid_bodies = rbtmp
2446  self.floppy_bodies = fbtmp
2447  self.super_rigid_bodies = srbtmp
2448 
2449  ps["Rigid_Bodies_SimplifiedModel"] = (
2450  self.rigid_bodies,
2451  self.maxtrans_rb,
2452  self.maxrot_rb)
2453  ps["Floppy_Bodies_SimplifiedModel"] = (
2454  self.floppy_bodies,
2455  self.maxtrans_fb)
2456  ps["SR_Bodies_SimplifiedModel"] = (
2457  self.super_rigid_bodies,
2458  self.maxtrans_srb,
2459  self.maxrot_srb)
2460  return ps
2461 
2462  def set_output_level(self, level):
2463  self.output_level = level
2464 
2465  def _evaluate(self, deriv):
2466  """Evaluate the total score of all added restraints"""
2467  r = IMP.pmi.tools.get_restraint_set(self.model)
2468  return r.evaluate(deriv)
2469 
2470  def get_output(self):
2471  output = {}
2472  score = 0.0
2473 
2474  output["SimplifiedModel_Total_Score_" +
2475  self.label] = str(self._evaluate(False))
2476  output["SimplifiedModel_Linker_Score_" +
2477  self.label] = str(self.linker_restraints.unprotected_evaluate(None))
2478  for name in self.sortedsegments_cr_dict:
2479  partialscore = self.sortedsegments_cr_dict[name].evaluate(False)
2480  score += partialscore
2481  output[
2482  "SimplifiedModel_Link_SortedSegments_" +
2483  name +
2484  "_" +
2485  self.label] = str(
2486  partialscore)
2487  partialscore = self.unmodeledregions_cr_dict[name].evaluate(False)
2488  score += partialscore
2489  output[
2490  "SimplifiedModel_Link_UnmodeledRegions_" +
2491  name +
2492  "_" +
2493  self.label] = str(
2494  partialscore)
2495  for rb in self.rigid_body_symmetries:
2496  name=rb.get_name()
2497  output[name +"_" +self.label]=str(rb.get_reference_frame().get_transformation_to())
2498  for name in self.linker_restraints_dict:
2499  output[
2500  name +
2501  "_" +
2502  self.label] = str(
2503  self.linker_restraints_dict[
2504  name].unprotected_evaluate(
2505  None))
2506 
2507  if len(self.reference_structures.keys()) != 0:
2508  rmsds = self.get_all_rmsds()
2509  for label in rmsds:
2510  output[
2511  "SimplifiedModel_" +
2512  label +
2513  "_" +
2514  self.label] = rmsds[
2515  label]
2516 
2517  if self.output_level == "high":
2518  for p in IMP.atom.get_leaves(self.prot):
2519  d = IMP.core.XYZR(p)
2520  output["Coordinates_" +
2521  p.get_name() + "_" + self.label] = str(d)
2522 
2523  output["_TotalScore"] = str(score)
2524  return output
2525 
2526  def get_test_output(self):
2527  # this method is called by test functions and return an enriched output
2528  output = self.get_output()
2529  for n, p in enumerate(self.get_particles_to_sample()):
2530  output["Particle_to_sample_" + str(n)] = str(p)
2531  output["Number_of_particles"] = len(IMP.atom.get_leaves(self.prot))
2532  output["Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2533  output["Number_of_floppy_bodies"] = len(self.floppy_bodies)
2534  output["Number_of_rigid_bodies"] = len(self.rigid_bodies)
2535  output["Number_of_super_bodies"] = len(self.super_rigid_bodies)
2536  output["Selection_resolution_1"] = len(
2537  IMP.pmi.tools.select(self, resolution=1))
2538  output["Selection_resolution_5"] = len(
2539  IMP.pmi.tools.select(self, resolution=5))
2540  output["Selection_resolution_7"] = len(
2541  IMP.pmi.tools.select(self, resolution=5))
2542  output["Selection_resolution_10"] = len(
2543  IMP.pmi.tools.select(self, resolution=10))
2544  output["Selection_resolution_100"] = len(
2545  IMP.pmi.tools.select(self, resolution=100))
2546  output["Selection_All"] = len(IMP.pmi.tools.select(self))
2547  output["Selection_resolution=1"] = len(
2548  IMP.pmi.tools.select(self, resolution=1))
2549  output["Selection_resolution=1,resid=10"] = len(
2550  IMP.pmi.tools.select(self, resolution=1, residue=10))
2551  for resolution in self.hier_resolution:
2552  output["Hier_resolution_dictionary" +
2553  str(resolution)] = len(self.hier_resolution[resolution])
2554  for name in self.hier_dict:
2555  output[
2556  "Selection_resolution=1,resid=10,name=" +
2557  name] = len(
2559  self,
2560  resolution=1,
2561  name=name,
2562  residue=10))
2563  output[
2564  "Selection_resolution=1,resid=10,name=" +
2565  name +
2566  ",ambiguous"] = len(
2568  self,
2569  resolution=1,
2570  name=name,
2571  name_is_ambiguous=True,
2572  residue=10))
2573  output[
2574  "Selection_resolution=1,resid=10,name=" +
2575  name +
2576  ",ambiguous"] = len(
2578  self,
2579  resolution=1,
2580  name=name,
2581  name_is_ambiguous=True,
2582  residue=10))
2583  output[
2584  "Selection_resolution=1,resrange=(10,20),name=" +
2585  name] = len(
2587  self,
2588  resolution=1,
2589  name=name,
2590  first_residue=10,
2591  last_residue=20))
2592  output[
2593  "Selection_resolution=1,resrange=(10,20),name=" +
2594  name +
2595  ",ambiguous"] = len(
2597  self,
2598  resolution=1,
2599  name=name,
2600  name_is_ambiguous=True,
2601  first_residue=10,
2602  last_residue=20))
2603  output[
2604  "Selection_resolution=10,resrange=(10,20),name=" +
2605  name] = len(
2607  self,
2608  resolution=10,
2609  name=name,
2610  first_residue=10,
2611  last_residue=20))
2612  output[
2613  "Selection_resolution=10,resrange=(10,20),name=" +
2614  name +
2615  ",ambiguous"] = len(
2617  self,
2618  resolution=10,
2619  name=name,
2620  name_is_ambiguous=True,
2621  first_residue=10,
2622  last_residue=20))
2623  output[
2624  "Selection_resolution=100,resrange=(10,20),name=" +
2625  name] = len(
2627  self,
2628  resolution=100,
2629  name=name,
2630  first_residue=10,
2631  last_residue=20))
2632  output[
2633  "Selection_resolution=100,resrange=(10,20),name=" +
2634  name +
2635  ",ambiguous"] = len(
2637  self,
2638  resolution=100,
2639  name=name,
2640  name_is_ambiguous=True,
2641  first_residue=10,
2642  last_residue=20))
2643 
2644  return output
def link_components_to_rmf
Load coordinates in the current representation.
Select non water and non hydrogen atoms.
Definition: pdb.h:245
Tools for handling Gaussian Mixture Models.
Definition: gmm_tools.py:1
Add mass to a particle.
Definition: Mass.h:23
double get_volume_from_residue_type(ResidueType rt)
Return an estimate for the volume of a given residue.
def add_metadata
Associate some metadata with this modeling.
Simple 3D transformation class.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:155
A decorator to associate a particle with a part of a protein/DNA/RNA.
Definition: Fragment.h:20
static Gaussian setup_particle(Model *m, ParticleIndex pi)
Definition: core/Gaussian.h:48
void show_molecular_hierarchy(Hierarchy h)
Print out the molecular hierarchy.
double get_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def get_restraint_set
Get a RestraintSet containing all PMI restraints added to the model.
Definition: tools.py:88
Apply a SingletonFunction to a SingletonContainer to maintain an invariant.
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
Set the coordinates of a particle to be a transformed version of a reference.
Definition: core/symmetry.h:76
Store the representations for a system.
Definition: tools.py:895
def shuffle_configuration
Shuffle configuration; used to restart the optimization.
A member of a rigid body, it has internal (local) coordinates.
Definition: rigid_bodies.h:575
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:576
Set of Python classes to create a multi-state, multi-resolution IMP hierarchy.
static Provenanced setup_particle(Model *m, ParticleIndex pi, Provenance p)
Definition: provenance.h:611
static Atom setup_particle(Model *m, ParticleIndex pi, Atom other)
Definition: atom/Atom.h:242
static Fragment setup_particle(Model *m, ParticleIndex pi)
Definition: Fragment.h:63
Upper bound harmonic function (non-zero when feature > mean)
static Symmetric setup_particle(Model *m, ParticleIndex pi, Float symmetric)
Definition: Symmetric.h:44
def set_coordinates_from_rmf
Read and replace coordinates from an RMF file.
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
A decorator for a particle which has bonds.
Select atoms which are selected by both selectors.
Definition: pdb.h:348
Rotation3D get_random_rotation_3d(const Rotation3D &center, double distance)
Pick a rotation at random near the provided one.
double get_mass(ResidueType c)
Get the mass from the residue type.
def create_non_modeled_component
Create a component that isn't used in the modeling.
Color get_rgb_color(double f)
Return the color for f from the RGB color map.
double get_mass_from_number_of_residues(unsigned int num_aa)
Estimate the mass of a protein from the number of amino acids.
Miscellaneous utilities.
Definition: tools.py:1
Provenance create_clone(Provenance p)
Clone provenance (including previous provenance)
def set_rigid_bodies_as_fixed
Fix rigid bodies in their actual position.
def get_particles_from_selection_tuples
selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....
Add symmetric attribute to a particle.
Definition: Symmetric.h:23
double get_ball_radius_from_volume_3d(double volume)
Return the radius of a sphere with a given volume.
Definition: Sphere3D.h:35
def setup_component_sequence_connectivity
Generate restraints between contiguous fragments in the hierarchy.
def create_rotational_symmetry
The copies must not contain rigid bodies.
Add uncertainty to a particle.
Definition: Uncertainty.h:24
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)
Definition: Residue.h:157
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:690
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: XYZ.h:49
def get_prot_name_from_particle
Return the component name provided a particle and a list of names.
Definition: tools.py:1043
static XYZ setup_particle(Model *m, ParticleIndex pi)
Definition: XYZ.h:51
void read_pdb(TextInput input, int model, Hierarchy h)
def get_residue_gaps_in_hierarchy
Return the residue index gaps and contiguous segments in the hierarchy.
Definition: tools.py:655
def remove_floppy_bodies
Remove leaves of hierarchies from the floppy bodies list.
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
static Uncertainty setup_particle(Model *m, ParticleIndex pi, Float uncertainty)
Definition: Uncertainty.h:45
def remove_floppy_bodies_from_component
Remove leaves of hierarchies from the floppy bodies list based on the component name.
A reference frame in 3D.
Add resolution to a particle.
Definition: Resolution.h:24
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.
Definition: RestraintSet.h:36
static Molecule setup_particle(Model *m, ParticleIndex pi)
Definition: Molecule.h:37
Rotation3D get_rotation_about_axis(const Vector3D &axis, double angle)
Generate a Rotation3D object from a rotation around an axis.
Definition: Rotation3D.h:367
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
def create_transformed_component
Create a transformed view of a component.
Select all CB ATOM records.
Definition: pdb.h:91
A Gaussian distribution in 3D.
Definition: Gaussian3D.h:24
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
static Hierarchy setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor children=ParticleIndexesAdaptor())
Create a Hierarchy of level t by adding the needed attributes.
def deprecated_pmi1_object
Mark a PMI1 class as deprecated.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
def set_rigid_bodies
Construct a rigid body from a list of subunits.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Store a list of ParticleIndexes.
def deprecated_method
Python decorator to mark a method as deprecated.
Definition: __init__.py:9703
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Definition: Mass.h:44
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:69
static Bonded setup_particle(Model *m, ParticleIndex pi)
def add_protocol_output
Capture details of the modeling protocol.
The type for a residue.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Symmetric.h:29
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def add_all_atom_densities
Decorates all specified particles as Gaussians directly.
def coarse_hierarchy
Generate all needed coarse grained layers.
def add_component_pdb
Add a component that has an associated 3D structure in a PDB file.
Track creation of a system fragment from a PDB file.
Definition: provenance.h:86
Basic utilities for handling cryo-electron microscopy 3D density maps.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Resolution.h:30
A base class for Keys.
Definition: Key.h:44
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
Definition: Colored.h:62
def set_chain_of_super_rigid_bodies
Make a chain of super rigid bodies from a list of hierarchies.
def add_component_sequence
Add the primary sequence for a single component.
def write_gmm_to_map
write density map from GMM.
Definition: gmm_tools.py:113
Tools for clustering and cluster analysis.
Definition: pmi/Analysis.py:1
static Resolution setup_particle(Model *m, ParticleIndex pi, Float resolution)
Definition: Resolution.h:45
def create_components_from_rmf
still not working.
def get_closest_residue_position
this function works with plain hierarchies, as read from the pdb, no multi-scale hierarchies ...
Definition: tools.py:571
3D rotation class.
Definition: Rotation3D.h:46
Transformation3D get_identity_transformation_3d()
Return a transformation that does not do anything.
Find all nearby pairs by testing all pairs.
Classes for writing output files and processing them.
Definition: output.py:1
Simple implementation of segments in 3D.
Definition: Segment3D.h:24
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:158
A decorator for a residue.
Definition: Residue.h:134
Sampling of the system.
Definition: samplers.py:1
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def add_component_density
Sets up a Gaussian Mixture Model for this component.
Sample using Monte Carlo.
Definition: samplers.py:36
The general base class for IMP exceptions.
Definition: exception.h:49
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
VectorD< 3 > Vector3D
Definition: VectorD.h:395
def get_file_dataset
Get the dataset associated with a filename, or None.
IMP::core::RigidBody create_rigid_body(Hierarchy h)
static Reference setup_particle(Model *m, ParticleIndex pi, ParticleIndexAdaptor reference)
Definition: core/symmetry.h:32
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:302
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual particles of a Model object.
Definition: Particle.h:41
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.
Definition: pdb.h:371
def get_hierarchies_at_given_resolution
Get the hierarchies at the given resolution.
Store info for a chain of a protein.
Definition: Chain.h:61
Applies a PairScore to a Pair.
Definition: PairRestraint.h:29
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:726
Select all CA ATOM records.
Definition: pdb.h:77
Python classes to represent, score, sample and analyze models.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Uncertainty.h:30
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.
A decorator for a rigid body.
Definition: rigid_bodies.h:82
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: core/Gaussian.h:47
A dictionary-like wrapper for reading and storing sequence data.
def write_gmm_to_text
write a list of gaussians to text.
Definition: gmm_tools.py:62
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.
Definition: provenance.h:586
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Definition: Chain.h:81
Hierarchies get_leaves(const Selection &h)
def add_component_necklace
Generates a string of beads with given length.
An exception for an invalid value being passed to IMP.
Definition: exception.h:137
def set_file_dataset
Associate a dataset with a filename.
def set_floppy_bodies_as_fixed
Fix floppy bodies in their actual position.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
static RigidBody setup_particle(Model *m, ParticleIndex pi, ParticleIndexesAdaptor ps)
Definition: rigid_bodies.h:169
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:189
Support for the RMF file format for storing hierarchical molecular data and markup.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: provenance.h:109
def set_rigid_body_from_hierarchies
Construct a rigid body from hierarchies (and optionally particles).
def decorate_gmm_from_text
read the output from write_gmm_to_text, decorate as Gaussian and Mass
Definition: gmm_tools.py:22
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: tools.py:1063
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:711
Transformation3D get_random_local_transformation(Vector3D origin, double max_translation=5., double max_angle_in_rad=0.26)
Get a local transformation.
def check_root
If the root hierarchy does not exist, construct it.
def add_component_beads
add beads to the representation
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
def add_component_hierarchy_clone
Make a copy of a hierarchy and append it to a component.
Harmonic function (symmetric about the mean)
Definition: core/Harmonic.h:24
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27