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