Following Daniel's suggestion I added a set of Movers to my MC optimizer (see code attached). Now the optimizer actually does something, though the results are very strange. It is much probable I poorly chose my parameters, and 1000 steps are not so much for a MC optimization. It also appears I have a problem with the cost function : it evaluates to 18.36 for the native configuration and to 1371 for an initial "noisy" configuration not so far from the initial configuration. It then evaluates to 1.03 for an "optimized" configuration where all four chains are placed far outside the sampled density. I did not discover yet if and how it is possible to set parameters for the global cost function, maybe my problem has something to do with it. Anyhow, I would be very grateful if someone would glance at my script, just to check I did use IMP in a proper manner. --Ben |
'''
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 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")
config = (10,2.5)
maxStep = 500
MCmaxTranslation, MCmaxRotation = (5,10)
print "*** START"
print " working on " + pdbFilePrefix
IMP.set_log_level(IMP.SILENT)
m = IMP.Model()
print "** read pdb file twice (work and reference copies)"+pdbFilePath
sel = IMP.atom.ATOMPDBSelector()
mh = IMP.atom.read_pdb(pdbFilePath,m,sel)
IMP.atom.add_radii(mh)
mhref = IMP.atom.read_pdb(pdbFilePath,m,sel)
refChains = {}
for chain in IMP.atom.Chains(IMP.atom.get_by_type(mhref, IMP.atom.CHAIN_TYPE)):
refChains[chain.get_id()]=chain
print "** create complex map"
resolution,voxelSize = config
print " - map for config resolution:"+str(resolution)+" apix:"+str(voxelSize)
dmap = IMP.em.SampledDensityMap(mh.get_leaves(),resolution,voxelSize)
mapName = 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 = IMP.atom.Chains(IMP.atom.get_by_type(mh, IMP.atom.CHAIN_TYPE))
rigidChains = {}
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
print " (",str(len(chains))," 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(chains.values()[0].get_traits())) )
print " - initial model evaluation : "+str(m.evaluate(False))
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,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)
print " - add movers to optimizer for all distinct chains"
for chain in rigidChains.values():
rb_mover = IMP.core.RigidBodyMover(chain, MCmaxTranslation, MCmaxRotation )
opti.add_mover(rb_mover)
print " - optimize"
opti.optimize(maxStep)
print "** Analysis"
print " - optimized model evaluation : "+str(m.evaluate(False))
print " - save optimized model"
optimizedModelFilePath = os.path.join(outputDataDirPath,pdbFilePrefix+"-optimized.pdb")
IMP.atom.write_pdb(mh,optimizedModelFilePath)
print "-- measure discrepency between native and optimized chains"
for chainID in chains.keys() :
rmsd=IMP.atom.get_rmsd(
IMP.core.XYZs( chains[chainID].get_leaves() ) ,
IMP.core.XYZs( refChains[chainID].get_leaves() )
)
print " "+chainID+" : "+str(rmsd)
resolution,voxelSize = config
dmap = IMP.em.SampledDensityMap(mh.get_leaves(),resolution,voxelSize)
mapName = pdbFilePrefix+"-optimized-r"+str(resolution)+"-a"+str(voxelSize)+".mrc"
mapFilePath = os.path.join(outputDataDirPath,mapName)
print " -- save optimized complex map for config resolution:"+str(resolution)+" apix:"+str(voxelSize)+" to "+mapFilePath
IMP.em.write_map(dmap,mapFilePath,IMP.em.MRCReaderWriter())Le 23 juil. 2010 à 19:44, Daniel Russel a écrit :
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 |