IMP  2.4.0
The Integrative Modeling Platform
display/displaying_ensembles.py

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