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