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