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