IMP  2.4.0
The Integrative Modeling Platform
em/local_fitting.py

Shows how to locally refine a fit of a protein inside its density using a MC/CG optimization protocol. This example does not necessarily converges to the global minimum as that may require more optimization steps. If one wishes to use this example as a template for real refinement purposes, please adjust the parameters of the function IMP.em.local_rigid_fitting accordingly.

1 ## \example em/local_fitting.py
2 # Shows how to locally refine a fit of a protein inside
3 # its density using a MC/CG optimization protocol.
4 # This example does not necessarily converges to the global minimum
5 # as that may require more optimization steps.
6 # If one wishes to use this example as a template for real refinement purposes,
7 # please adjust the parameters of the function IMP.em.local_rigid_fitting
8 # accordingly.
9 
10 from __future__ import print_function
11 import IMP.em
12 import IMP.core
13 import IMP.atom
14 import random
15 import math
16 
17 IMP.base.set_log_level(IMP.base.SILENT)
18 IMP.base.set_check_level(IMP.base.NONE)
19 m = IMP.kernel.Model()
20 # 1. setup the input protein
21 # 1.1 select a selector.
22 # using NonWater selector is more accurate but slower
23 # sel=IMP.atom.NonWaterPDBSelector()
24 sel = IMP.atom.CAlphaPDBSelector()
25 # 1.2 read the protein
26 mh = IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"), m, sel)
27 mh_ref = IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"), m, sel)
28 # 1.3 add radius info to each atom, otherwise the resampling would fail.
29 IMP.atom.add_radii(mh)
30 IMP.atom.add_radii(mh_ref)
31 ps = IMP.core.get_leaves(mh)
32 ps_ref = IMP.core.get_leaves(mh_ref)
33 # 2. read the density map of the protein
34 resolution = 8.
35 voxel_size = 1.5
36 dmap = IMP.em.read_map(
37  IMP.em.get_example_path("input.mrc"), IMP.em.MRCReaderWriter())
38 dmap.get_header_writable().set_resolution(resolution)
39 # 3. The protein is now fitted correctly in the density. We can validate
40 # that by making sure that the cross-correlation score is close to 1.
41 
42 # 3.1 generate a sampled density map to the same resolution and spacing as
43 # the target density map. Note that the function we are going to use
44 # (cross_correlation_coefficient) expect to get the same map dimensions as
45 # the target density map.
46 sampled_input_density = IMP.em.SampledDensityMap(dmap.get_header())
47 sampled_input_density.set_particles(ps)
48 sampled_input_density.resample()
49 sampled_input_density.calcRMS()
50 IMP.em.write_map(sampled_input_density, "vv0.mrc", IMP.em.MRCReaderWriter())
51 # 3.2 calculate the cross-correlation score, which should be close to 1
52 best_score = IMP.em.CoarseCC.cross_correlation_coefficient(
53  dmap, sampled_input_density, sampled_input_density.get_header().dmin)
54 print("The CC score of the native transformation is:", best_score)
55 
56 # 4. To denostrate local fitting we locally rotate and translate the
57 # protein and show how we can go back to the correct placement.
58 
59 # 4.1 define a local transformatione
60 translation = IMP.algebra.get_random_vector_in(
61  IMP.algebra.get_unit_bounding_box_3d())
62 axis = IMP.algebra.get_random_vector_on(IMP.algebra.get_unit_sphere_3d())
63 rand_angle = random.uniform(-70. / 180 * math.pi, 70. / 180 * math.pi)
64 r = IMP.algebra.get_rotation_about_axis(axis, rand_angle)
65 local_trans = IMP.algebra.Transformation3D(r, translation)
66 # 4.2 rotate the protein
67 # prot_xyz=IMP.core.XYZs(IMP.core.get_leaves(mh))
68 # for xyz in prot_xyz:
69 # xyz.set_coordinates(local_trans.get_transformed(xyz.get_coordinates()))
70 # 4.2 set the protein as a rigid body
71 IMP.atom.create_rigid_body(mh)
72 prot_rb = IMP.core.RigidMember(IMP.core.get_leaves(mh)[0]).get_rigid_body()
73 # 4.3 apply the trasnformation to the protein
74 IMP.core.transform(prot_rb, local_trans)
75 m.evaluate(False) # to make sure the transformation was applied
76 # 4.4 print the new correlation score, should be lower than before
77 print(len(IMP.core.get_leaves(mh)))
78 IMP.atom.write_pdb(mh, "input2.pdb")
79 sampled_input_density.resample()
80 sampled_input_density.calcRMS()
81 IMP.em.write_map(sampled_input_density, "vv.mrc", IMP.em.MRCReaderWriter())
82 start_score = IMP.em.CoarseCC.cross_correlation_coefficient(
83  dmap, sampled_input_density, sampled_input_density.get_header().dmin)
84 start_rmsd = IMP.atom.get_rmsd(IMP.core.XYZs(ps), IMP.core.XYZs(ps_ref))
85 print("The start score is:", start_score, "with rmsd of:", start_rmsd)
86 # 5. apply local fitting
87 # 5.1 run local fitting
88 print("preforming local refinement, may run for 3-4 minutes")
89 # translate the molecule to the center of the density
90 IMP.core.transform(prot_rb, IMP.algebra.Transformation3D(
91  IMP.algebra.get_identity_rotation_3d(), dmap.get_centroid() - IMP.core.get_centroid(ps)))
92 m.evaluate(False) # to make sure the transformation was applied
93 sampled_input_density.resample()
94 sampled_input_density.calcRMS()
95 rmsd = IMP.atom.get_rmsd(IMP.core.XYZs(ps), IMP.core.XYZs(ps_ref))
96 score2 = IMP.em.CoarseCC.cross_correlation_coefficient(
97  dmap, sampled_input_density, sampled_input_density.get_header().dmin)
98 print("The score after centering is:", score2, "with rmsd of:", rmsd)
99 # IMP.em.local_rigid_fitting_grid_search(
100 # ps,IMP.core.XYZR.get_radius_key(),
101 # IMP.atom.Mass.get_mass_key(),
102 # dmap,fitting_sols)
103 
104 refiner = IMP.core.LeavesRefiner(IMP.atom.Hierarchy.get_traits())
105 fitting_sols = IMP.em.local_rigid_fitting(
106  mh.get_particle(), refiner,
107  IMP.atom.Mass.get_mass_key(),
108  dmap, [], 2, 10, 10)
109 
110 # 5.2 report best result
111 # 5.2.1 transform the protein to the preferred transformation
112 print("The start score is:", start_score, "with rmsd of:", start_rmsd)
113 for i in range(fitting_sols.get_number_of_solutions()):
114  IMP.core.transform(prot_rb, fitting_sols.get_transformation(i))
115  # prot_rb.set_reference_frame(IMP.algebra.ReferenceFrame3D(fitting_sols.get_transformation(i)))
116  m.evaluate(False) # to make sure the transformation was applied
117 # 5.2.2 calc rmsd to native configuration
118  rmsd = IMP.atom.get_rmsd(
119  IMP.core.XYZs(ps), IMP.core.XYZs(IMP.core.get_leaves(mh_ref)))
120  IMP.atom.write_pdb(mh, "temp_" + str(i) + ".pdb")
121  print("Fit with index:", i, " with cc: ", 1. - fitting_sols.get_score(i), " and rmsd to native of:", rmsd)
122  IMP.atom.write_pdb(mh, "sol_" + str(i) + ".pdb")
123  IMP.core.transform(
124  prot_rb, fitting_sols.get_transformation(i).get_inverse())
125 print("done")