 |
Analyze the convergence of the IMP.em.FitRestraint. The script build a simple model and then displays the derivatives, em score and how well conjugate gradients converges under various displacements of the model.
import IMP.display
import IMP.em
use_rigid_bodies=True
bd= 10
radius=10
m= IMP.Model()
p= IMP.Particle(m)
IMP.atom.Mass.setup_particle(p, 10000)
d.set_radius(radius)
# Set up the particle as either a rigid body or a simple ball
if use_rigid_bodies:
prb= IMP.Particle(m)
prb.set_name("rigid body")
d.set_coordinates(IMP.algebra.Vector3D(0,0,0))
drb.add_member(p)
print "initial frame", drb.get_reference_frame()
fp= prb
drb.set_coordinates_are_optimized(True)
refiner= IMP.core.TableRefiner()
refiner.add_particle(prb, [p])
to_move= drb
print [p.get_name() for p in refiner.get_refined(prb)]
fp=d
else:
fp= d
to_move=d
d.set_coordinates_are_optimized(True)
refiner=None
bb= IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(-bd-radius, -bd-radius, -bd-radius),
IMP.algebra.Vector3D( bd+radius, bd+radius, bd+radius))
dheader = IMP.em.create_density_header(bb,1)
dheader.set_resolution(1)
dmap = IMP.em.SampledDensityMap(dheader)
dmap.set_particles([p])
dmap.resample()
# computes statistic stuff about the map and insert it in the header
dmap.calcRMS()
rs= IMP.RestraintSet()
m.add_restraint(rs)
#rs.set_weight(.003)
# if rigid bodies are used, we need to define a refiner as
# FitRestraint doesn't support just passing all the geometry
r= IMP.em.FitRestraint([fp], dmap)
rs.add_restraint(r)
g.set_name("deriv")
w= IMP.display.PymolWriter("derivatives.pym")
# kind of abusive
steps=4
m.set_log_level(IMP.SILENT)
opt= IMP.core.ConjugateGradients(m)
def try_point(i, j, k):
print "trying", i,j,k
vc=IMP.algebra.Vector3D(i,j,k)
to_move.set_coordinates(vc)
# display the score at this position
cg= IMP.display.SphereGeometry(IMP.algebra.Sphere3D(vc, 1))
cg.set_name("score")
v=m.evaluate(True)
cg.set_color(IMP.display.get_hot_color(v))
w.add_geometry(cg)
print "score and derivatives", v, to_move.get_derivatives()
w.add_geometry(g)
opt.optimize(10)
print "after", d.get_coordinates()
mag= to_move.get_coordinates().get_magnitude()
converge_color= IMP.display.get_grey_color(1.0/(1.0+mag))
# display the distance after optimization at this position
sg= IMP.display.SphereGeometry(IMP.algebra.Sphere3D(vc, 1))
sg.set_color(converge_color)
sg.set_name("converge")
w.add_geometry(sg)
try_point(-bd,-bd,-bd)
# For a more informative (but much slower) test, use the following instead:
#for i in range(-bd, bd+1, 2*bd/steps):
# for j in range(-bd, bd+1, 2*bd/steps):
# for k in range(-bd, bd+1, 2*bd/steps):
# try_point(i, j, k)
|
|
The example creates a simple mrc file that is filled uniformly with 1s. In addition, a chimera file is written with a marker at each corner of the density. When both files are opened, the density should be centered among the markers. This can be used for testing for registration errors when reading and writing density maps.
import IMP.em
import IMP.display
bb= IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(1,1,1),
IMP.algebra.Vector3D(11,11,11))
dm= IMP.em.create_density_map(bb, 2.5)
dmbb= IMP.em.get_bounding_box(dm)
print dmbb, bb
for i in range(0, dm.get_number_of_voxels()):
dm.set_value(i,1)
nm=IMP.create_temporary_file_name("cube", ".mrc")
print nm
IMP.em.write_map(dm, nm)
nm=IMP.create_temporary_file_name("cube", ".py")
print nm
w= IMP.display.ChimeraWriter(nm)
g= IMP.display.SphereGeometry(IMP.algebra.Sphere3D(v, 1))
w.add_geometry(g)
|
A simple example showing how to set up a fit restraint. The number of spheres and resolution are randomly chosen and so should not be considered significant.
import IMP.em
import IMP.core
import IMP.atom
IMP.set_log_level(IMP.SILENT)
m= IMP.Model()
#1. setup the input protein
##1.1 select a selector.
##1.2 read the protein
mh=IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"),m,sel)
ps=IMP.core.get_leaves(mh)
#2. read the density map
resolution=8.
voxel_size=1.5
dmap.get_header_writable().set_resolution(resolution)
#3. calculate the cross correlation between the density and the map
#4. add a fitting restraint
r= IMP.em.FitRestraint(ps, dmap)
m.add_restraint(r)
print "The fit of the particles in the density is:",r.evaluate(False)
|
Shows how to generate a density map of fixed dimension and how to sample particles within this density map.
import IMP
import IMP.atom
import IMP.em
m= IMP.Model()
#1. read a protein and get its bounding box dimension
# read protein
mh=IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"),m,sel)
apix=1.
resolution=6.
# compute bbox, and map size in voxels
bbox = IMP.atom.get_bounding_box(mh)
#diag = bbox.get_corner(0) - bbox.get_corner(1)
#nx = int(bboxCoverage * diag[0] / apix)
#ny = int(bboxCoverage * diag[1] / apix)
#nz = int(bboxCoverage * diag[2] / apix)
#create a density header of the requested size
dheader = IMP.em.create_density_header(bbox,apix)
dheader.set_resolution(resolution)
dmap = IMP.em.SampledDensityMap(dheader)
ps = IMP.core.get_leaves(mh)
dmap.set_particles(ps)
#dmap.get_header_writable().set_number_of_voxels(nx,ny,nz)
dmap.resample()
dmap.calcRMS() # computes statistic stuff about the map and insert it in the header
|
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.
import IMP.em
import IMP.core
import IMP.atom
import random,math
IMP.set_log_level(IMP.SILENT)
IMP.set_check_level(IMP.NONE)
m= IMP.Model()
#1. setup the input protein
##1.1 select a selector.
#using NonWater selector is more accurate but slower
#sel=IMP.atom.NonWaterPDBSelector()
##1.2 read the protein
mh=IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"),m,sel)
mh_ref=IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"),m,sel)
##1.3 add radius info to each atom, otherwise the resampling would fail.
IMP.atom.add_radii(mh_ref)
ps= IMP.core.get_leaves(mh)
ps_ref= IMP.core.get_leaves(mh_ref)
#2. read the density map of the protein
resolution=8.
voxel_size=1.5
dmap.get_header_writable().set_resolution(resolution)
#3. The protein is now fitted correctly in the density. We can validate
#that by making sure that the cross-correlation score is close to 1.
##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.
sampled_input_density = IMP.em.SampledDensityMap(dmap.get_header())
sampled_input_density.set_particles(ps)
sampled_input_density.resample()
sampled_input_density.calcRMS()
#3.2 calculate the cross-correlation score, which should be close to 1
best_score=IMP.em.CoarseCC.cross_correlation_coefficient(
dmap,sampled_input_density,sampled_input_density.get_header().dmin)
print "The CC score of the native transformation is:",best_score
#4. To denostrate local fitting we locally rotate and translate the protein and show how we can go back to the correct placement.
##4.1 define a local transformatione
translation = IMP.algebra.get_random_vector_in(IMP.algebra.get_unit_bounding_box_3d())
axis = IMP.algebra.get_random_vector_on(IMP.algebra.get_unit_sphere_3d())
rand_angle = random.uniform(-70./180*math.pi,70./180*math.pi)
r= IMP.algebra.get_rotation_about_axis(axis, rand_angle);
local_trans=IMP.algebra.Transformation3D(r,translation)
##4.2 rotate the protein
# prot_xyz=IMP.core.XYZs(IMP.core.get_leaves(mh))
# for xyz in prot_xyz:
# xyz.set_coordinates(local_trans.get_transformed(xyz.get_coordinates()))
##4.2 set the protein as a rigid body
prot_rb=IMP.core.RigidMember(IMP.core.get_leaves(mh)[0]).get_rigid_body()
##4.3 apply the trasnformation to the protein
IMP.core.transform(prot_rb,local_trans)
m.evaluate(None)#to make sure the transformation was applied
##4.4 print the new correlation score, should be lower than before
print len(IMP.core.get_leaves(mh))
IMP.atom.write_pdb(mh,"input2.pdb")
sampled_input_density.resample()
sampled_input_density.calcRMS()
start_score=IMP.em.CoarseCC.cross_correlation_coefficient(
dmap,sampled_input_density,sampled_input_density.get_header().dmin)
start_rmsd=IMP.atom.get_rmsd(IMP.core.XYZs(ps),IMP.core.XYZs(ps_ref))
print "The start score is:",start_score, "with rmsd of:",start_rmsd
##5. apply local fitting
## 5.1 run local fitting
print "preforming local refinement, may run for 3-4 minutes"
## translate the molecule to the center of the density
IMP.core.transform(prot_rb,IMP.algebra.Transformation3D(IMP.algebra.get_identity_rotation_3d(),dmap.get_centroid()-IMP.core.get_centroid(ps)))
m.evaluate(None)#to make sure the transformation was applied
sampled_input_density.resample()
sampled_input_density.calcRMS()
rmsd=IMP.atom.get_rmsd(IMP.core.XYZs(ps),IMP.core.XYZs(ps_ref))
dmap,sampled_input_density,sampled_input_density.get_header().dmin)
print "The score after centering is:",score2, "with rmsd of:",rmsd
# IMP.em.local_rigid_fitting_grid_search(
# ps,IMP.core.XYZR.get_radius_key(),
# IMP.atom.Mass.get_mass_key(),
# dmap,fitting_sols)
refiner = IMP.core.LeavesRefiner(IMP.atom.Hierarchy.get_traits())
fitting_sols=IMP.em.local_rigid_fitting(
mh.get_particle(),refiner,
dmap,[],2,10,10)
## 5.2 report best result
### 5.2.1 transform the protein to the preferred transformation
print "The start score is:",start_score, "with rmsd of:",start_rmsd
for i in range(fitting_sols.get_number_of_solutions()):
IMP.core.transform(prot_rb,fitting_sols.get_transformation(i))
#prot_rb.set_reference_frame(IMP.algebra.ReferenceFrame3D(fitting_sols.get_transformation(i)))
m.evaluate(None)#to make sure the transformation was applied
## 5.2.2 calc rmsd to native configuration
rmsd=IMP.atom.get_rmsd(IMP.core.XYZs(ps),IMP.core.XYZs(IMP.core.get_leaves(mh_ref)))
print "Fit with index:",i," with cc: ",1.-fitting_sols.get_score(i), " and rmsd to native of:",rmsd
IMP.core.transform(prot_rb,fitting_sols.get_transformation(i).get_inverse())
print "done"
|
|
A simple example showing how to simulate density from a protein. IMP uses a Gaussian smoothing kernel. see SampledDensityMap::resample for documentation.
import IMP.em
import IMP.core
import IMP.atom
m= IMP.Model()
#read protein
mh=IMP.atom.read_pdb(IMP.em.get_example_path("input.pdb"),m,sel)
#add radius info to each atom, otherwise the resampling would fail.
ps= IMP.core.get_leaves(mh)
#decide on resolution and spacing you would like to simulate to
resolution=10.
apix=1.5
dmap=IMP.em.particles2density(ps,resolution,apix)
#write out the map in the favorite format (xplor, mrc, em and spider are supported)
|