Hi list, In order to have a better insight I modified my script, keeping only the em.FitRestraint in my model (script attached). The score for the initial model (5.96046447754e-08) is reassuringly near 0; which is reasonable when comparing two maps sampled from the same structure. Nevertheless, when chains are randomly moved, the score does not radically improve, even when the transformations move them far away from their initial position. Moreover, I have an example where the score is minimized from 0.37 for the initial (moved chains) position to 0.49 for the optimized solution !? I may have missed something, and I would need some technical support to fix the problem quickly. --Ben |
#!/usr/bin/python
'''
Created on 22 juil. 2010
@author: schwarz
A first serious attempt to use IMP to model a complex
- a complex is forst loaded (1RDT : 2 domains and 2 peptides)
- a map of the complex is sampled at 10A with apix 2.5 A
- the chains are extracted and rigidified
- the chains are randomely moved
- two restraits are added to the model :
- fit sampled map with all atomic particles
- rigid bodies don't clash
- optimize m
'''
import IMP
import IMP.atom
import IMP.em
import IMP.core
import IMP.container
import IMP.algebra
import math,random
import datetime
#import time
import os
#pdbFilePrefix = "1E1Q"
pdbFilePrefix = "1RDT"
inputDataDirPath = "/Users/schwarz/Dev/IMPstuff/src/tests/fit_with_restraints_test/input"
outputDataDirPath = "/Users/schwarz/Dev/IMPstuff/src/tests/fit_with_restraints_test/output"
pdbFilePath = os.path.join(inputDataDirPath,pdbFilePrefix+".pdb")
params = {
"optimizer" : { "maxSteps": 100},
"densityMap" : { "resolution": 10, "apix": 2.5},
"MC" : { "maxTranslation": 5, "maxRotation": 10}
}
# Generation of a unique fileName prefix
d=datetime.datetime.now()
genericFilePrefix = "FitRestraint--" + "-".join([str(d.year),str(d.month),str(d.day),str(d.hour),str(d.minute)]) + "--"
print "*** START experiments "+ genericFilePrefix +" on structure "+ pdbFilePrefix
print " parameters are :"
for k,v in params.iteritems():
print k+":",v
IMP.set_log_level(IMP.SILENT)
m = IMP.Model()
print "** read pdb file twice (work and reference copies)"+pdbFilePath
# load and store a structure to work on
sel = IMP.atom.ATOMPDBSelector()
mh = IMP.atom.read_pdb(pdbFilePath,m,sel)
IMP.atom.add_radii(mh)
# also load and store a reference structure
mhref = IMP.atom.read_pdb(pdbFilePath,m,sel)
print "** Initial complex map"
print " - compute map"
resolution,voxelSize = params["densityMap"]["resolution"],params["densityMap"]["apix"]
print " - map for config resolution:"+str(resolution)+" apix:"+str(voxelSize)
dmap = IMP.em.SampledDensityMap(mh.get_leaves(),resolution,voxelSize)
mapName = genericFilePrefix + pdbFilePrefix+"-r"+str(resolution)+"-a"+str(voxelSize)+".mrc"
mapFilePath = os.path.join(outputDataDirPath,mapName)
print " - writing map to "+mapFilePath
IMP.em.write_map(dmap,mapFilePath,IMP.em.MRCReaderWriter())
print "** extract and rigidify protein chains from structure"
chains = {} # Will store chains on which I work
rigidChains = {} # Will store rigidBody for chains on which I work
refChains = {} # Will store reference chains
for chain in IMP.atom.Chains(IMP.atom.get_by_type(mh, IMP.atom.CHAIN_TYPE)) :
rigidChains[chain.get_id()] = IMP.atom.setup_as_rigid_body(chain)
chains[chain.get_id()] = chain
for chain in IMP.atom.Chains(IMP.atom.get_by_type(mhref, IMP.atom.CHAIN_TYPE)) :
refChains[chain.get_id()] = chain
print " (",str(len(refChains))," chains)"
print "** add restraints"
print " - FitRestraint between map and all particles"
m.add_restraint( IMP.em.FitRestraint(IMP.Particles(IMP.core.get_leaves(mh)),dmap) )
#print " - Volume exclusion restraint"
#print " - build singleton container"
#lsc = IMP.container.ListSingletonContainer()
#for rb in rigidChains.values() :
# lsc.add_particle(rb.get_particle())
#print " - add restraint"
##m.add_restraint( IMP.core.ExcludedVolumeRestraint(lsc, IMP.core.LeavesRefiner(rigidChains.values()[0].get_traits())) )
#m.add_restraint( IMP.core.ExcludedVolumeRestraint(lsc, IMP.core.LeavesRefiner(refChains.values()[0].get_traits())) )
print " - initial model evaluation : "+str(m.evaluate(False))
print "** initial structure"
print " - move chains randomly"
for chain in rigidChains.values():
translation = IMP.algebra.get_random_vector_in(IMP.algebra.get_unit_bounding_box_3d())
translation = translation * random.uniform( 1.0 , 5.0 )
axis = IMP.algebra.get_random_vector_on(IMP.algebra.get_unit_sphere_3d())
rand_angle = random.uniform(-50./180*math.pi,50./180*math.pi)
rotation = IMP.algebra.get_rotation_in_radians_about_axis(axis, rand_angle);
transformation = IMP.algebra.Transformation3D(rotation,translation)
IMP.core.transform( chain , transformation )
print " - save noisy structure"
noisyModelFilePath = os.path.join(outputDataDirPath,genericFilePrefix + pdbFilePrefix+"-noisy.pdb")
IMP.atom.write_pdb(mh,noisyModelFilePath)
print " - noisy model evaluation : "+str(m.evaluate(False))
print "** Optimization"
print " - create an optimizer"
#opti = IMP.core.SteepestDescent(m)
#opti = IMP.core.ConjugateGradients(m)
opti = IMP.core.MonteCarlo(m)
#opti.set_log_level(IMP.TERSE)
opti.set_log_level(IMP.WARNING)
print " - add movers to MC optimizer for all distinct chains"
for chain in rigidChains.values():
rb_mover = IMP.core.RigidBodyMover(chain, params["MC"]["maxTranslation"], params["MC"]["maxRotation"] )
opti.add_mover(rb_mover)
print " - optimize"
opti.optimize(params["optimizer"]["maxSteps"])
print "** Analysis"
print " - optimized model evaluation : "+str(m.evaluate(False))
print " - save optimized model"
optimizedModelFilePath = os.path.join(outputDataDirPath,genericFilePrefix + pdbFilePrefix+"-optimized.pdb")
IMP.atom.write_pdb(mh,optimizedModelFilePath)
print " - measure RMSD discrepency between native and optimized chains"
for chainID in refChains.keys() :
rmsd=IMP.atom.get_rmsd(
IMP.core.XYZs( chains[chainID].get_leaves() ) ,
IMP.core.XYZs( refChains[chainID].get_leaves() )
)
print " "+chainID+" : "+str(rmsd)
print " - save optimized complex map"
dmapOptimized = IMP.em.SampledDensityMap(mh.get_leaves(),resolution,voxelSize)
mapName = genericFilePrefix + pdbFilePrefix+"-optimized-r"+str(resolution)+"-a"+str(voxelSize)+".mrc"
mapFilePath = os.path.join(outputDataDirPath,mapName)
IMP.em.write_map(dmapOptimized,mapFilePath,IMP.em.MRCReaderWriter())
Dr. Benjamin SCHWARZ Biocomputing group Email : "> Voice : +33 (0)3 68 85 47 30 FAX : +33 (0)3 68 85 47 18 |
Structural Biology & Genomics Dept. - IGBMC 1 rue Laurent Fries
BP 10142 F - 67404 Illkirch CEDEX |