IMP logo
IMP Reference Guide  2.5.0
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.pmi
14 import IMP.isd
15 from math import pi, sqrt
16 
17 from operator import itemgetter
18 import IMP.pmi.tools
19 import IMP.pmi.output
20 import IMP.rmf
21 import RMF
22 import os
23 
24 
25 class Representation(object):
26  # Authors: Peter Cimermancic, Riccardo Pellarin, Charles Greenberg
27 
28  '''
29  Set up the representation of all proteins and nucleic acid macromolecules.
30 
31  Create the molecular hierarchies, representation,
32  sequence connectivity for the various involved proteins and
33  nucleic acid macromolecules:
34 
35  Create a protein, DNA or RNA, represent it as a set of connected balls of appropriate
36  radii and number of residues, PDB at given resolution(s), or ideal helices.
37 
38  How to use the SimplifiedModel class (typical use):
39 
40  see test/test_hierarchy_contruction.py
41 
42  examples:
43 
44  1) Create a chain of helices and flexible parts
45 
46  c_1_119 =self.add_component_necklace("prot1",1,119,20)
47  c_120_131 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(120,131))
48  c_132_138 =self.add_component_beads("prot1",[(132,138)])
49  c_139_156 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(139,156))
50  c_157_174 =self.add_component_beads("prot1",[(157,174)])
51  c_175_182 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(175,182))
52  c_183_194 =self.add_component_beads("prot1",[(183,194)])
53  c_195_216 =self.add_component_ideal_helix("prot1",resolutions=[1,10],resrange=(195,216))
54  c_217_250 =self.add_component_beads("prot1",[(217,250)])
55 
56 
57  self.set_rigid_body_from_hierarchies(c_120_131)
58  self.set_rigid_body_from_hierarchies(c_139_156)
59  self.set_rigid_body_from_hierarchies(c_175_182)
60  self.set_rigid_body_from_hierarchies(c_195_216)
61 
62  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,
63  c_217_250]
64 
65  self.set_chain_of_super_rigid_bodies(clist,2,3)
66 
67  self.set_super_rigid_bodies(["prot1"])
68 
69  '''
70 
71  def __init__(self, m, upperharmonic=True, disorderedlength=True):
72  """Constructor.
73  @param m the model
74  @param upperharmonic This flag uses either harmonic (False)
75  or upperharmonic (True) in the intra-pair
76  connectivity restraint.
77  @param disorderedlength This flag uses either disordered length
78  calculated for random coil peptides (True) or zero
79  surface-to-surface distance between beads (False)
80  as optimal distance for the sequence connectivity
81  restraint.
82  """
83 
84  # this flag uses either harmonic (False) or upperharmonic (True)
85  # in the intra-pair connectivity restraint. Harmonic is used whe you want to
86  # remove the intra-ev term from energy calculations, e.g.:
87  # upperharmonic=False
88  # ip=self.get_connected_intra_pairs()
89  # ev.add_excluded_particle_pairs(ip)
90 
91  self.upperharmonic = upperharmonic
92  self.disorderedlength = disorderedlength
93  self.rigid_bodies = []
94  self.fixed_rigid_bodies = []
95  self.fixed_floppy_bodies = []
96  self.floppy_bodies = []
97  # self.super_rigid_bodies is a list of tuples.
98  # each tuple, corresponding to a particular super rigid body
99  # the tuple is (super_rigid_xyzs,super_rigid_rbs)
100  # where super_rigid_xyzs are the flexible xyz particles
101  # and super_rigid_rbs is the list of rigid bodies.
102  self.super_rigid_bodies = []
103  self.output_level = "low"
104  self.label = "None"
105 
106  self.maxtrans_rb = 2.0
107  self.maxrot_rb = 0.04
108  self.maxtrans_srb = 2.0
109  self.maxrot_srb = 0.2
110  self.rigidbodiesarefixed = False
111  self.maxtrans_fb = 3.0
112  self.resolution = 10.0
113  self.bblenght = 100.0
114  self.kappa = 100.0
115  self.m = m
116 
117  self.representation_is_modified = False
118  self.unmodeledregions_cr_dict = {}
119  self.sortedsegments_cr_dict = {}
121  self.connected_intra_pairs = []
122  self.hier_dict = {}
123  self.color_dict = {}
124  self.sequence_dict = {}
125  self.hier_geometry_pairs = {}
126  self.hier_db = IMP.pmi.tools.HierarchyDatabase()
127  # this dictionary stores the hierarchies by component name and representation type
128  # self.hier_representation[name][representation_type]
129  # where representation type is Res:X, Beads, Densities, Representation,
130  # etc...
131  self.hier_representation = {}
132  self.hier_resolution = {}
133  # reference structures is a dictionary that contains the coordinates of
134  # structures that are used to calculate the rmsd
135  self.reference_structures = {}
136  self.elements = {}
137  self.linker_restraints = IMP.RestraintSet(self.m, "linker_restraints")
138  self.linker_restraints_dict = {}
139  self.threetoone = {'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D',
140  'CYS': 'C', 'GLU': 'E', 'GLN': 'Q', 'GLY': 'G',
141  'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'LYS': 'K',
142  'MET': 'M', 'PHE': 'F', 'PRO': 'P', 'SER': 'S',
143  'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V', 'UNK': 'X'}
144 
145  self.onetothree = dict((v, k) for k, v in self.threetoone.items())
146 
147  self.residuenamekey = IMP.StringKey("ResidueName")
148 
149  def set_label(self, label):
150  self.label = label
151 
152  def create_component(self, name, color=0.0):
154  protein_h.set_name(name)
155  self.hier_dict[name] = protein_h
156  self.hier_representation[name] = {}
157  self.hier_db.add_name(name)
158  self.prot.add_child(protein_h)
159  self.color_dict[name] = color
160  self.elements[name] = []
161 
162  # Deprecation warning
163 
164  @IMP.deprecated_method("2.5", "Use create_component() instead.")
165  def add_component_name(self, *args, **kwargs):
166  self.create_component(*args, **kwargs)
167 
168  def get_component_names(self):
169  return list(self.hier_dict.keys())
170 
171  def add_component_sequence(self, name, filename, id=None, offs=None,
172  format="FASTA"):
173  '''
174  Add the primary sequence for a single component.
175 
176  @param name Human-readable name of the component
177  @param filename Name of the FASTA file
178  @param id Identifier of the sequence in the FASTA file header
179  (if not provided, use `name` instead)
180  '''
181  from Bio import SeqIO
182  handle = open(filename, "rU")
183  record_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta"))
184  handle.close()
185  if id is None:
186  id = name
187  try:
188  length = len(record_dict[id].seq)
189  except KeyError:
190  raise KeyError("id %s not found in fasta file" % id)
191 
192  self.sequence_dict[name] = str(record_dict[id].seq).replace("*", "")
193  if offs is not None:
194  offs_str="-"*offs
195  self.sequence_dict[name]=offs_str+self.sequence_dict[name]
196 
197  self.elements[name].append((length, length, " ", "end"))
198 
199  def autobuild_model(self, name, pdbname, chain,
200  resolutions=None, resrange=None,
201  missingbeadsize=20,
202  beadsize=None, # this is deprecated
203  color=None, pdbresrange=None, offset=0,
204  show=False, isnucleicacid=False,
205  attachbeads=False):
206 
207  if not beadsize is None:
208  IMP.pmi.tools.print_deprecation_warning(
209  "beadsize",
210  "missingbeadsize")
211 
212  self.representation_is_modified = True
213  outhiers = []
214 
215  if color is None:
216  color = self.color_dict[name]
217  else:
218  self.color_dict[name] = color
219 
220  if resolutions is None:
221  resolutions = [1]
222  print("autobuild_model: constructing %s from pdb %s and chain %s" % (name, pdbname, str(chain)))
223 
224  # get the initial and end residues of the pdb
225  t = IMP.atom.read_pdb(pdbname, self.m,
227 
228  # find start and end indexes
229 
230  start = IMP.atom.Residue(
231  t.get_children()[0].get_children()[0]).get_index()
232  end = IMP.atom.Residue(
233  t.get_children()[0].get_children()[-1]).get_index()
234 
235  # check if resrange was defined, otherwise
236  # use the sequence, or the pdb resrange
237 
238  if resrange is None:
239  if name in self.sequence_dict:
240  resrange = (1, len(self.sequence_dict[name]))
241  else:
242  resrange = (start + offset, end + offset)
243  else:
244  start = resrange[0] - offset
245  end = resrange[1] - offset
246 
248  t,
249  resrange[0],
250  resrange[1])
251 
253  t,
254  resrange[0],
255  terminus="N")
257  t,
258  resrange[1],
259  terminus="C")
260 
261  # construct pdb fragments and intervening beads
262  for n, g in enumerate(gaps):
263  first = g[0]
264  last = g[1]
265  if g[2] == "cont":
266  print("autobuild_model: constructing fragment %s from pdb" % (str((first, last))))
267  outhiers += self.add_component_pdb(name, pdbname,
268  chain, resolutions=resolutions,
269  color=color, cacenters=True,
270  resrange=(first, last),
271  offset=offset, isnucleicacid=isnucleicacid)
272  elif g[2] == "gap" and n > 0:
273  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
274  parts = self.hier_db.get_particles_at_closest_resolution(
275  name,
276  first -
277  1,
278  1)
279  xyz = IMP.core.XYZ(parts[0]).get_coordinates()
280  outhiers += self.add_component_necklace(name,
281  first+offset, last+offset, missingbeadsize, incoord=xyz)
282 
283  elif g[2] == "gap" and n == 0:
284  # add pre-beads
285  print("autobuild_model: constructing fragment %s as a bead" % (str((first, last))))
286  outhiers += self.add_component_necklace(name,
287  first+offset, last+offset, missingbeadsize, incoord=xyznter)
288 
289  return outhiers
290 
291  # Deprecation warning
292 
293  @IMP.deprecated_method("2.5", "Use autobuild_model() instead.")
294  def autobuild_pdb_and_intervening_beads(self, *args, **kwargs):
295  r = self.autobuild_model(*args, **kwargs)
296  return r
297 
298  def add_component_pdb(self, name, pdbname, chain, resolutions, color=None,
299  resrange=None, offset=0, cacenters=True, show=False,
300  isnucleicacid=False, readnonwateratoms=False,
301  read_ca_cb_only=False):
302  '''
303  Add a component that has an associated 3D structure in a PDB file.
304 
305  Reads the PDB, and constructs the fragments corresponding to contiguous
306  senquence stretches.
307 
308  @return a list of hierarchies.
309 
310  @param name (string) the name of the component
311  @param pdbname (string) the name of the PDB file
312  @param chain (string or integer) can be either a string (eg, "A")
313  or an integer (eg, 0 or 1) in case you want
314  to get the corresponding chain number in the PDB.
315  @param resolutions (integers) a list of integers that corresponds
316  to the resolutions that have to be generated
317  @param color (float from 0 to 1) the color applied to the
318  hierarchies generated
319  @param resrange (tuple of integers): the residue range to extract
320  from the PDB. It is a tuple (beg,end). If not specified,
321  all residues belonging to the specified chain are read.
322  @param offset (integer) specifies the residue index offset to be
323  applied when reading the PDB (the FASTA sequence is
324  assumed to start from residue 1, so use this parameter
325  if the PDB file does not start at residue 1)
326  @param cacenters (boolean) if True generates resolution=1 beads
327  centered on C-alpha atoms.
328  @param show (boolean) print out the molecular hierarchy at the end.
329  @param isnucleicacid (boolean) use True if you're reading a PDB
330  with nucleic acids.
331  @param readnonwateratoms (boolean) if True fixes some pathological PDB.
332  @param read_ca_cb_only (boolean) if True, only reads CA/CB
333  '''
334 
335  self.representation_is_modified = True
336  if color is None:
337  # if the color is not passed, then get the stored color
338  color = self.color_dict[name]
339  protein_h = self.hier_dict[name]
340  outhiers = []
341 
342  # determine selector
344  if read_ca_cb_only:
345  cacbsel = IMP.atom.OrPDBSelector(
348  sel = IMP.atom.AndPDBSelector(cacbsel, sel)
349  if type(chain) == str:
352  sel)
353  t = IMP.atom.read_pdb(pdbname, self.m, sel)
354 
355  # get the first and last residue
356  start = IMP.atom.Residue(
357  t.get_children()[0].get_children()[0]).get_index()
358  end = IMP.atom.Residue(
359  t.get_children()[0].get_children()[-1]).get_index()
360  c = IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
361  else:
362  t = IMP.atom.read_pdb(pdbname, self.m, sel)
363  c = IMP.atom.Chain(
364  IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[chain])
365 
366  # get the first and last residue
367  start = IMP.atom.Residue(c.get_children()[0]).get_index()
368  end = IMP.atom.Residue(c.get_children()[-1]).get_index()
369  chain = c.get_id()
370 
371  if not resrange is None:
372  if resrange[0] > start:
373  start = resrange[0]
374  if resrange[1] < end:
375  end = resrange[1]
376 
377  if not isnucleicacid:
378  # do what you have to do for proteins
379  sel = IMP.atom.Selection(
380  c,
381  residue_indexes=list(range(
382  start,
383  end + 1)),
384  atom_type=IMP.atom.AT_CA)
385 
386  else:
387  # do what you have to do for nucleic-acids
388  sel = IMP.atom.Selection(
389  c,
390  residue_indexes=list(range(
391  start,
392  end + 1)),
393  atom_type=IMP.atom.AT_P)
394 
395  ps = sel.get_selected_particles()
396  if len(ps) == 0:
397  raise ValueError("%s no residue found in pdb %s chain %s that overlaps with the queried stretch %s-%s" \
398  % (name, pdbname, str(chain), str(resrange[0]), str(resrange[1])))
400 
401  for p in ps:
402  par = IMP.atom.Atom(p).get_parent()
403  ri = IMP.atom.Residue(par).get_index()
404  IMP.atom.Residue(par).set_index(ri + offset)
405  c0.add_child(par)
406  start = start + offset
407  end = end + offset
408 
409  self.elements[name].append(
410  (start, end, pdbname.split("/")[-1] + ":" + chain, "pdb"))
411 
412  outhiers += self.coarse_hierarchy(name, start, end,
413  resolutions, isnucleicacid, c0, protein_h, "pdb", color)
414 
415  if show:
417 
418  del c
419  del c0
420  del t
421 
422  return outhiers
423 
424  def add_component_ideal_helix(
425  self,
426  name,
427  resolutions,
428  resrange,
429  color=None,
430  show=False):
431 
432  self.representation_is_modified = True
433  from math import pi, cos, sin
434 
435  protein_h = self.hier_dict[name]
436  outhiers = []
437  if color is None:
438  color = self.color_dict[name]
439 
440  start = resrange[0]
441  end = resrange[1]
442  self.elements[name].append((start, end, " ", "helix"))
444  for n, res in enumerate(range(start, end + 1)):
445  if name in self.sequence_dict:
446  try:
447  rtstr = self.onetothree[
448  self.sequence_dict[name][res-1]]
449  except:
450  rtstr = "UNK"
451  rt = IMP.atom.ResidueType(rtstr)
452  else:
453  rt = IMP.atom.ResidueType("ALA")
454 
455  # get the residue volume
456  try:
458  # mass=IMP.atom.get_mass_from_residue_type(rt)
459  except IMP.ValueException:
461  IMP.atom.ResidueType("ALA"))
462  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
464 
465  r = IMP.atom.Residue.setup_particle(IMP.Particle(self.m), rt, res)
466  p = IMP.Particle(self.m)
468  x = 2.3 * cos(n * 2 * pi / 3.6)
469  y = 2.3 * sin(n * 2 * pi / 3.6)
470  z = 6.2 / 3.6 / 2 * n * 2 * pi / 3.6
471  d.set_coordinates(IMP.algebra.Vector3D(x, y, z))
472  d.set_radius(radius)
473  # print d
474  a = IMP.atom.Atom.setup_particle(p, IMP.atom.AT_CA)
475  r.add_child(a)
476  c0.add_child(r)
477 
478  outhiers += self.coarse_hierarchy(name, start, end,
479  resolutions, False, c0, protein_h, "helix", color)
480 
481  if show:
483 
484  del c0
485  return outhiers
486 
487  def add_component_beads(self, name, ds, colors=None, incoord=None):
488  """ add beads to the representation
489  @param name the component name
490  @param ds a list of tuples corresponding to the residue ranges
491  of the beads
492  @param colors a list of colors associated to the beads
493  @param incoord the coordinate tuple correspoding to the position
494  of the beads
495  """
496 
497  from math import pi
498  self.representation_is_modified = True
499 
500  protein_h = self.hier_dict[name]
501  outhiers = []
502  if colors is None:
503  colors = [self.color_dict[name]]
504 
505 
506  for n, dss in enumerate(ds):
507  ds_frag = (dss[0], dss[1])
508  self.elements[name].append((dss[0], dss[1], " ", "bead"))
509  prt = IMP.Particle(self.m)
510  if ds_frag[0] == ds_frag[1]:
511  # if the bead represent a single residue
512  if name in self.sequence_dict:
513  try:
514  rtstr = self.onetothree[
515  self.sequence_dict[name][ds_frag[0]-1]]
516  except:
517  rtstr = "UNK"
518  rt = IMP.atom.ResidueType(rtstr)
519  else:
520  rt = IMP.atom.ResidueType("ALA")
521  h = IMP.atom.Residue.setup_particle(prt, rt, ds_frag[0])
522  h.set_name(name + '_%i_bead' % (ds_frag[0]))
523  prt.set_name(name + '_%i_bead' % (ds_frag[0]))
524  resolution = 1
525  else:
527  h.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
528  prt.set_name(name + '_%i-%i_bead' % (ds_frag[0], ds_frag[1]))
529  h.set_residue_indexes(list(range(ds_frag[0], ds_frag[1] + 1)))
530  resolution = len(h.get_residue_indexes())
531  if "Beads" not in self.hier_representation[name]:
533  root.set_name("Beads")
534  self.hier_representation[name]["Beads"] = root
535  protein_h.add_child(root)
536  self.hier_representation[name]["Beads"].add_child(h)
537 
538  for kk in range(ds_frag[0], ds_frag[1] + 1):
539  self.hier_db.add_particles(name, kk, resolution, [h])
540 
541  try:
542  clr = IMP.display.get_rgb_color(colors[n])
543  except:
544  clr = IMP.display.get_rgb_color(colors[0])
545 
547 
548  # decorate particles according to their resolution
549  IMP.pmi.Resolution.setup_particle(prt, resolution)
550 
552  ptem = IMP.core.XYZR(prt)
554  if resolution == 1:
555  try:
557  except IMP.ValueException:
559  IMP.atom.ResidueType("ALA"))
562  ptem.set_radius(radius)
563  else:
564  volume = IMP.atom.get_volume_from_mass(mass)
565  radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
567  ptem.set_radius(radius)
568  try:
569  if not tuple(incoord) is None:
570  ptem.set_coordinates(incoord)
571  except TypeError:
572  pass
575  self.floppy_bodies.append(prt)
576  IMP.core.XYZ(prt).set_coordinates_are_optimized(True)
577  outhiers += [h]
578 
579  return outhiers
580 
581  def add_component_necklace(self, name, begin, end, length, color=None,
582  incoord=None):
583  '''
584  Generates a string of beads with given length.
585  '''
586  self.representation_is_modified = True
587  outhiers = []
588  if color is None:
589  colors=None
590  else:
591  colors=[color]
592  for chunk in IMP.pmi.tools.list_chunks_iterator(range(begin, end + 1), length):
593  outhiers += self.add_component_beads(name,
594  [(chunk[0], chunk[-1])], colors=colors,incoord=incoord)
595  return outhiers
596 
598  self, name, hierarchies=None, selection_tuples=None,
599  particles=None,
600  resolution=0.0, num_components=10,
601  inputfile=None, outputfile=None,
602  outputmap=None,
603  kernel_type=None,
604  covariance_type='full', voxel_size=1.0,
605  out_hier_name='',
606  sampled_points=1000000, num_iter=100,
607  simulation_res=1.0,
608  multiply_by_total_mass=True,
609  transform=None,
610  intermediate_map_fn=None,
611  density_ps_to_copy=None,
612  use_precomputed_gaussians=False):
613  '''
614  Sets up a Gaussian Mixture Model for this component.
615  Can specify input GMM file or it will be computed.
616  @param name component name
617  @param hierarchies set up GMM for some hierarchies
618  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
619  @param particles set up GMM for particles directly
620  @param resolution usual PMI resolution for selecting particles from the hierarchies
621  @param inputfile read the GMM from this file
622  @param outputfile fit and write the GMM to this file (text)
623  @param outputmap after fitting, create GMM density file (mrc)
624  @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
625  @param covariance_type for fitting the GMM. options are 'full', 'diagonal' and 'spherical'
626  @param voxel_size for creating the intermediate density map and output map.
627  lower number increases accuracy but also rasterizing time grows
628  @param out_hier_name name of the output density hierarchy
629  @param sampled_points number of points to sample. more will increase accuracy and fitting time
630  @param num_iter num GMM iterations. more will increase accuracy and fitting time
631  @param multiply_by_total_mass multiply the weights of the GMM by this value (only works on creation!)
632  @param transform for input file only, apply a transformation (eg for multiple copies same GMM)
633  @param intermediate_map_fn for debugging, this will write the itermediate (simulated) map
634  @param density_ps_to_copy in case you already created the appropriate GMM (eg, for beads)
635  @param use_precomputed_gaussians Set this flag and pass fragments - will use roughly spherical Gaussian setup
636  '''
637  import numpy as np
638  import sys
639  import IMP.em
640  import IMP.isd.gmm_tools
641 
642  # prepare output
643  self.representation_is_modified = True
644  out_hier = []
645  protein_h = self.hier_dict[name]
646  if "Densities" not in self.hier_representation[name]:
648  root.set_name("Densities")
649  self.hier_representation[name]["Densities"] = root
650  protein_h.add_child(root)
651 
652  # gather passed particles
653  fragment_particles = []
654  if not particles is None:
655  fragment_particles += particles
656  if not hierarchies is None:
657  fragment_particles += IMP.pmi.tools.select(
658  self, resolution=resolution,
659  hierarchies=hierarchies)
660  if not selection_tuples is None:
661  for st in selection_tuples:
662  fragment_particles += IMP.pmi.tools.select_by_tuple(
663  self, tupleselection=st,
664  resolution=resolution,
665  name_is_ambiguous=False)
666 
667  # compute or read gaussians
668  density_particles = []
669  if inputfile:
671  inputfile, density_particles,
672  self.m, transform)
673  elif density_ps_to_copy:
674  for ip in density_ps_to_copy:
675  p = IMP.Particle(self.m)
676  shape = IMP.core.Gaussian(ip).get_gaussian()
677  mass = IMP.atom.Mass(ip).get_mass()
680  density_particles.append(p)
681  elif use_precomputed_gaussians:
682  if len(fragment_particles) == 0:
683  print("add_component_density: no particle was selected")
684  return out_hier
685  for p in fragment_particles:
686  if not (IMP.atom.Fragment.get_is_setup(self.m,p.get_particle_index()) and
687  IMP.core.XYZ.get_is_setup(self.m,p.get_particle_index())):
688  raise Exception("The particles you selected must be Fragments and XYZs")
689  nres=len(IMP.atom.Fragment(self.m,p.get_particle_index()).get_residue_indexes())
690  pos=IMP.core.XYZ(self.m,p.get_particle_index()).get_coordinates()
691  density_particles=[]
692  try:
693  IMP.isd.get_data_path("beads/bead_%i.txt"%nres)
694  except:
695  raise Exception("We haven't created a bead file for",nres,"residues")
696  transform = IMP.algebra.Transformation3D(pos)
698  IMP.isd.get_data_path("beads/bead_%i.txt"%nres), density_particles,
699  self.m, transform)
700  else:
701  #compute the gaussians here
702  if len(fragment_particles) == 0:
703  print("add_component_density: no particle was selected")
704  return out_hier
705 
706  density_particles = IMP.isd.gmm_tools.sample_and_fit_to_particles(
707  self.m,
708  fragment_particles,
709  num_components,
710  sampled_points,
711  simulation_res,
712  voxel_size,
713  num_iter,
714  covariance_type,
715  multiply_by_total_mass,
716  outputmap,
717  outputfile)
718 
719  # prepare output hierarchy
721  s0.set_name(out_hier_name)
722  self.hier_representation[name]["Densities"].add_child(s0)
723  out_hier.append(s0)
724  for nps, p in enumerate(density_particles):
725  s0.add_child(p)
726  p.set_name(s0.get_name() + '_gaussian_%i' % nps)
727  return out_hier
728 
729  def get_component_density(self, name):
730  return self.hier_representation[name]["Densities"]
731 
732  def add_all_atom_densities(self, name, hierarchies=None,
733  selection_tuples=None,
734  particles=None,
735  resolution=0,
736  output_map=None,
737  voxel_size=1.0):
738  '''Decorates all specified particles as Gaussians directly.
739  @param name component name
740  @param hierarchies set up GMM for some hierarchies
741  @param selection_tuples (list of tuples) example (first_residue,last_residue,component_name)
742  @param particles set up GMM for particles directly
743  @param resolution usual PMI resolution for selecting particles from the hierarchies
744  '''
745 
746  import IMP.em
747  import numpy as np
748  import sys
749  from math import sqrt
750  self.representation_is_modified = False
751 
752  if particles is None:
753  fragment_particles = []
754  else:
755  fragment_particles = particles
756 
757  if not hierarchies is None:
758  fragment_particles += IMP.pmi.tools.select(
759  self, resolution=resolution,
760  hierarchies=hierarchies)
761 
762  if not selection_tuples is None:
763  for st in selection_tuples:
764  fragment_particles += IMP.pmi.tools.select_by_tuple(
765  self, tupleselection=st,
766  resolution=resolution,
767  name_is_ambiguous=False)
768 
769  if len(fragment_particles) == 0:
770  print("add all atom densities: no particle was selected")
771  return
772 
773  # create a spherical gaussian for each particle based on atom type
774  print('setting up all atom gaussians num_particles',len(fragment_particles))
775  for n,p in enumerate(fragment_particles):
776  center=IMP.core.XYZ(p).get_coordinates()
777  rad=IMP.core.XYZR(p).get_radius()
778  mass=IMP.atom.Mass(p).get_mass()
780  center)
783  print('setting up particle',p.get_name(), " as individual gaussian particle")
784 
785  if not output_map is None:
786  print('writing map to', output_map)
788  fragment_particles,
789  output_map,
790  voxel_size)
791 
792  def add_component_hierarchy_clone(self, name, hierarchy):
793  '''
794  Make a copy of a hierarchy and append it to a component.
795  '''
796  outhier = []
797  self.representation_is_modified = True
798  protein_h = self.hier_dict[name]
799  hierclone = IMP.atom.create_clone(hierarchy)
800  hierclone.set_name(hierclone.get_name() + "_clone")
801  protein_h.add_child(hierclone)
802  outhier.append(hierclone)
803 
804  psmain = IMP.atom.get_leaves(hierarchy)
805  psclone = IMP.atom.get_leaves(hierclone)
806 
807  # copying attributes
808  for n, pmain in enumerate(psmain):
809  pclone = psclone[n]
811  resolution = IMP.pmi.Resolution(pmain).get_resolution()
812  IMP.pmi.Resolution.setup_particle(pclone, resolution)
813  for kk in IMP.pmi.tools.get_residue_indexes(pclone):
814  self.hier_db.add_particles(
815  name,
816  kk,
818  [pclone])
819 
821  uncertainty = IMP.pmi.Uncertainty(pmain).get_uncertainty()
822  IMP.pmi.Uncertainty.setup_particle(pclone, uncertainty)
823 
825  symmetric = IMP.pmi.Symmetric(pmain).get_symmetric()
826  IMP.pmi.Symmetric.setup_particle(pclone, symmetric)
827 
828  return outhier
829 
830  def set_coordinates_from_rmf(self, component_name, rmf_file_name,
831  rmf_frame_number,
832  rmf_component_name=None,
833  check_number_particles=True,
834  state_number=0):
835  '''Read and replace coordinates from an RMF file.
836  Replace the coordinates of particles with the same name.
837  It assumes that the RMF and the representation have the particles
838  in the same order.
839  @param component_name Component name
840  @param rmf_component_name Name of the component in the RMF file
841  (if not specified, use `component_name`)
842 
843  '''
844 
845  import IMP.pmi.analysis
846 
847  prots = IMP.pmi.analysis.get_hiers_from_rmf(
848  self.m,
849  rmf_frame_number,
850  rmf_file_name)
851 
852  if not prots:
853  raise ValueError("cannot read hiearchy from rmf")
854 
855  prot=prots[0]
856 
857  # if len(self.rigid_bodies)!=0:
858  # print "set_coordinates_from_rmf: cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies"
859  # exit()
860 
861  allpsrmf = IMP.atom.get_leaves(prot)
862 
863  psrmf = []
864  for p in allpsrmf:
865  (protname, is_a_bead) = IMP.pmi.tools.get_prot_name_from_particle(
866  p, self.hier_dict.keys())
867  if not rmf_component_name is None and protname == rmf_component_name:
868  psrmf.append(p)
869  elif rmf_component_name is None and protname == component_name:
870  psrmf.append(p)
871 
872  psrepr = IMP.atom.get_leaves(self.hier_dict[component_name])
873 
874  if check_number_particles:
875  if len(psrmf) != len(psrepr):
876  raise ValueError("%s cannot proceed the rmf and the representation don't have the same number of particles; "
877  "particles in rmf: %s particles in the representation: %s" % (str(component_name), str(len(psrmf)), str(len(psrepr))))
878 
879  for n, prmf in enumerate(psrmf):
880  prmfname = prmf.get_name()
881  preprname = psrepr[n].get_name()
882  if IMP.core.RigidMember.get_is_setup(psrepr[n]):
883  raise ValueError("component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
885  raise ValueError("component %s cannot proceed if rigid bodies were initialized. Use the function before defining the rigid bodies" % component_name)
886 
887  if prmfname != preprname:
888  print("set_coordinates_from_rmf: WARNING rmf particle and representation particles have not the same name %s %s " % (prmfname, preprname))
890  xyz = IMP.core.XYZ(prmf).get_coordinates()
891  IMP.core.XYZ(psrepr[n]).set_coordinates(xyz)
892  else:
893  print("set_coordinates_from_rmf: WARNING particles are not XYZ decorated %s %s " % (str(IMP.core.XYZ.get_is_setup(prmf)), str(IMP.core.XYZ.get_is_setup(psrepr[n]))))
894 
895  def check_root(self, name, protein_h, resolution):
896  '''
897  If the root hierarchy does not exist, construct it.
898  '''
899  if "Res:" + str(int(resolution)) not in self.hier_representation[name]:
901  root.set_name(name + "_Res:" + str(int(resolution)))
902  self.hier_representation[name][
903  "Res:" + str(int(resolution))] = root
904  protein_h.add_child(root)
905 
906  def coarse_hierarchy(self, name, start, end, resolutions, isnucleicacid,
907  input_hierarchy, protein_h, type, color):
908  '''
909  Generate all needed coarse grained layers.
910 
911  @param name name of the protein
912  @param resolutions list of resolutions
913  @param protein_h root hierarchy
914  @param input_hierarchy hierarchy to coarse grain
915  @param type a string, typically "pdb" or "helix"
916  '''
917  outhiers = []
918 
919  if (1 in resolutions) or (0 in resolutions):
920  # in that case create residues and append atoms
921 
922  if 1 in resolutions:
923  self.check_root(name, protein_h, 1)
925  s1.set_name('%s_%i-%i_%s' % (name, start, end, type))
926  # s1.set_residue_indexes(range(start,end+1))
927  self.hier_representation[name]["Res:1"].add_child(s1)
928  outhiers += [s1]
929  if 0 in resolutions:
930  self.check_root(name, protein_h, 0)
932  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
933  # s0.set_residue_indexes(range(start,end+1))
934  self.hier_representation[name]["Res:0"].add_child(s0)
935  outhiers += [s0]
936 
937  if not isnucleicacid:
938  sel = IMP.atom.Selection(
939  input_hierarchy,
940  atom_type=IMP.atom.AT_CA)
941  else:
942  sel = IMP.atom.Selection(
943  input_hierarchy,
944  atom_type=IMP.atom.AT_P)
945 
946  for p in sel.get_selected_particles():
947  resobject = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
948  if 0 in resolutions:
949  # if you ask for atoms
950  resclone0 = IMP.atom.create_clone(resobject)
951  resindex = IMP.atom.Residue(resclone0).get_index()
952  s0.add_child(resclone0)
953  self.hier_db.add_particles(
954  name,
955  resindex,
956  0,
957  resclone0.get_children())
958 
959  chil = resclone0.get_children()
960  for ch in chil:
962  try:
963  clr = IMP.display.get_rgb_color(color)
964  except:
965  clr = IMP.display.get_rgb_color(1.0)
967 
968  if 1 in resolutions:
969  # else clone the residue
970  resclone1 = IMP.atom.create_clone_one(resobject)
971  resindex = IMP.atom.Residue(resclone1).get_index()
972  s1.add_child(resclone1)
973  self.hier_db.add_particles(name, resindex, 1, [resclone1])
974 
975  rt = IMP.atom.Residue(resclone1).get_residue_type()
976  xyz = IMP.core.XYZ(p).get_coordinates()
977  prt = resclone1.get_particle()
978  prt.set_name('%s_%i_%s' % (name, resindex, type))
979  IMP.core.XYZ.setup_particle(prt).set_coordinates(xyz)
980 
981  try:
983  # mass=IMP.atom.get_mass_from_residue_type(rt)
984  except IMP.ValueException:
986  IMP.atom.ResidueType("ALA"))
987  # mass=IMP.atom.get_mass_from_residue_type(IMP.atom.ResidueType("ALA"))
989  IMP.core.XYZR.setup_particle(prt).set_radius(radius)
994  try:
995  clr = IMP.display.get_rgb_color(color)
996  except:
997  clr = IMP.display.get_rgb_color(1.0)
999 
1000  for r in resolutions:
1001  if r != 0 and r != 1:
1002  self.check_root(name, protein_h, r)
1004  input_hierarchy,
1005  r)
1006 
1007  chil = s.get_children()
1009  s0.set_name('%s_%i-%i_%s' % (name, start, end, type))
1010  # s0.set_residue_indexes(range(start,end+1))
1011  for ch in chil:
1012  s0.add_child(ch)
1013  self.hier_representation[name][
1014  "Res:" + str(int(r))].add_child(s0)
1015  outhiers += [s0]
1016  del s
1017 
1018  for prt in IMP.atom.get_leaves(s0):
1019  ri = IMP.atom.Fragment(prt).get_residue_indexes()
1020  first = ri[0]
1021  last = ri[-1]
1022  if first == last:
1023  prt.set_name('%s_%i_%s' % (name, first, type))
1024  else:
1025  prt.set_name('%s_%i_%i_%s' % (name, first, last, type))
1026  for kk in ri:
1027  self.hier_db.add_particles(name, kk, r, [prt])
1028 
1029  radius = IMP.core.XYZR(prt).get_radius()
1033 
1034  # setting up color for each particle in the
1035  # hierarchy, if colors missing in the colors list set it to
1036  # red
1037  try:
1038  clr = IMP.display.get_rgb_color(color)
1039  except:
1040  colors.append(1.0)
1041  clr = IMP.display.get_rgb_color(colors[pdb_part_count])
1043 
1044  return outhiers
1045 
1047  '''
1048  Get the hierarchies at the given resolution.
1049 
1050  The map between resolution and hierarchies is cached to accelerate
1051  the selection algorithm. The cache is invalidated when the
1052  representation was changed by any add_component_xxx.
1053  '''
1054 
1055  if self.representation_is_modified:
1056  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1057  resolution)
1058  self.hier_resolution[resolution] = rhbr
1059  self.representation_is_modified = False
1060  return rhbr
1061  else:
1062  if resolution in self.hier_resolution:
1063  return self.hier_resolution[resolution]
1064  else:
1065  rhbr = self.hier_db.get_all_root_hierarchies_by_resolution(
1066  resolution)
1067  self.hier_resolution[resolution] = rhbr
1068  return rhbr
1069 
1071  self, max_translation=300., max_rotation=2.0 * pi,
1072  avoidcollision=True, cutoff=10.0, niterations=100,
1073  bounding_box=None,
1074  excluded_rigid_bodies=None,
1075  hierarchies_excluded_from_collision=None):
1076  '''
1077  Shuffle configuration; used to restart the optimization.
1078 
1079  The configuration of the system is initialized by placing each
1080  rigid body and each bead randomly in a box with a side of
1081  `max_translation` angstroms, and far enough from each other to
1082  prevent any steric clashes. The rigid bodies are also randomly rotated.
1083 
1084  @param avoidcollision check if the particle/rigid body was
1085  placed close to another particle; uses the optional
1086  arguments cutoff and niterations
1087  @param bounding_box defined by ((x1,y1,z1),(x2,y2,z2))
1088  '''
1089 
1090  if excluded_rigid_bodies is None:
1091  excluded_rigid_bodies = []
1092  if hierarchies_excluded_from_collision is None:
1093  hierarchies_excluded_from_collision = []
1094 
1095  if len(self.rigid_bodies) == 0:
1096  print("shuffle_configuration: rigid bodies were not intialized")
1097 
1099  gcpf.set_distance(cutoff)
1100  ps = []
1101  hierarchies_excluded_from_collision_indexes = []
1102  for p in IMP.atom.get_leaves(self.prot):
1104  ps.append(p)
1106  # remove the densities particles out of the calculation
1107  hierarchies_excluded_from_collision_indexes += IMP.get_indexes([p])
1108  allparticleindexes = IMP.get_indexes(ps)
1109 
1110  if not bounding_box is None:
1111  ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1112  ub = IMP.algebra.Vector3D(x1, y1, z1)
1113  lb = IMP.algebra.Vector3D(x2, y2, z2)
1114  bb = IMP.algebra.BoundingBox3D(ub, lb)
1115 
1116  for h in hierarchies_excluded_from_collision:
1117  hierarchies_excluded_from_collision_indexes += IMP.get_indexes(IMP.atom.get_leaves(h))
1118 
1119 
1120  allparticleindexes = list(
1121  set(allparticleindexes) - set(hierarchies_excluded_from_collision_indexes))
1122 
1123  print(hierarchies_excluded_from_collision)
1124  print(len(allparticleindexes),len(hierarchies_excluded_from_collision_indexes))
1125 
1126  for rb in self.rigid_bodies:
1127  if rb not in excluded_rigid_bodies:
1128  if avoidcollision:
1129 
1130  rbindexes = rb.get_member_particle_indexes()
1131 
1132  rbindexes = list(
1133  set(rbindexes) - set(hierarchies_excluded_from_collision_indexes))
1134  otherparticleindexes = list(
1135  set(allparticleindexes) - set(rbindexes))
1136 
1137  if len(otherparticleindexes) is None:
1138  continue
1139 
1140  niter = 0
1141  while niter < niterations:
1142  rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1143 
1144  if not bounding_box is None:
1145  # overrides the perturbation
1146  translation = IMP.algebra.get_random_vector_in(bb)
1148  transformation = IMP.algebra.Transformation3D(
1149  rotation,
1150  translation)
1151  else:
1153  rbxyz,
1154  max_translation,
1155  max_rotation)
1156 
1157  IMP.core.transform(rb, transformation)
1158 
1159  if avoidcollision:
1160  self.m.update()
1161  npairs = len(
1162  gcpf.get_close_pairs(
1163  self.m,
1164  otherparticleindexes,
1165  rbindexes))
1166  if npairs == 0:
1167  niter = niterations
1168  else:
1169  niter += 1
1170  print("shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1171  print("shuffle_configuration: rigid body name: " + rb.get_name())
1172  if niter == niterations:
1173  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1174  else:
1175  break
1176 
1177  print('shuffling', len(self.floppy_bodies), 'floppy bodies')
1178  for fb in self.floppy_bodies:
1179  if avoidcollision:
1182  if rm and nrm:
1183  fbindexes = IMP.get_indexes([fb])
1184  otherparticleindexes = list(
1185  set(allparticleindexes) - set(fbindexes))
1186  if len(otherparticleindexes) is None:
1187  continue
1188  else:
1189  continue
1190 
1191  niter = 0
1192  while niter < niterations:
1193  fbxyz = IMP.core.XYZ(fb).get_coordinates()
1194  if not bounding_box is None:
1195  # overrides the perturbation
1196  translation = IMP.algebra.get_random_vector_in(bb)
1197  transformation = IMP.algebra.Transformation3D(translation)
1198  else:
1200  fbxyz,
1201  max_translation,
1202  max_rotation)
1203 
1204  IMP.core.transform(IMP.core.XYZ(fb), transformation)
1205 
1206  if avoidcollision:
1207  self.m.update()
1208  npairs = len(
1209  gcpf.get_close_pairs(
1210  self.m,
1211  otherparticleindexes,
1212  fbindexes))
1213  if npairs == 0:
1214  niter = niterations
1215  else:
1216  niter += 1
1217  print("shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1218  if niter == niterations:
1219  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1220  else:
1221  break
1222 
1223  def set_current_coordinates_as_reference_for_rmsd(self, label="None"):
1224  # getting only coordinates from pdb
1225  ps = IMP.pmi.tools.select(self, resolution=1.0)
1226  # storing the reference coordinates and the particles
1227  self.reference_structures[label] = (
1228  [IMP.core.XYZ(p).get_coordinates() for p in ps],
1229  ps)
1230 
1231  def get_all_rmsds(self):
1232  rmsds = {}
1233 
1234  for label in self.reference_structures:
1235 
1236  current_coordinates = [IMP.core.XYZ(p).get_coordinates()
1237  for p in self.reference_structures[label][1]]
1238  reference_coordinates = self.reference_structures[label][0]
1239  if len(reference_coordinates) != len(current_coordinates):
1240  print("calculate_all_rmsds: reference and actual coordinates are not the same")
1241  continue
1243  current_coordinates,
1244  reference_coordinates)
1245  rmsd_global = IMP.algebra.get_rmsd(
1246  reference_coordinates,
1247  current_coordinates)
1248  # warning: temporary we are calculating the drms, and not the rmsd,
1249  # for the relative distance
1250  rmsd_relative = IMP.atom.get_drms(
1251  reference_coordinates,
1252  current_coordinates)
1253  rmsds[label + "_GlobalRMSD"] = rmsd_global
1254  rmsds[label + "_RelativeDRMS"] = rmsd_relative
1255  return rmsds
1256 
1257  def setup_component_geometry(self, name, color=None, resolution=1.0):
1258  if color is None:
1259  color = self.color_dict[name]
1260  # this function stores all particle pairs
1261  # ordered by residue number, to be used
1262  # to construct backbone traces
1263  self.hier_geometry_pairs[name] = []
1264  protein_h = self.hier_dict[name]
1265  pbr = IMP.pmi.tools.select(self, name=name, resolution=resolution)
1266  pbr = IMP.pmi.tools.sort_by_residues(pbr)
1267 
1268  for n in range(len(pbr) - 1):
1269  self.hier_geometry_pairs[name].append((pbr[n], pbr[n + 1], color))
1270 
1272  self, name, resolution=10, scale=1.0):
1273  '''
1274  Generate restraints between contiguous fragments in the hierarchy.
1275  The linkers are generated at resolution 10 by default.
1276 
1277  '''
1278 
1279  unmodeledregions_cr = IMP.RestraintSet(self.m, "unmodeledregions")
1280  sortedsegments_cr = IMP.RestraintSet(self.m, "sortedsegments")
1281 
1282  protein_h = self.hier_dict[name]
1283  SortedSegments = []
1284  frs = self.hier_db.get_preroot_fragments_by_resolution(
1285  name,
1286  resolution)
1287 
1288  for fr in frs:
1289  try:
1290  start = fr.get_children()[0]
1291  except:
1292  start = fr
1293 
1294  try:
1295  end = fr.get_children()[-1]
1296  except:
1297  end = fr
1298 
1299  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
1300  endres = IMP.pmi.tools.get_residue_indexes(end)[-1]
1301  SortedSegments.append((start, end, startres))
1302  SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1303 
1304  # connect the particles
1305  for x in range(len(SortedSegments) - 1):
1306  last = SortedSegments[x][1]
1307  first = SortedSegments[x + 1][0]
1308 
1309  nreslast = len(IMP.pmi.tools.get_residue_indexes(last))
1310  lastresn = IMP.pmi.tools.get_residue_indexes(last)[-1]
1311  nresfirst = len(IMP.pmi.tools.get_residue_indexes(first))
1312  firstresn = IMP.pmi.tools.get_residue_indexes(first)[0]
1313 
1314  residuegap = firstresn - lastresn - 1
1315  if self.disorderedlength and (nreslast / 2 + nresfirst / 2 + residuegap) > 20.0:
1316  # calculate the distance between the sphere centers using Kohn
1317  # PNAS 2004
1318  optdist = sqrt(5 / 3) * 1.93 * \
1319  (nreslast / 2 + nresfirst / 2 + residuegap) ** 0.6
1320  # optdist2=sqrt(5/3)*1.93*((nreslast)**0.6+(nresfirst)**0.6)/2
1321  if self.upperharmonic:
1322  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1323  else:
1324  hu = IMP.core.Harmonic(optdist, self.kappa)
1325  dps = IMP.core.DistancePairScore(hu)
1326  else: # default
1327  optdist = (0.0 + (float(residuegap) + 1.0) * 3.6) * scale
1328  if self.upperharmonic: # default
1329  hu = IMP.core.HarmonicUpperBound(optdist, self.kappa)
1330  else:
1331  hu = IMP.core.Harmonic(optdist, self.kappa)
1333 
1334  pt0 = last.get_particle()
1335  pt1 = first.get_particle()
1336  r = IMP.core.PairRestraint(self.m, dps, (pt0.get_index(), pt1.get_index()))
1337 
1338  print("Adding sequence connectivity restraint between", pt0.get_name(), " and ", pt1.get_name(), 'of distance', optdist)
1339  sortedsegments_cr.add_restraint(r)
1340  self.linker_restraints_dict[
1341  "LinkerRestraint-" + pt0.get_name() + "-" + pt1.get_name()] = r
1342  self.connected_intra_pairs.append((pt0, pt1))
1343  self.connected_intra_pairs.append((pt1, pt0))
1344 
1345  IMP.pmi.tools.add_restraint_to_model(self.m, sortedsegments_cr)
1346  IMP.pmi.tools.add_restraint_to_model(self.m, unmodeledregions_cr)
1347  self.linker_restraints.add_restraint(sortedsegments_cr)
1348  self.linker_restraints.add_restraint(unmodeledregions_cr)
1349  self.sortedsegments_cr_dict[name] = sortedsegments_cr
1350  self.unmodeledregions_cr_dict[name] = unmodeledregions_cr
1351 
1352  def optimize_floppy_bodies(self, nsteps, temperature=1.0):
1353  import IMP.pmi.samplers
1354  pts = IMP.pmi.tools.ParticleToSampleList()
1355  for n, fb in enumerate(self.floppy_bodies):
1356  pts.add_particle(fb, "Floppy_Bodies", 1.0, "Floppy_Body_" + str(n))
1357  if len(pts.get_particles_to_sample()) > 0:
1358  mc = IMP.pmi.samplers.MonteCarlo(self.m, [pts], temperature)
1359  print("optimize_floppy_bodies: optimizing %i floppy bodies" % len(self.floppy_bodies))
1360  mc.optimize(nsteps)
1361  else:
1362  print("optimize_floppy_bodies: no particle to optimize")
1363 
1364  def create_rotational_symmetry(self, maincopy, copies):
1365  from math import pi
1366  self.representation_is_modified = True
1367  ncopies = len(copies) + 1
1368 
1369  sel = IMP.atom.Selection(self.prot, molecule=maincopy)
1370  mainparticles = sel.get_selected_particles()
1371 
1372  for k in range(len(copies)):
1374  IMP.algebra.Vector3D(0, 0, 1), 2 * pi / ncopies * (k + 1))
1375  sm = IMP.core.TransformationSymmetry(rotation3D)
1376 
1377  sel = IMP.atom.Selection(self.prot, molecule=copies[k])
1378  copyparticles = sel.get_selected_particles()
1379 
1380  mainpurged = []
1381  copypurged = []
1382  for n, p in enumerate(mainparticles):
1383 
1384  pc = copyparticles[n]
1385 
1386  mainpurged.append(p)
1388 
1389  copypurged.append(pc)
1391 
1393  for n, p in enumerate(mainpurged):
1394 
1395  pc = copypurged[n]
1396  print("setting " + p.get_name() + " as reference for " + pc.get_name())
1397 
1399  lc.add_particle(pc)
1400 
1401  c = IMP.container.SingletonsConstraint(sm, None, lc)
1402  self.m.add_score_state(c)
1403 
1404  self.m.update()
1405 
1406  def create_amyloid_fibril_symmetry(self, maincopy, axial_copies,
1407  longitudinal_copies, axis=(0, 0, 1), translation_value=4.8):
1408 
1409  from math import pi
1410  self.representation_is_modified = True
1411 
1412  outhiers = []
1413  protein_h = self.hier_dict[maincopy]
1414  mainparts = IMP.atom.get_leaves(protein_h)
1415 
1416  for ilc in range(-longitudinal_copies, longitudinal_copies + 1):
1417  for iac in range(axial_copies):
1418  copyname = maincopy + "_a" + str(ilc) + "_l" + str(iac)
1419  self.add_component_name(copyname, 0.0)
1420  for hier in protein_h.get_children():
1421  self.add_component_hierarchy_clone(copyname, hier)
1422  copyhier = self.hier_dict[copyname]
1423  outhiers.append(copyhier)
1424  copyparts = IMP.atom.get_leaves(copyhier)
1426  IMP.algebra.Vector3D(axis),
1427  2 * pi / axial_copies * (float(iac)))
1428  translation_vector = tuple(
1429  [translation_value * float(ilc) * x for x in axis])
1430  print(translation_vector)
1431  translation = IMP.algebra.Vector3D(translation_vector)
1433  IMP.algebra.Transformation3D(rotation3D, translation))
1435  for n, p in enumerate(mainparts):
1436  pc = copyparts[n]
1439  if not IMP.pmi.Symmetric.get_is_setup(pc):
1442  lc.add_particle(pc)
1443  c = IMP.container.SingletonsConstraint(sm, None, lc)
1444  self.m.add_score_state(c)
1445  self.m.update()
1446  return outhiers
1447 
1448  def link_components_to_rmf(self, rmfname, frameindex):
1449  '''
1450  Load coordinates in the current representation.
1451  This should be done only if the hierarchy self.prot is identical
1452  to the one as stored in the rmf i.e. all components were added.
1453  '''
1454  import IMP.rmf
1455  import RMF
1456 
1457  rh = RMF.open_rmf_file_read_only(rmfname)
1458  IMP.rmf.link_hierarchies(rh, [self.prot])
1459  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1460  del rh
1461 
1462  def create_components_from_rmf(self, rmfname, frameindex):
1463  '''
1464  still not working.
1465  create the representation (i.e. hierarchies) from the rmf file.
1466  it will be stored in self.prot, which will be overwritten.
1467  load the coordinates from the rmf file at frameindex.
1468  '''
1469  rh = RMF.open_rmf_file(rmfname)
1470  self.prot = IMP.rmf.create_hierarchies(rh, self.m)[0]
1472  IMP.rmf.link_hierarchies(rh, [self.prot])
1473  IMP.rmf.load_frame(rh, RMF.FrameID(frameindex))
1474  del rh
1475  for p in self.prot.get_children():
1476  self.add_component_name(p.get_name())
1477  self.hier_dict[p.get_name()] = p
1478  '''
1479  still missing: save rigid bodies contained in the rmf in self.rigid_bodies
1480  save floppy bodies in self.floppy_bodies
1481  get the connectivity restraints
1482  '''
1483 
1484  def set_rigid_body_from_hierarchies(self, hiers, particles=None):
1485  '''
1486  Construct a rigid body from hierarchies (and optionally particles).
1487 
1488  @param hiers list of hierarchies to make rigid
1489  @param particles list of particles to add to the rigid body
1490  '''
1491 
1492  if particles is None:
1493  rigid_parts = set()
1494  else:
1495  rigid_parts = set(particles)
1496 
1497  name = ""
1498  print("set_rigid_body_from_hierarchies> setting up a new rigid body")
1499  for hier in hiers:
1500  ps = IMP.atom.get_leaves(hier)
1501 
1502  for p in ps:
1504  rb = IMP.core.RigidMember(p).get_rigid_body()
1505  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1506  else:
1507  rigid_parts.add(p)
1508  name += hier.get_name() + "-"
1509  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1510 
1511  if len(list(rigid_parts)) != 0:
1512  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1513  rb.set_coordinates_are_optimized(True)
1514  rb.set_name(name + "rigid_body")
1515  self.rigid_bodies.append(rb)
1516  return rb
1517 
1518  else:
1519  print("set_rigid_body_from_hierarchies> rigid body could not be setup")
1520 
1521  def set_rigid_bodies(self, subunits):
1522  '''
1523  Construct a rigid body from a list of subunits.
1524 
1525  Each subunit is a tuple that identifies the residue ranges and the
1526  component name (as used in create_component()).
1527 
1528  subunits: [(name_1,(first_residue_1,last_residue_1)),(name_2,(first_residue_2,last_residue_2)),.....]
1529  or
1530  [name_1,name_2,(name_3,(first_residue_3,last_residue_3)),.....]
1531 
1532  example: ["prot1","prot2",("prot3",(1,10))]
1533 
1534  sometimes, we know about structure of an interaction
1535  and here we make such PPIs rigid
1536  '''
1537 
1538  rigid_parts = set()
1539  for s in subunits:
1540  if type(s) == type(tuple()) and len(s) == 2:
1541  sel = IMP.atom.Selection(
1542  self.prot,
1543  molecule=s[0],
1544  residue_indexes=list(range(s[1][0],
1545  s[1][1] + 1)))
1546  if len(sel.get_selected_particles()) == 0:
1547  print("set_rigid_bodies: selected particle does not exists")
1548  for p in sel.get_selected_particles():
1549  # if not p in self.floppy_bodies:
1551  rb = IMP.core.RigidMember(p).get_rigid_body()
1552  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1553  else:
1554  rigid_parts.add(p)
1555 
1556  elif type(s) == type(str()):
1557  sel = IMP.atom.Selection(self.prot, molecule=s)
1558  if len(sel.get_selected_particles()) == 0:
1559  print("set_rigid_bodies: selected particle does not exists")
1560  for p in sel.get_selected_particles():
1561  # if not p in self.floppy_bodies:
1563  rb = IMP.core.RigidMember(p).get_rigid_body()
1564  print("set_rigid_body_from_hierarchies> WARNING particle %s already belongs to rigid body %s" % (p.get_name(), rb.get_name()))
1565  else:
1566  rigid_parts.add(p)
1567 
1568  rb = IMP.atom.create_rigid_body(list(rigid_parts))
1569  rb.set_coordinates_are_optimized(True)
1570  rb.set_name(''.join(str(subunits)) + "_rigid_body")
1571  self.rigid_bodies.append(rb)
1572  return rb
1573 
1574  def set_super_rigid_body_from_hierarchies(
1575  self,
1576  hiers,
1577  axis=None,
1578  min_size=0):
1579  # axis is the rotation axis for 2D rotation
1580  super_rigid_xyzs = set()
1581  super_rigid_rbs = set()
1582  name = ""
1583  print("set_super_rigid_body_from_hierarchies> setting up a new SUPER rigid body")
1584 
1585  for hier in hiers:
1586  ps = IMP.atom.get_leaves(hier)
1587  for p in ps:
1589  rb = IMP.core.RigidMember(p).get_rigid_body()
1590  super_rigid_rbs.add(rb)
1591  else:
1592  super_rigid_xyzs.add(p)
1593  print("set_rigid_body_from_hierarchies> adding %s to the rigid body" % hier.get_name())
1594  if len(super_rigid_rbs) < min_size:
1595  return
1596  if axis is None:
1597  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
1598  else:
1599  # these will be 2D rotation SRB
1600  self.super_rigid_bodies.append(
1601  (super_rigid_xyzs, super_rigid_rbs, axis))
1602 
1603  def fix_rigid_bodies(self, rigid_bodies):
1604  self.fixed_rigid_bodies += rigid_bodies
1605 
1606 
1608  self, hiers, min_length=None, max_length=None, axis=None):
1609  '''
1610  Make a chain of super rigid bodies from a list of hierarchies.
1611 
1612  Takes a linear list of hierarchies (they are supposed to be
1613  sequence-contiguous) and produces a chain of super rigid bodies
1614  with given length range, specified by min_length and max_length.
1615  '''
1616  try:
1617  hiers = IMP.pmi.tools.flatten_list(hiers)
1618  except:
1619  pass
1620  for hs in IMP.pmi.tools.sublist_iterator(hiers, min_length, max_length):
1621  self.set_super_rigid_body_from_hierarchies(hs, axis, min_length)
1622 
1623  def set_super_rigid_bodies(self, subunits, coords=None):
1624  super_rigid_xyzs = set()
1625  super_rigid_rbs = set()
1626 
1627  for s in subunits:
1628  if type(s) == type(tuple()) and len(s) == 3:
1629  sel = IMP.atom.Selection(
1630  self.prot,
1631  molecule=s[2],
1632  residue_indexes=list(range(s[0],
1633  s[1] + 1)))
1634  if len(sel.get_selected_particles()) == 0:
1635  print("set_rigid_bodies: selected particle does not exists")
1636  for p in sel.get_selected_particles():
1638  rb = IMP.core.RigidMember(p).get_rigid_body()
1639  super_rigid_rbs.add(rb)
1640  else:
1641  super_rigid_xyzs.add(p)
1642  elif type(s) == type(str()):
1643  sel = IMP.atom.Selection(self.prot, molecule=s)
1644  if len(sel.get_selected_particles()) == 0:
1645  print("set_rigid_bodies: selected particle does not exists")
1646  for p in sel.get_selected_particles():
1647  # if not p in self.floppy_bodies:
1649  rb = IMP.core.RigidMember(p).get_rigid_body()
1650  super_rigid_rbs.add(rb)
1651  else:
1652  super_rigid_xyzs.add(p)
1653  self.super_rigid_bodies.append((super_rigid_xyzs, super_rigid_rbs))
1654 
1655  def remove_floppy_bodies(self, hierarchies):
1656  '''
1657  Remove leaves of hierarchies from the floppy bodies list.
1658 
1659  Given a list of hierarchies, get the leaves and remove the
1660  corresponding particles from the floppy bodies list. We need this
1661  function because sometimes
1662  we want to constrain the floppy bodies in a rigid body - for instance
1663  when you want to associate a bead with a density particle.
1664  '''
1665  for h in hierarchies:
1666  ps = IMP.atom.get_leaves(h)
1667  for p in ps:
1668  if p in self.floppy_bodies:
1669  print("remove_floppy_bodies: removing %s from floppy body list" % p.get_name())
1670  self.floppy_bodies.remove(p)
1671 
1672 
1673  def set_floppy_bodies(self):
1674  for p in self.floppy_bodies:
1675  name = p.get_name()
1676  p.set_name(name + "_floppy_body")
1678  print("I'm trying to make this particle flexible although it was assigned to a rigid body", p.get_name())
1679  rb = IMP.core.RigidMember(p).get_rigid_body()
1680  try:
1681  rb.set_is_rigid_member(p.get_index(), False)
1682  except:
1683  # some IMP versions still work with that
1684  rb.set_is_rigid_member(p.get_particle_index(), False)
1685  p.set_name(p.get_name() + "_rigid_body_member")
1686 
1687  def set_floppy_bodies_from_hierarchies(self, hiers):
1688  for hier in hiers:
1689  ps = IMP.atom.get_leaves(hier)
1690  for p in ps:
1691  IMP.core.XYZ(p).set_coordinates_are_optimized(True)
1692  self.floppy_bodies.append(p)
1693 
1695  self,
1696  selection_tuples,
1697  resolution=None):
1698  '''
1699  selection tuples must be [(r1,r2,"name1"),(r1,r2,"name2"),....]
1700  @return the particles
1701  '''
1702  particles = []
1703  print(selection_tuples)
1704  for s in selection_tuples:
1705  ps = IMP.pmi.tools.select_by_tuple(
1706  representation=self, tupleselection=tuple(s),
1707  resolution=None, name_is_ambiguous=False)
1708  particles += ps
1709  return particles
1710 
1711  def get_connected_intra_pairs(self):
1712  return self.connected_intra_pairs
1713 
1714  def set_rigid_bodies_max_trans(self, maxtrans):
1715  self.maxtrans_rb = maxtrans
1716 
1717  def set_rigid_bodies_max_rot(self, maxrot):
1718  self.maxrot_rb = maxrot
1719 
1720  def set_super_rigid_bodies_max_trans(self, maxtrans):
1721  self.maxtrans_srb = maxtrans
1722 
1723  def set_super_rigid_bodies_max_rot(self, maxrot):
1724  self.maxrot_srb = maxrot
1725 
1726  def set_floppy_bodies_max_trans(self, maxtrans):
1727  self.maxtrans_fb = maxtrans
1728 
1729  def set_rigid_bodies_as_fixed(self, rigidbodiesarefixed=True):
1730  '''
1731  Fix rigid bodies in their actual position.
1732  The get_particles_to_sample() function will return
1733  just the floppy bodies.
1734  '''
1735  self.rigidbodiesarefixed = rigidbodiesarefixed
1736 
1737  def draw_hierarchy_graph(self):
1738  for c in IMP.atom.Hierarchy(self.prot).get_children():
1739  print("Drawing hierarchy graph for " + c.get_name())
1740  IMP.pmi.output.get_graph_from_hierarchy(c)
1741 
1742  def get_geometries(self):
1743  # create segments at the lowest resolution
1744  seggeos = []
1745  for name in self.hier_geometry_pairs:
1746  for pt in self.hier_geometry_pairs[name]:
1747  p1 = pt[0]
1748  p2 = pt[1]
1749  color = pt[2]
1750  try:
1751  clr = IMP.display.get_rgb_color(color)
1752  except:
1753  clr = IMP.display.get_rgb_color(1.0)
1754  coor1 = IMP.core.XYZ(p1).get_coordinates()
1755  coor2 = IMP.core.XYZ(p2).get_coordinates()
1756  seg = IMP.algebra.Segment3D(coor1, coor2)
1757  seggeos.append(IMP.display.SegmentGeometry(seg, clr))
1758  return seggeos
1759 
1760  def setup_bonds(self):
1761  # create segments at the lowest resolution
1762  seggeos = []
1763  for name in self.hier_geometry_pairs:
1764  for pt in self.hier_geometry_pairs[name]:
1765 
1766  p1 = pt[0]
1767  p2 = pt[1]
1768  if not IMP.atom.Bonded.get_is_setup(p1):
1770  if not IMP.atom.Bonded.get_is_setup(p2):
1772 
1774  IMP.atom.Bonded(p1),
1775  IMP.atom.Bonded(p2),
1776  1)
1777 
1778  def show_component_table(self, name):
1779  if name in self.sequence_dict:
1780  lastresn = len(self.sequence_dict[name])
1781  firstresn = 1
1782  else:
1783  residues = self.hier_db.get_residue_numbers(name)
1784  firstresn = min(residues)
1785  lastresn = max(residues)
1786 
1787  for nres in range(firstresn, lastresn + 1):
1788  try:
1789  resolutions = self.hier_db.get_residue_resolutions(name, nres)
1790  ps = []
1791  for r in resolutions:
1792  ps += self.hier_db.get_particles(name, nres, r)
1793  print("%20s %7s" % (name, nres), " ".join(["%20s %7s" % (str(p.get_name()),
1794  str(IMP.pmi.Resolution(p).get_resolution())) for p in ps]))
1795  except:
1796  print("%20s %20s" % (name, nres), "**** not represented ****")
1797 
1798  def draw_hierarchy_composition(self):
1799 
1800  ks = list(self.elements.keys())
1801  ks.sort()
1802 
1803  max = 0
1804  for k in self.elements:
1805  for l in self.elements[k]:
1806  if l[1] > max:
1807  max = l[1]
1808 
1809  for k in ks:
1810  self.draw_component_composition(k, max)
1811 
1812  def draw_component_composition(self, name, max=1000, draw_pdb_names=False):
1813  from matplotlib import pyplot
1814  import matplotlib as mpl
1815  k = name
1816  tmplist = sorted(self.elements[k], key=itemgetter(0))
1817  try:
1818  endres = tmplist[-1][1]
1819  except:
1820  print("draw_component_composition: missing information for component %s" % name)
1821  return
1822  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
1823  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
1824 
1825  # Set the colormap and norm to correspond to the data for which
1826  # the colorbar will be used.
1827  cmap = mpl.cm.cool
1828  norm = mpl.colors.Normalize(vmin=5, vmax=10)
1829  bounds = [1]
1830  colors = []
1831 
1832  for n, l in enumerate(tmplist):
1833  firstres = l[0]
1834  lastres = l[1]
1835  if l[3] != "end":
1836  if bounds[-1] != l[0]:
1837  colors.append("white")
1838  bounds.append(l[0])
1839  if l[3] == "pdb":
1840  colors.append("#99CCFF")
1841  if l[3] == "bead":
1842  colors.append("#FFFF99")
1843  if l[3] == "helix":
1844  colors.append("#33CCCC")
1845  if l[3] != "end":
1846  bounds.append(l[1] + 1)
1847  else:
1848  if l[3] == "pdb":
1849  colors.append("#99CCFF")
1850  if l[3] == "bead":
1851  colors.append("#FFFF99")
1852  if l[3] == "helix":
1853  colors.append("#33CCCC")
1854  if l[3] != "end":
1855  bounds.append(l[1] + 1)
1856  else:
1857  if bounds[-1] - 1 == l[0]:
1858  bounds.pop()
1859  bounds.append(l[0])
1860  else:
1861  colors.append("white")
1862  bounds.append(l[0])
1863 
1864  bounds.append(bounds[-1])
1865  colors.append("white")
1866  cmap = mpl.colors.ListedColormap(colors)
1867  cmap.set_over('0.25')
1868  cmap.set_under('0.75')
1869 
1870  norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
1871  cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
1872  norm=norm,
1873  # to use 'extend', you must
1874  # specify two extra boundaries:
1875  boundaries=bounds,
1876  ticks=bounds, # optional
1877  spacing='proportional',
1878  orientation='horizontal')
1879 
1880  extra_artists = []
1881  npdb = 0
1882 
1883  if draw_pdb_names:
1884  for l in tmplist:
1885  if l[3] == "pdb":
1886  npdb += 1
1887  mid = 1.0 / endres * float(l[0])
1888  # t =ax.text(mid, float(npdb-1)/2.0+1.5, l[2], ha="left", va="center", rotation=0,
1889  # size=10)
1890  # t=ax.annotate(l[0],2)
1891  t = ax.annotate(
1892  l[2], xy=(mid, 1), xycoords='axes fraction',
1893  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords='axes fraction',
1894  arrowprops=dict(arrowstyle="->",
1895  connectionstyle="angle,angleA=0,angleB=90,rad=10"),
1896  )
1897  extra_artists.append(t)
1898 
1899  # set the title of the bar
1900  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
1901  size=20)
1902 
1903  extra_artists.append(title)
1904  # changing the xticks labels
1905  labels = len(bounds) * [" "]
1906  ax.set_xticklabels(labels)
1907  mid = 1.0 / endres * float(bounds[0])
1908  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
1909  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
1910  extra_artists.append(t)
1911  offsets = [0, -0.5, -1.0]
1912  nclashes = 0
1913  for n in range(1, len(bounds)):
1914  if bounds[n] == bounds[n - 1]:
1915  continue
1916  mid = 1.0 / endres * float(bounds[n])
1917  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
1918  nclashes += 1
1919  offset = offsets[nclashes % 3]
1920  else:
1921  nclashes = 0
1922  offset = offsets[0]
1923  if offset > -0.75:
1924  t = ax.annotate(
1925  bounds[n], xy=(mid, 0), xycoords='axes fraction',
1926  xytext=(mid, -0.5 + offset), textcoords='axes fraction')
1927  else:
1928  t = ax.annotate(
1929  bounds[n], xy=(mid, 0), xycoords='axes fraction',
1930  xytext=(mid, -0.5 + offset), textcoords='axes fraction', arrowprops=dict(arrowstyle="-"))
1931  extra_artists.append(t)
1932 
1933  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
1934  # cb2.set_label(k)
1935 
1936  pyplot.savefig(
1937  k + "structure.pdf",
1938  dpi=150,
1939  transparent="True",
1940  bbox_extra_artists=(extra_artists),
1941  bbox_inches='tight')
1942  pyplot.show()
1943 
1944  def draw_coordinates_projection(self):
1945  import matplotlib.pyplot as pp
1946  xpairs = []
1947  ypairs = []
1948  for name in self.hier_geometry_pairs:
1949  for pt in self.hier_geometry_pairs[name]:
1950  p1 = pt[0]
1951  p2 = pt[1]
1952  color = pt[2]
1953  coor1 = IMP.core.XYZ(p1).get_coordinates()
1954  coor2 = IMP.core.XYZ(p2).get_coordinates()
1955  x1 = coor1[0]
1956  x2 = coor2[0]
1957  y1 = coor1[1]
1958  y2 = coor2[1]
1959  xpairs.append([x1, x2])
1960  ypairs.append([y1, y2])
1961  xlist = []
1962  ylist = []
1963  for xends, yends in zip(xpairs, ypairs):
1964  xlist.extend(xends)
1965  xlist.append(None)
1966  ylist.extend(yends)
1967  ylist.append(None)
1968  pp.plot(xlist, ylist, 'b-', alpha=0.1)
1969  pp.show()
1970 
1971  def get_prot_name_from_particle(self, particle):
1972  names = self.get_component_names()
1973  particle0 = particle
1974  name = None
1975  while not name in names:
1976  h = IMP.atom.Hierarchy(particle0).get_parent()
1977  name = h.get_name()
1978  particle0 = h.get_particle()
1979  return name
1980 
1981 
1982  def get_particles_to_sample(self):
1983  # get the list of samplable particles with their type
1984  # and the mover displacement. Everything wrapped in a dictionary,
1985  # to be used by samplers modules
1986  ps = {}
1987 
1988  # remove symmetric particles: they are not sampled
1989  rbtmp = []
1990  fbtmp = []
1991  srbtmp = []
1992  if not self.rigidbodiesarefixed:
1993  for rb in self.rigid_bodies:
1995  if IMP.pmi.Symmetric(rb).get_symmetric() != 1:
1996  rbtmp.append(rb)
1997  else:
1998  if rb not in self.fixed_rigid_bodies:
1999  rbtmp.append(rb)
2000 
2001  for fb in self.floppy_bodies:
2003  if IMP.pmi.Symmetric(fb).get_symmetric() != 1:
2004  fbtmp.append(fb)
2005  else:
2006  fbtmp.append(fb)
2007 
2008  for srb in self.super_rigid_bodies:
2009  # going to prune the fixed rigid bodies out
2010  # of the super rigid body list
2011  rigid_bodies = list(srb[1])
2012  filtered_rigid_bodies = []
2013  for rb in rigid_bodies:
2014  if rb not in self.fixed_rigid_bodies:
2015  filtered_rigid_bodies.append(rb)
2016  srbtmp.append((srb[0], filtered_rigid_bodies))
2017 
2018  self.rigid_bodies = rbtmp
2019  self.floppy_bodies = fbtmp
2020  self.super_rigid_bodies = srbtmp
2021 
2022  ps["Rigid_Bodies_SimplifiedModel"] = (
2023  self.rigid_bodies,
2024  self.maxtrans_rb,
2025  self.maxrot_rb)
2026  ps["Floppy_Bodies_SimplifiedModel"] = (
2027  self.floppy_bodies,
2028  self.maxtrans_fb)
2029  ps["SR_Bodies_SimplifiedModel"] = (
2030  self.super_rigid_bodies,
2031  self.maxtrans_srb,
2032  self.maxrot_srb)
2033  return ps
2034 
2035  def set_output_level(self, level):
2036  self.output_level = level
2037 
2038  def _evaluate(self, deriv):
2039  """Evaluate the total score of all added restraints"""
2041  return r.evaluate(deriv)
2042 
2043  def get_output(self):
2044  output = {}
2045  score = 0.0
2046 
2047  output["SimplifiedModel_Total_Score_" +
2048  self.label] = str(self._evaluate(False))
2049  output["SimplifiedModel_Linker_Score_" +
2050  self.label] = str(self.linker_restraints.unprotected_evaluate(None))
2051  for name in self.sortedsegments_cr_dict:
2052  partialscore = self.sortedsegments_cr_dict[name].evaluate(False)
2053  score += partialscore
2054  output[
2055  "SimplifiedModel_Link_SortedSegments_" +
2056  name +
2057  "_" +
2058  self.label] = str(
2059  partialscore)
2060  partialscore = self.unmodeledregions_cr_dict[name].evaluate(False)
2061  score += partialscore
2062  output[
2063  "SimplifiedModel_Link_UnmodeledRegions_" +
2064  name +
2065  "_" +
2066  self.label] = str(
2067  partialscore)
2068  for name in self.linker_restraints_dict:
2069  output[
2070  name +
2071  "_" +
2072  self.label] = str(
2073  self.linker_restraints_dict[
2074  name].unprotected_evaluate(
2075  None))
2076 
2077  if len(self.reference_structures.keys()) != 0:
2078  rmsds = self.get_all_rmsds()
2079  for label in rmsds:
2080  output[
2081  "SimplifiedModel_" +
2082  label +
2083  "_" +
2084  self.label] = rmsds[
2085  label]
2086 
2087  if self.output_level == "high":
2088  for p in IMP.atom.get_leaves(self.prot):
2089  d = IMP.core.XYZR(p)
2090  output["Coordinates_" +
2091  p.get_name() + "_" + self.label] = str(d)
2092 
2093  output["_TotalScore"] = str(score)
2094  return output
2095 
2096  def get_test_output(self):
2097  # this method is called by test functions and return an enriched output
2098  output = self.get_output()
2099  for n, p in enumerate(self.get_particles_to_sample()):
2100  output["Particle_to_sample_" + str(n)] = str(p)
2101  output["Number_of_particles"] = len(IMP.atom.get_leaves(self.prot))
2102  output["Hierarchy_Dictionary"] = list(self.hier_dict.keys())
2103  output["Number_of_floppy_bodies"] = len(self.floppy_bodies)
2104  output["Number_of_rigid_bodies"] = len(self.rigid_bodies)
2105  output["Number_of_super_bodies"] = len(self.super_rigid_bodies)
2106  output["Selection_resolution_1"] = len(
2107  IMP.pmi.tools.select(self, resolution=1))
2108  output["Selection_resolution_5"] = len(
2109  IMP.pmi.tools.select(self, resolution=5))
2110  output["Selection_resolution_7"] = len(
2111  IMP.pmi.tools.select(self, resolution=5))
2112  output["Selection_resolution_10"] = len(
2113  IMP.pmi.tools.select(self, resolution=10))
2114  output["Selection_resolution_100"] = len(
2115  IMP.pmi.tools.select(self, resolution=100))
2116  output["Selection_All"] = len(IMP.pmi.tools.select(self))
2117  output["Selection_resolution=1"] = len(
2118  IMP.pmi.tools.select(self, resolution=1))
2119  output["Selection_resolution=1,resid=10"] = len(
2120  IMP.pmi.tools.select(self, resolution=1, residue=10))
2121  for resolution in self.hier_resolution:
2122  output["Hier_resolution_dictionary" +
2123  str(resolution)] = len(self.hier_resolution[resolution])
2124  for name in self.hier_dict:
2125  output[
2126  "Selection_resolution=1,resid=10,name=" +
2127  name] = len(
2129  self,
2130  resolution=1,
2131  name=name,
2132  residue=10))
2133  output[
2134  "Selection_resolution=1,resid=10,name=" +
2135  name +
2136  ",ambiguous"] = len(
2138  self,
2139  resolution=1,
2140  name=name,
2141  name_is_ambiguous=True,
2142  residue=10))
2143  output[
2144  "Selection_resolution=1,resid=10,name=" +
2145  name +
2146  ",ambiguous"] = len(
2148  self,
2149  resolution=1,
2150  name=name,
2151  name_is_ambiguous=True,
2152  residue=10))
2153  output[
2154  "Selection_resolution=1,resrange=(10,20),name=" +
2155  name] = len(
2157  self,
2158  resolution=1,
2159  name=name,
2160  first_residue=10,
2161  last_residue=20))
2162  output[
2163  "Selection_resolution=1,resrange=(10,20),name=" +
2164  name +
2165  ",ambiguous"] = len(
2167  self,
2168  resolution=1,
2169  name=name,
2170  name_is_ambiguous=True,
2171  first_residue=10,
2172  last_residue=20))
2173  output[
2174  "Selection_resolution=10,resrange=(10,20),name=" +
2175  name] = len(
2177  self,
2178  resolution=10,
2179  name=name,
2180  first_residue=10,
2181  last_residue=20))
2182  output[
2183  "Selection_resolution=10,resrange=(10,20),name=" +
2184  name +
2185  ",ambiguous"] = len(
2187  self,
2188  resolution=10,
2189  name=name,
2190  name_is_ambiguous=True,
2191  first_residue=10,
2192  last_residue=20))
2193  output[
2194  "Selection_resolution=100,resrange=(10,20),name=" +
2195  name] = len(
2197  self,
2198  resolution=100,
2199  name=name,
2200  first_residue=10,
2201  last_residue=20))
2202  output[
2203  "Selection_resolution=100,resrange=(10,20),name=" +
2204  name +
2205  ",ambiguous"] = len(
2207  self,
2208  resolution=100,
2209  name=name,
2210  name_is_ambiguous=True,
2211  first_residue=10,
2212  last_residue=20))
2213 
2214  return output
2215 
2216 
2217 # Deprecation warning for the old SimplifiedModel class
2218 
2219 @IMP.deprecated_object("2.5", "Use Representation instead.")
2220 class SimplifiedModel(Representation):
2221 
2222  def __init__(self, *args, **kwargs):
2223  Representation.__init__(self, *args, **kwargs)
def link_components_to_rmf
Load coordinates in the current representation.
Select non water and non hydrogen atoms.
Definition: pdb.h:197
def list_chunks_iterator
Yield successive length-sized chunks from a list.
Definition: tools.py:1074
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.
Simple 3D transformation class.
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: 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:34
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:75
Store the representations for a system.
Definition: tools.py:791
def shuffle_configuration
Shuffle configuration; used to restart the optimization.
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:275
Hierarchy create_clone(Hierarchy d)
Clone the Hierarchy.
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.
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
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.
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:469
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:939
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:585
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
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.
Select all CB ATOM records.
Definition: pdb.h:90
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.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
def set_rigid_bodies
Construct a rigid body from a list of subunits.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
Store a list of ParticleIndexes.
def deprecated_method
Python decorator to mark a method as deprecated.
Definition: __init__.py:9383
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Definition: Mass.h:44
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:22
Hierarchies get_by_type(Hierarchy mhd, GetByType t)
Gather all the molecular particles of a certain level in the hierarchy.
static Bonded setup_particle(Model *m, ParticleIndex pi)
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.
Basic utilities for handling cryo-electron microscopy 3D density maps.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Hierarchy create_clone_one(Hierarchy d)
Clone the node in the Hierarchy.
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:46
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:62
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:522
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
def deprecated_object
Python decorator to mark a class as deprecated.
Definition: __init__.py:9367
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:35
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:31
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:236
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual model particles.
Definition: Particle.h:37
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 selector.
Definition: pdb.h:293
def get_hierarchies_at_given_resolution
Get the hierarchies at the given resolution.
Store info for a chain of a protein.
Definition: Chain.h:21
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:662
Select all CA ATOM records.
Definition: pdb.h:76
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.
Set up the representation of all proteins and nucleic acid macromolecules.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Gaussian.h:47
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.
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Definition: Chain.h:41
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
Select hierarchy particles identified by the biological name.
Definition: Selection.h:65
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:143
Support for the RMF file format for storing hierarchical molecular data and markup.
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:959
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:490
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:1055
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