IMP  2.4.0
The Integrative Modeling Platform
em2d/optimize_em2d_with_montecarlo.py

Example of optimizing an EM2D restraint using Monte Carlo.

1 ## \example em2d/optimize_em2d_with_montecarlo.py
2 # Example of optimizing an EM2D restraint using Monte Carlo.
3 #
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.base
8 import IMP.core
9 import IMP.atom
10 import IMP.em2d
11 import IMP.algebra
12 import IMP.container
13 import random
14 
15 
16 # An Optimizer score to get the values of the statistics after a given set
17 # of evaluations
18 class WriteStatisticsOptimizerScore(IMP.OptimizerState):
19 
20  def __init__(self, m, restraints):
21  IMP.OptimizerState.__init__(self, m, "WriteStats")
22  self.count = 0
23  self.restraints = restraints
24 
25  def update(self):
26  if (self.count != 10):
27  self.count += 1
28  return
29  else:
30  self.count = 0
31  o = self.get_optimizer()
32  m = o.get_model()
33  for r in self.restraints:
34  print(r.get_name(), r.get_last_score())
35  # for i in range(0,m.get_number_of_restraints()):
36  # r=m.get_restraint(i)
37  # print "restraint",r.get_name(),"value",r.evaluate(False)
38 
39 
40 # Get model from PDB file
41 IMP.base.set_log_level(IMP.base.TERSE)
42 m = IMP.kernel.Model()
43 prot = IMP.atom.read_pdb(
44  IMP.em2d.get_example_path("1z5s.pdb"), m, IMP.atom.ATOMPDBSelector())
45 IMP.atom.add_radii(prot)
46 
47 # get the chains
48 chains = IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
49 print("there are", len(chains), "chains in 1z5s.pdb")
50 
51 # set the chains as rigid bodies
52 native_chain_centers = []
53 rigid_bodies = []
54 for c in chains:
55  atoms = IMP.core.get_leaves(c)
56  rbd = IMP.core.RigidBody.setup_particle(c, atoms)
57  rbd.set_coordinates_are_optimized(True)
58  rigid_bodies.append(rbd)
59  print("chain has", rbd.get_number_of_members(), \
60  "atoms", "coordinates: ", rbd.get_coordinates())
61  native_chain_centers.append(rbd.get_coordinates())
62 
63 bb = IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(-25, -40, -60),
64  IMP.algebra.Vector3D(25, 40, 60))
65 # rotate and translate the chains
66 for rbd in rigid_bodies:
67  # rotation
68  rotation = IMP.algebra.get_random_rotation_3d()
69  transformation1 = IMP.algebra.get_rotation_about_point(
70  rbd.get_coordinates(), rotation)
71  # translation
72  transformation2 = IMP.algebra.Transformation3D(
73  IMP.algebra.get_random_vector_in(bb))
74  # Apply
75  final_transformation = IMP.algebra.compose(
76  transformation1, transformation2)
77  IMP.core.transform(rbd, final_transformation)
78 print("Writing transformed assembly")
79 IMP.atom.write_pdb(prot, "1z5s-transformed.pdb")
80 
81 # set distance restraints measusring some distances between rigid bodies
82 # for the solution.
83 d01 = IMP.algebra.get_distance(
84  native_chain_centers[0], native_chain_centers[1])
85 r01 = IMP.core.DistanceRestraint(
86  IMP.core.Harmonic(d01, 1), chains[0], chains[1])
87 r01.set_name("distance 0-1")
88 d12 = IMP.algebra.get_distance(
89  native_chain_centers[1], native_chain_centers[2])
90 r12 = IMP.core.DistanceRestraint(
91  IMP.core.Harmonic(d12, 1), chains[1], chains[2])
92 r12.set_name("distance 1-2")
93 d23 = IMP.algebra.get_distance(
94  native_chain_centers[2], native_chain_centers[3])
95 r23 = IMP.core.DistanceRestraint(
96  IMP.core.Harmonic(d23, 1), chains[2], chains[3])
97 r23.set_name("distance 2-3")
98 d30 = IMP.algebra.get_distance(
99  native_chain_centers[3], native_chain_centers[0])
100 r30 = IMP.core.DistanceRestraint(
101  IMP.core.Harmonic(d30, 1), chains[3], chains[0])
102 r30.set_name("distance 3-0")
103 print("Distances in the solution: d01 =", \
104  d01, "d12 =", d12, "d23 =", d23, "d30 =", d30)
105 
106 # set distance restraints
107 print("adding distance restraints ")
108 for r in [r01, r12, r23, r30]:
109  m.add_restraint(r)
110 print("model has ", m.get_number_of_restraints(), "restraints")
111 
112 # set em2D restraint
113 srw = IMP.em2d.SpiderImageReaderWriter()
114 selection_file = IMP.em2d.get_example_path("all-1z5s-projections.sel")
115 images_to_read_names = [IMP.base.get_relative_path(selection_file, x) for x in
116  IMP.em2d.read_selection_file(selection_file)]
117 em_images = IMP.em2d.read_images(images_to_read_names, srw)
118 print(len(em_images), "images read")
119 
120 em2d_restraint = IMP.em2d.Em2DRestraint(m)
121 apix = 1.5 # sampling rate of the available EM images
122 # resolution at which you want to generate the projections of the model
123 # In principle you want "perfect" projections, so use the highest resolution
124 resolution = 1
125 # Number of projections to use for the initial registration
126 # (coarse registration) to estimate the registration parameters
127 n_projections = 20
128 params = IMP.em2d.Em2DRestraintParameters(apix, resolution, n_projections)
129 
130 # This method (recommended) uses preprocessing of the images and projections
131 # to speed-up the registration
132 params.coarse_registration_method = IMP.em2d.ALIGN2D_PREPROCESSING
133 # use true if you want to save the projections from the model that best
134 # match the Em images
135 params.save_match_images = False
136 
137 score_function = IMP.em2d.EM2DScore()
138 em2d_restraint = IMP.em2d.Em2DRestraint(m)
139 em2d_restraint.setup(score_function, params)
140 em2d_restraint.set_images(em_images)
141 em2d_restraint.set_name("em2d restraint")
142 container = IMP.container.ListSingletonContainer(IMP.core.get_leaves(prot))
143 em2d_restraint.set_particles(container)
144 em2d_restraints_set = IMP.kernel.RestraintSet(m)
145 
146 # The next two lines are commented, because the optimization of the example
147 # is expensive. To run the full example, uncomment them (It can take a few
148 # hours).
149 # em2d_restraints_set.add_restraint(em2d_restraint)
150 # em2d_restraints_set.set_weight(1000) # weight for the em2D restraint
151 
152 print("adding em2d restraint ")
153 m.add_restraint(em2d_restraints_set)
154 # Add all restraints to a model
155 print("model has ", m.get_number_of_restraints(), "restraints")
156 
157 
158 # MONTECARLO OPTIMIZATION
159 s = IMP.core.MonteCarlo(m)
160 # Add movers for the rigid bodies
161 movers = []
162 for rbd in rigid_bodies:
163  movers.append(IMP.core.RigidBodyMover(rbd, 5, 2))
164 s.add_movers(movers)
165 print("MonteCarlo sampler has", s.get_number_of_movers(), "movers")
166 # Add an optimizer state to save intermediate configurations of the hierarchy
167 o_state = IMP.atom.WritePDBOptimizerState(chains, "intermediate-step-%1%.pdb")
168 o_state.set_period(10)
169 s.add_optimizer_state(o_state)
170 
171 ostate2 = WriteStatisticsOptimizerScore(m, m.get_restraints())
172 s.add_optimizer_state(ostate2)
173 
174 # Perform optimization
175 temperatures = [200, 100, 60, 40, 20, 5]
176 # 200 steps recommended for accurate optimization; a smaller number is used
177 # here for demonstration purposes
178 optimization_steps = 10
179 for T in temperatures:
180  s.set_kt(T)
181  s.optimize(optimization_steps)
182 IMP.atom.write_pdb(prot, "solution.pdb")
183 
184 
185 # Check that the optimization achieves distances close to the solution
186 print("*** End optimization ***")
187 new_centers = []
188 for rbd in rigid_bodies:
189  print("chain has", rbd.get_number_of_members(), \
190  "atoms", "coordinates: ", rbd.get_coordinates())
191  new_centers.append(rbd.get_coordinates())
192 
193 d01 = IMP.algebra.get_distance(new_centers[0], new_centers[1])
194 d12 = IMP.algebra.get_distance(new_centers[1], new_centers[2])
195 d23 = IMP.algebra.get_distance(new_centers[2], new_centers[3])
196 d30 = IMP.algebra.get_distance(new_centers[3], new_centers[0])
197 print("Distances at the end of the optimization: d01 =", \
198  d01, "d12 =", d12, "d23 =", d23, "d30 =", d30)