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