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