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