Re: [IMP-users] problems with a sample script (basic IMP optimization)

[Benjamin SCHWARZ <>]

…And it is even better with the script.

'''
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 = 100

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 " - 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 "** 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 "** Optimization"
print " - create an optimizer"
#opti = IMP.core.SteepestDescent(m)
#opti = IMP.core.ConjugateGradients(m)
opti = IMP.core.MonteCarlo(m)
print " - optimize"
opti.optimize(maxStep)
print " - optimized model evaluation : "+str(m.evaluate(False))



print "** Analysis"
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)

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