IMP logo
IMP Reference Guide  2.6.1
The Integrative Modeling Platform
Edit on GitHub version 2.6.1
display/displaying_ensembles.py

The script shows a couple experiments with trying to visualize an ensemble of structures. The ensemble is fairly tight on the assembly scale, but there is significant variation between the location and orientation of the individual proteins (which were modeled as rigid bodies). To save space, the models have had their sidechain atoms removed.

1 ## \example display/displaying_ensembles.py
2 # The script shows a couple experiments with trying to visualize an
3 # ensemble of structures. The ensemble is fairly tight on the assembly
4 # scale, but there is significant variation between the location and
5 # orientation of the individual proteins (which were modeled as rigid
6 # bodies). To save space, the models have had their sidechain atoms
7 # removed.
8 
9 from __future__ import print_function
10 import IMP.display
11 import IMP.atom
12 import sys
13 
14 IMP.setup_from_argv(sys.argv,
15  "Experiments with trying to visualize an ensemble of structures")
16 
17 Segment = IMP.algebra.Segment3D
18 Cylinder = IMP.algebra.Cylinder3D
19 
20 # turn off internal checks to speed things up
21 IMP.set_check_level(IMP.USAGE)
22 
23 
24 def read(m, beyond_file):
25  print("reading")
26  hs = []
27  for i in range(0, beyond_file):
28  # create a simplified version for each chain to speed up computations
29  name = IMP.display.get_example_path(
30  "ensemble/aligned-" + str(i) + ".pdb")
31  h = IMP.atom.read_pdb(name, m, IMP.atom.CAlphaPDBSelector())
32  hr = IMP.atom.Hierarchy.setup_particle(IMP.Particle(m))
33  hs.append(hr)
34  for c in IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE):
35  simp = IMP.atom.create_simplified_along_backbone(
36  IMP.atom.Chain(c), 4)
37  hr.add_child(simp)
38  IMP.atom.destroy(h)
39  if i == 0:
40  base = IMP.atom.get_leaves(hr)
41  print(" ", i)
42  return hs
43 
44 
45 def add_markers(h, c, w):
46  """Add markers to a the passed conformation. The marker locations are chosen
47  pretty thoughtlessly and don't really illustrate the technique well."""
48  def add_marker(s, name):
49  g = IMP.core.XYZRGeometry(s.get_selected_particles()[0])
50  g.set_name(name)
51  g.set_color(c)
52  w.add_geometry(g)
53  s = IMP.atom.Selection(h, chain='B', residue_index=317)
54  add_marker(s, "m0")
55  s = IMP.atom.Selection(h, chain='G', residue_index=212)
56  add_marker(s, "m1")
57  s = IMP.atom.Selection(h, chain='I', residue_index=237)
58  add_marker(s, "m2")
59  s = IMP.atom.Selection(h, chain='F', residue_index=101)
60  add_marker(s, "m3")
61 
62 
63 def get_nice_name(h):
64  nm = h.get_name()
65  return nm[nm.find('-') + 1:nm.rfind('.')]
66 
67 
68 def add_axis(h, c, w, chain_colors):
69  """Add a coordinate axis to show the relative orientation of the protein"""
70  for hc in IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE):
71  rb = IMP.core.RigidMember(hc).get_rigid_body()
72  g = IMP.display.ReferenceFrameGeometry(rb.get_reference_frame())
73  g.set_name(get_nice_name(h) + "_orient")
74  if c:
75  g.set_color(c)
76  else:
77  g.set_color(chain_colors[IMP.atom.Chain(hc).get_id()])
78  w.add_geometry(g)
79 
80 
81 def add_skeleton(h, c, r, w, chain_colors):
82  """Show the connectivity skeleton of the conformation to give an idea of
83  how things are laid out"""
84  for hc0 in IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE):
85  for hc1 in IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE):
86  if hc1 <= hc0:
87  continue
88  d = ps.evaluate_index(h.get_model(),
89  (hc0.get_particle_index(),
90  hc1.get_particle_index()), None)
91  if d < 1:
92  d0 = IMP.core.XYZ(hc0)
93  d1 = IMP.core.XYZ(hc1)
94  mp = .5 * (d0.get_coordinates() + d1.get_coordinates())
95  g = IMP.display.CylinderGeometry(
96  Cylinder(Segment(d0.get_coordinates(), mp), r))
97  if c:
98  g.set_color(c)
99  else:
100  g.set_color(chain_colors[IMP.atom.Chain(d0).get_id()])
101  g.set_name(get_nice_name(h) + "_skel")
102  w.add_geometry(g)
103  g = IMP.display.CylinderGeometry(
104  Cylinder(Segment(d1.get_coordinates(), mp), r))
105  if c:
106  g.set_color(c)
107  else:
108  g.set_color(chain_colors[IMP.atom.Chain(d1).get_id()])
109  g.set_name(get_nice_name(h) + "_skel")
110  w.add_geometry(g)
111 
112 IMP.set_log_level(IMP.TERSE)
113 m = IMP.Model()
114 
115 # change to 46 to display all of them
116 hs = read(m, 3)
117 
118 # used to test of two molecules are touching one another
119 ps = IMP.core.KClosePairsPairScore(
120  IMP.core.SphereDistancePairScore(IMP.core.HarmonicUpperBound(10, 1)),
121  IMP.core.LeavesRefiner(IMP.atom.Hierarchy.get_traits()))
122 ps.set_log_level(IMP.SILENT)
123 
124 
125 print("creating rigid bodies")
126 base_chains = {}
127 for hc in IMP.atom.get_by_type(hs[0], IMP.atom.CHAIN_TYPE):
128  c = IMP.atom.Chain(hc)
129  base_chains[c.get_id()] = c
130 
131 for i, h in enumerate(hs):
132  for hc in IMP.atom.get_by_type(h, IMP.atom.CHAIN_TYPE):
133  c = IMP.atom.Chain(hc)
134  if h == hs[0]:
135  crb = IMP.atom.create_rigid_body(hc)
136  else:
137  # make sure the rigid bodies have equivalent defining reference frames
138  # if we just used IMP.atom.create_rigid_body, globular proteins are likely
139  # to have different axis computed when starting in different
140  # orientations
141  crb = IMP.atom.create_compatible_rigid_body(
142  hc, base_chains[c.get_id()])
143  print(" ", i)
144 
145 chains = IMP.atom.get_by_type(hs[0], IMP.atom.CHAIN_TYPE)
146 chains.sort(key = lambda x: IMP.core.XYZ(x).get_x() + IMP.core.XYZ(x).get_y())
147 chain_colors = {}
148 for i, c in enumerate(chains):
149  id = IMP.atom.Chain(c).get_id()
150  #f= i/float(len(chains))
151  color = IMP.display.get_display_color(i)
152  # IMP.display.get_jet_color(f)
153  chain_colors[id] = color
154 
155 w = IMP.display.PymolWriter("markers.pym")
156 add_markers(hs[0], IMP.display.Color(1, 1, 1), w)
157 hso = hs[1:]
158 
159 
160 # sort them spatially so the colors are nicely arranged and allow one to visually connect
161 # the position of one end with that of the other
162 hso.sort(key=lambda h: IMP.core.XYZ(IMP.atom.Selection(h, chain='I',
163  residue_index=237).get_selected_particles()[0]).get_z())
164 print("adding markers", end=' ')
165 for i, h in enumerate(hso):
166  c = IMP.display.get_interpolated_rgb(
167  IMP.display.Color(1, 0, 0), IMP.display.Color(0, 0, 1), i / 50.)
168  add_markers(h, c, w)
169  print(" ", i)
170 w = IMP.display.PymolWriter("axis.pym")
171 print("adding axis", end=' ')
172 add_axis(hs[0], IMP.display.Color(1, 1, 1), w, chain_colors)
173 for i, h in enumerate(hs[1:]):
174  add_axis(h, None, w, chain_colors)
175  print(i, end=' ')
176 
177 w = IMP.display.PymolWriter("skeletons.pym")
178 add_skeleton(hs[0], IMP.display.Color(1, 1, 1), 5, w, chain_colors)
179 print("adding skeleton", end=' ')
180 for i, h in enumerate(hs[1:]):
181  add_skeleton(h, None, 1, w, chain_colors)
182  print(" ", i)