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