IMP  2.2.1
The Integrative Modeling Platform
local_fitting.py
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
14 import math
15 
16 IMP.base.set_log_level(IMP.base.SILENT)
17 IMP.base.set_check_level(IMP.base.NONE)
18 m = IMP.kernel.Model()
19 # 1. setup the input protein
20 # 1.1 select a selector.
21 # using NonWater selector is more accurate but slower
22 # sel=IMP.atom.NonWaterPDBSelector()
24 # 1.2 read the protein
25 mh = IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"), m, sel)
26 mh_ref = IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"), m, sel)
27 # 1.3 add radius info to each atom, otherwise the resampling would fail.
29 IMP.atom.add_radii(mh_ref)
30 ps = IMP.core.get_leaves(mh)
31 ps_ref = IMP.core.get_leaves(mh_ref)
32 # 2. read the density map of the protein
33 resolution = 8.
34 voxel_size = 1.5
35 dmap = IMP.em.read_map(
37 dmap.get_header_writable().set_resolution(resolution)
38 # 3. The protein is now fitted correctly in the density. We can validate
39 # that by making sure that the cross-correlation score is close to 1.
40 
41 # 3.1 generate a sampled density map to the same resolution and spacing as
42 # the target density map. Note that the function we are going to use
43 # (cross_correlation_coefficient) expect to get the same map dimensions as
44 # the target density map.
45 sampled_input_density = IMP.em.SampledDensityMap(dmap.get_header())
46 sampled_input_density.set_particles(ps)
47 sampled_input_density.resample()
48 sampled_input_density.calcRMS()
49 IMP.em.write_map(sampled_input_density, "vv0.mrc", IMP.em.MRCReaderWriter())
50 # 3.2 calculate the cross-correlation score, which should be close to 1
52  dmap, sampled_input_density, sampled_input_density.get_header().dmin)
53 print "The CC score of the native transformation is:", best_score
54 
55 # 4. To denostrate local fitting we locally rotate and translate the
56 # protein and show how we can go back to the correct placement.
57 
58 # 4.1 define a local transformatione
60  IMP.algebra.get_unit_bounding_box_3d())
61 axis = IMP.algebra.get_random_vector_on(IMP.algebra.get_unit_sphere_3d())
62 rand_angle = random.uniform(-70. / 180 * math.pi, 70. / 180 * math.pi)
63 r = IMP.algebra.get_rotation_about_axis(axis, rand_angle)
64 local_trans = IMP.algebra.Transformation3D(r, translation)
65 # 4.2 rotate the protein
66 # prot_xyz=IMP.core.XYZs(IMP.core.get_leaves(mh))
67 # for xyz in prot_xyz:
68 # xyz.set_coordinates(local_trans.get_transformed(xyz.get_coordinates()))
69 # 4.2 set the protein as a rigid body
71 prot_rb = IMP.core.RigidMember(IMP.core.get_leaves(mh)[0]).get_rigid_body()
72 # 4.3 apply the trasnformation to the protein
73 IMP.core.transform(prot_rb, local_trans)
74 m.evaluate(False) # to make sure the transformation was applied
75 # 4.4 print the new correlation score, should be lower than before
76 print len(IMP.core.get_leaves(mh))
77 IMP.atom.write_pdb(mh, "input2.pdb")
78 sampled_input_density.resample()
79 sampled_input_density.calcRMS()
80 IMP.em.write_map(sampled_input_density, "vv.mrc", IMP.em.MRCReaderWriter())
82  dmap, sampled_input_density, sampled_input_density.get_header().dmin)
83 start_rmsd = IMP.atom.get_rmsd(IMP.core.XYZs(ps), IMP.core.XYZs(ps_ref))
84 print "The start score is:", start_score, "with rmsd of:", start_rmsd
85 # 5. apply local fitting
86 # 5.1 run local fitting
87 print "preforming local refinement, may run for 3-4 minutes"
88 # translate the molecule to the center of the density
90  IMP.algebra.get_identity_rotation_3d(), dmap.get_centroid() - IMP.core.get_centroid(ps)))
91 m.evaluate(False) # to make sure the transformation was applied
92 sampled_input_density.resample()
93 sampled_input_density.calcRMS()
96  dmap, sampled_input_density, sampled_input_density.get_header().dmin)
97 print "The score after centering is:", score2, "with rmsd of:", rmsd
98 # IMP.em.local_rigid_fitting_grid_search(
99 # ps,IMP.core.XYZR.get_radius_key(),
100 # IMP.atom.Mass.get_mass_key(),
101 # dmap,fitting_sols)
102 
104 fitting_sols = IMP.em.local_rigid_fitting(
105  mh.get_particle(), refiner,
107  dmap, [], 2, 10, 10)
108 
109 # 5.2 report best result
110 # 5.2.1 transform the protein to the preferred transformation
111 print "The start score is:", start_score, "with rmsd of:", start_rmsd
112 for i in range(fitting_sols.get_number_of_solutions()):
113  IMP.core.transform(prot_rb, fitting_sols.get_transformation(i))
114  # prot_rb.set_reference_frame(IMP.algebra.ReferenceFrame3D(fitting_sols.get_transformation(i)))
115  m.evaluate(False) # to make sure the transformation was applied
116 # 5.2.2 calc rmsd to native configuration
117  rmsd = IMP.atom.get_rmsd(
119  IMP.atom.write_pdb(mh, "temp_" + str(i) + ".pdb")
120  print "Fit with index:", i, " with cc: ", 1. - fitting_sols.get_score(i), " and rmsd to native of:", rmsd
121  IMP.atom.write_pdb(mh, "sol_" + str(i) + ".pdb")
123  prot_rb, fitting_sols.get_transformation(i).get_inverse())
124 print "done"
Simple 3D transformation class.
void write_pdb(const Selection &mhd, base::TextOutput out, unsigned int model=1)
static FloatKey get_mass_key()
Return the hierarchy leaves under a particle.
Definition: LeavesRefiner.h:25
void set_log_level(LogLevel l)
Set the current global log level.
void add_radii(Hierarchy d, const ForceFieldParameters *ffp=get_all_atom_CHARMM_parameters(), FloatKey radius_key=FloatKey("radius"))
void set_check_level(CheckLevel tf)
Control runtime checks in the code.
algebra::Vector3D get_centroid(const XYZs &ps)
Get the centroid.
double get_rmsd(const core::XYZs &s0, const core::XYZs &s1, const IMP::algebra::Transformation3D &tr_for_second)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
void write_map(DensityMap *m, std::string filename)
Rotation3D get_rotation_about_axis(const Vector3D &axis, double angle)
Generate a Rotation3D object from a rotation around an axis.
Class for sampling a density map from particles.
static double cross_correlation_coefficient(const DensityMap *grid1, const DensityMap *grid2, float grid2_voxel_data_threshold, bool allow_padding=false, FloatPair norm_factors=FloatPair(0., 0.))
Calculates the cross correlation coefficient between two maps.
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
See IMP.em for more information.
Definition: CoarseCC.h:23
DensityMap * read_map(std::string filename)
static const IMP::core::HierarchyTraits & get_traits()
Get the molecular hierarchy HierararchyTraits.
std::string get_example_path(std::string file_name)
Return the path to installed example data for this module.
See IMP.core for more information.
IMP::core::RigidBody create_rigid_body(Hierarchy h)
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:228
VectorD< D > get_random_vector_on(const SphereD< D > &s)
Generate a random vector on a sphere with uniform density.
Select all CA ATOM records.
Definition: pdb.h:76
See IMP.atom for more information.
void read_pdb(base::TextInput input, int model, Hierarchy h)
FittingSolutions local_rigid_fitting(kernel::Particle *p, Refiner *refiner, const FloatKey &weight_key, DensityMap *dmap, OptimizerStates display_log, Int number_of_optimization_runs=5, Int number_of_mc_steps=10, Int number_of_cg_steps=100, Float max_translation=2., Float max_rotation=.3, bool fast=true)
Local rigid fitting of a rigid body.
Class for storing model, its restraints, constraints, and particles.
Definition: kernel/Model.h:72