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