 |
This illustrates a basic main loop for optimization and searching for the best scoring conformation.
import IMP.example
import IMP.statistics
(m,c)=IMP.example.create_model_and_particles()
r= IMP.container.PairsRestraint(ps, IMP.container.ClosePairContainer(c, 2.0))
m.add_restraint(r)
# we don't want to see lots of log messages about restraint evaluation
m.set_log_level(IMP.WARNING)
# the container (c) stores a list of particles, which are alse XYZR particles
# we can construct a list of all the decorated particles
xyzrs= c.get_particles()
s= IMP.core.MCCGSampler(m)
s.set_number_of_attempts(10)
# but we do want something to watch
s.set_log_level(IMP.TERSE)
s.set_number_of_monte_carlo_steps(10)
# find some configurations which move the particles far apart
configs= s.get_sample();
for i in range(0, configs.get_number_of_configurations()):
configs.load_configuration(i)
# print out the sphere containing the point set
# - Why? - Why not?
sphere= IMP.core.get_enclosing_sphere(xyzrs)
print sphere
# cluster the solutions based on their coordinates
e= IMP.statistics.ConfigurationSetXYZEmbedding(configs, c)
# of course, this doesn't return anything of interest since the points are
# randomly distributed, but, again, why not?
clustering = IMP.statistics.create_lloyds_kmeans(e, 3, 1000)
for i in range(0,clustering.get_number_of_clusters()):
# load the configuration for a central point
configs.load_configuration(clustering.get_cluster_representative(i))
sphere= IMP.core.get_enclosing_sphere(xyzrs)
print sphere
|
This example shows how to set up an optimization involving several particles constrained to be connected in a loop. It uses non bonded lists and a variety of restraints.
import IMP
import IMP.core
import random
import IMP.display
# A trivial example that constructs a set of particles which are restrained
# to form a chain via bonds between successive particles. In addition
# the head and the tail of the chain are restrained to be close to one
# another.
IMP.set_log_level(IMP.TERSE)
m= IMP.Model()
m.set_log_level(IMP.SILENT)
# The particles in the chain
ps=IMP.core.create_xyzr_particles(m, 10, 1.0)
chain= IMP.container.ListSingletonContainer(ps, "chain")
# create a spring between successive particles
bonds= IMP.container.ConsecutivePairContainer(ps)
hdps=IMP.core.HarmonicDistancePairScore(2,1)
chainr= IMP.container.PairsRestraint(hdps,bonds)
chainr.set_name("The chain restraint")
m.add_restraint(chainr)
# If you want to inspect the particles
# Notice that each bond is a particle
for p in m.get_particles():
p.show()
# Prevent non-bonded particles from penetrating one another
nbl= IMP.container.ClosePairContainer(chain, 0,2)
bpc=IMP.container.InContainerPairFilter(bonds) # exclude existing bonds
nbl.add_pair_filter(bpc)
"excluded volume")
m.add_restraint(lr)
# Tie the ends of the chain
(ps[0], ps[-1]))
tie.set_name("tie ends")
m.add_restraint(tie)
s= IMP.core.MCCGSampler(m) # sample using MC and CG
s.set_number_of_attempts(10)
m.set_maximum_score(1)
confs= s.get_sample()
print "Found", confs.get_number_of_configurations(), "configurations"
for i in range(0, confs.get_number_of_configurations()):
confs.load_configuration(i)
for p in chain.get_particles():
d.add_geometry(IMP.core.XYZRGeometry(p))
# print out info about used modules so that the versions are known
#IMP.show_used_modules()
|
|
When trying to understand what is going on in IMP, it can often be useful to view the dependency graph, that is, the graph showing how various entities relate to one another. In it, an arrow leads from an IMP::Container or IMP::Particle to an IMP::Restraint if the IMP::Restraint reads from that container or particle. Similarly, an arrow leads from an IMP::Container or IMP::Particle to an IMP::ScoreState if the score state reads from it, and an arrow leads from an IMP::ScoreState to an IMP::Container or IMP::Particle if the score state updates the particle.
import IMP
import IMP.atom
import IMP.container
def create_representation():
m= IMP.Model()
all.set_name("the universe")
def create_protein(name, ds):
h=IMP.atom.create_protein(m, name, 10, ds)
leaves= IMP.atom.get_leaves(h)
all.add_child(h)
for c in h.get_children()],
1)
if r:
m.add_restraint(r)
def create_protein_from_pdbs(name, files):
def create_from_pdb(file):
sls=IMP.base.SetLogState(IMP.NONE)
del sls
#IMP.atom.show_molecular_hierarchy(t)
c=IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
if c.get_number_of_children()==0:
IMP.atom.show_molecular_hierarchy(t)
# there is no reason to use all atoms, just approximate the pdb shape instead
10.0/2.0)
# make the simplified structure rigid
rb.set_coordinates_are_optimized(True)
return s
if len(files) >1:
p= IMP.Particle(m)
h.set_name(name)
for i, f in enumerate(files):
c=create_from_pdb(f)
h.add_child(c)
c.set_name(name+" chain "+str(i))
for c in h.get_children()],
1)
if r:
m.add_restraint(r)
else:
h= create_from_pdb(files[0])
h.set_name(name)
all.add_child(h)
create_protein("Nup85", 570)
create_protein("Nup84", 460)
create_protein("Nup145C", 442)
create_protein("Nup120", [0, 500, 761])
create_protein("Nup133", [0, 450, 778, 1160])
create_protein_from_pdbs("Seh1", ["seh1.pdb"])
create_protein_from_pdbs("Sec13", ["sec13.pdb"])
return (m, all)
def create_restraints(m, all):
def add_connectivity_restraint(s):
m.add_restraint(r)
def add_distance_restraint(s0, s1):
r=IMP.atom.create_distance_restraint(s0,s1, 0, 1)
m.add_restraint(r)
evr=IMP.atom.create_excluded_volume_restraint([all])
m.add_restraint(evr)
s0=S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)])
s1=S(hierarchy=all, molecule="Nup84")
s2=S(hierarchy=all, molecule="Sec13")
add_connectivity_restraint([s0,s1,s2])
add_distance_restraint(S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Nup85"))
add_distance_restraint(S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Nup120",
residue_indexes= [(500, 762)]))
add_distance_restraint(S(hierarchy=all, molecule="Nup84"),
S(hierarchy=all, molecule="Nup133",
residue_indexes=[(778, 1160)]))
add_distance_restraint(S(hierarchy=all, molecule="Nup85"),
S(hierarchy=all, molecule="Seh1"))
add_distance_restraint(S(hierarchy=all, molecule="Nup145C",
residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Sec13"))
# now do the actual work
(m,all)= create_representation()
create_restraints(m, all)
# we can get the full dependency graph for the whole model with all the restraints
# but it is pretty complex
dg= IMP.get_dependency_graph(m)
IMP.show_graphviz(dg);
# better thing to do is to get the "pruned" graph
pdg= IMP.get_pruned_dependency_graph(m)
try:
# these all open new windows which must be closed to continue
# also, the graph is no where near as nice as displayed by
# IMP.show_graphviz below
#import matplotlib
# the engine to be used must be selected before pyplot is imported
#matplotlib.use("macosx")
#import matplotlib.pyplot as plt
# the method below requires the altgraph python package
#xg=IMP.get_networkx_graph(pdg)
#import networkx
#networkx.draw_spectral(xg)
#plt.show()
pass
except:
try:
IMP.show_altgraph(pdg)
except:
print 'Need networkx and matplotlib or altgraph to display graphs interactively'
IMP.show_graphviz(pdg)
|
|
This example docks several proteins using excluded volume and crosslinking terms. To set the pdbs and crosslinks to use, edit the data at the start of the python script.
import IMP
import IMP.core
import IMP.algebra
import IMP.atom
import IMP.container
import random
import sys
try:
import networkx as nx
import networkx.algorithms
except:
print "Script requires networkx to run"
sys.exit()
# remove internal checks
IMP.set_check_level(IMP.USAGE)
pdbs=[IMP.get_example_path('dock_data/chainf.pdb'),
IMP.get_example_path('dock_data/chaind.pdb')]
xlinks=[
(256, "D", 89, "F", 17, 0.1),
(259, "D", 84, "F", 17, 0.1),
(259, "D", 89, "F", 17, 0.1),
(269, "D", 89, "F", 17, 0.1),
(128, "F", 275, "D", 17, 0.1),
(130, "F", 275, "D", 17, 0.1),
(144, "F", 269, "D", 17, 0.1),
(150, "F", 256, "D", 17, 0.1),
(150, "F", 259, "D", 17, 0.1),
(150, "F", 269, "D", 17, 0.1),
(150, "F", 179, "D", 17, 0.1),
(166, "F", 186, "D", 17, 0.1),
(170, "D", 89, "F", 17, 0.1),
(179, "D", 84, "F", 17, 0.1)]
def read_pdbs(list_pdb_file):
"""read pdbs from an external list file
create a simplified representation"""
chains=[]
chain_id={}
for pdb in pdbs:
prot = IMP.atom.read_pdb(pdb, m,
ch=IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
chains+=ch
print ch
for c in chains:
chain_id[c]=IMP.atom.Chain(c).get_id()
return chains,chain_id
def setup_radii(chains):
"setup the radii for each residue"
for c in chains:
for res in c.get_children():
dres=IMP.core.XYZR(res)
rt=IMP.atom.Residue(res).get_residue_type()
dres.set_radius(rg)
def add_excluded_volume(chains):
"""add excluded volume score on the coarse grained c-alpha model
residues are represented by beads with conveniently defined excluded volume radius"""
rs = IMP.RestraintSet('excluded_volume')
kappa_=0.1
IMP.atom.get_by_type
for c in chains:
lsc.add_particles(IMP.atom.get_leaves(c))
evr=IMP.core.ExcludedVolumeRestraint(lsc,kappa_)
rs.add_restraint(evr)
m.add_restraint(rs)
restraints['excluded_volume']=rs
def setup_move_them_all_MonteCarlo_internal(chains,mc_dx=0.3,mc_dang=0.1,mc_kt=1.0):
"""setup rigid body monte carlo move set
mc_dx is the traslation move
mc_dang is the rotation angle"""
mc=IMP.core.MonteCarlo(m)
mc.set_return_best(False)
mc.set_kt(10)
for c in chains:
p=IMP.atom.get_leaves(c)[0]
rb=IMP.core.RigidMember(p).get_rigid_body()
mv= IMP.core.RigidBodyMover(rb, mc_dx, mc_dang)
mc.add_mover(mv)
return mc
def setup_move_them_all_MonteCarlo_external(connected_chain_list,rb_ext_list=[],
mc_dx=0.3,mc_dang=0.1,mc_kt=1.0,return_best=False):
"""setup rigid body monte carlo move set
mc_dx is the traslation move
mc_dang is the rotation angle"""
print 'RIGID EXT BODY LIST', [r.get_name() for r in rb_ext_list]
for rb in rb_ext_list:
rb_ext_list=[]
mc=IMP.core.MonteCarlo(m)
mc.set_return_best(return_best)
mc.set_kt(mc_kt)
print connected_chain_list
for c in connected_chain_list:
rb_tmp_list=[]
for s in c:
p=IMP.atom.get_leaves(s)[0]
rba=IMP.core.RigidMember(p).get_rigid_body()
rb_tmp_list.append(rba)
print 'RIGID BODY LIST', [r.get_name( ) for r in rb_tmp_list]
p=IMP.Particle(m)
p.set_name("root rigid body "+chain_id[s])
rb=IMP.core.RigidBody.setup_particle(p,rb_tmp_list)
rb_ext_list.append(rb)
mv= IMP.core.RigidBodyMover(rb, mc_dx, mc_dang)
mc.add_mover(mv)
return mc,rb_ext_list
def setup_move_them_all_MonteCarlo_internal_2(chains,rb_ext_list=[],
mc_dx=0.3,mc_dang=0.1,mc_kt=1.0):
"""setup rigid body monte carlo move set
mc_dx is the traslation move
mc_dang is the rotation angle"""
print 'RIGID EXT BODY LIST', [r.get_name() for r in rb_ext_list]
for rb in rb_ext_list:
rb_ext_list=[]
mc=IMP.core.MonteCarlo(m)
mc.set_return_best(False)
mc.set_kt(10)
for c in chains:
p=IMP.atom.get_leaves(c)[0]
rb=IMP.core.RigidMember(p).get_rigid_body()
rb_ext_list.append(rb)
mv= IMP.core.RigidBodyMover(rb, mc_dx, mc_dang)
mc.add_mover(mv)
return mc,rb_ext_list
def shuffle_configuration(chains,bounding_box_length=200.0):
"shuffle configuration, used to restart the optimization"
hbbl=bounding_box_length/2
for c in chains:
p=IMP.atom.get_leaves(c)[0]
rb=IMP.core.RigidMember(p).get_rigid_body()
ub = IMP.algebra.Vector3D(-hbbl,-hbbl,-hbbl)
lb = IMP.algebra.Vector3D( hbbl, hbbl, hbbl)
bb = IMP.algebra.BoundingBox3D(ub, lb)
translation = IMP.algebra.get_random_vector_in(bb)
rotation = IMP.algebra.get_random_rotation_3d()
transformation = IMP.algebra.Transformation3D(rotation, translation)
rb.set_reference_frame(IMP.algebra.ReferenceFrame3D(transformation))
def add_restraints():
"""Handle the restraints defined above
sintax: resid1 chain1 resid2 chain2 max_distance strength_of_potential
example: 226 H 272 N 15.0 10.0"""
rs=IMP.RestraintSet('cross_links')
restraints={}
restraints_map={}
for tokens in xlinks:
r1=tokens[0]
c1=tokens[1]
r2=tokens[2]
c2=tokens[3]
distance=tokens[4]
strength=tokens[5]
try:
s1=IMP.atom.Selection(chains, residue_index=r1, chains=c1)
p1=s1.get_selected_particles()[0]
except:
print "residue %d of chain %s not found" % (r1,c1)
continue
try:
s2=IMP.atom.Selection(chains, residue_index=r2, chains=c2)
p2=s2.get_selected_particles()[0]
except:
print "residue %d of chain %s not found" % (r2,c2)
continue
d=distance
k=strength
hub= IMP.core.HarmonicUpperBound(d,k)
df= IMP.core.DistancePairScore(hub)
dr= IMP.core.PairRestraint(df, (p1, p2))
rs.add_restraint(dr)
a='cross_link-'+`r1`+'.'+c1+'-'+`r2`+'.'+c2+'.lower'
restraints_map
restraints[a]=dr
restraints_map[a]=(c1,c2)
m.add_restraint(rs)
restraints['cross_links']=rs
return restraints,restraints_map
def setup_md(temp=300.0, tau=0.01):
md.set_model(m)
md.assign_velocities(temp)
md.set_time_step(1.0)
st = IMP.atom.LangevinThermostatOptimizerState(md.get_simulation_particles(), temp,tau)
md.add_optimizer_state(st)
return md
def init_pdb(prefix='models'):
"append models to a pdb file"
pdbfile = prefix+'.pdb'
flpdb=open(pdbfile,'w')
flpdb.close()
def write_pdb(prefix='models'):
"append zero score models (solutions) to a pdb file"
pdbfile = prefix+'.pdb'
flpdb=open(pdbfile,'a')
tmpfile=".tmp.pdb"
IMP.atom.write_pdb(chains,tmpfile)
fltmp=open(tmpfile,'r')
pdb=fltmp.read()
flpdb.write(pdb)
flpdb.close()
fltmp.close()
def chain_pair_list(chains,nop_cutoff=0):
pair_list=[]
for c in chains:
sc1.add_particles(IMP.atom.get_leaves(c))
for s in chains:
if chain_id[s]>chain_id[c]:
sc2.add_particles(IMP.atom.get_leaves(s))
aaa= IMP.container.CloseBipartitePairContainer(sc1,sc2,0.0)
m.update()
nop=aaa.get_number_of_particle_pairs()
if nop > nop_cutoff :
pair_list.append((c,s))
else:
pair_list.append((c,c))
pair_list.append((s,s))
return pair_list
def chain_pair_list_based_on_restraint(chains,nop_cutoff=0):
pair_list=[]
id_chain=invert_dict(chain_id)
for i in restraints.keys():
(chaina_id,chainb_id)=restraints_map[i]
chaina=invert_dict(chaina_id)
chainb=invert_dict(chainb_id)
for c in chains:
sc1.add_particles(IMP.atom.get_leaves(c))
for s in chains:
if chain_id[s]>chain_id[c]:
sc2.add_particles(IMP.atom.get_leaves(s))
aaa= IMP.container.CloseBipartitePairContainer(sc1,sc2,0.0)
m.update()
nop=aaa.get_number_of_particle_pairs()
if nop > nop_cutoff :
pair_list.append((c,s))
else:
pair_list.append((c,c))
pair_list.append((s,s))
return pair_list
def create_graph(chains,nop_cutoff=0):
connected=[]
while len(connected)<=1:
G=nx.Graph()
pair_list=chain_pair_list(chains,nop_cutoff)
G.add_edges_from(pair_list)
connected=nx.algorithms.connected_components(G)
G.clear()
nop_cutoff+=10
return connected,nop_cutoff
def invert_dict(d):
inv = {}
for k, v in d.iteritems():
keys = inv.setdefault(v, [])
keys.append(k)
return inv
m = IMP.Model()
chain_id={}
rb_ext_list=[]
restraints={}
restraints_map={}
init_pdb("models")
init_pdb("solutions")
init_pdb("best")
chains,chain_id=read_pdbs("pdb.list")
write_pdb("best")
setup_radii(chains)
add_excluded_volume(chains)
restraints,restraints_map=add_restraints()
shuffle_configuration(chains,300.0)
delta_r=5.0
delta_a=0.2
delta_r_best=0.5
delta_a_best=0.1
temp=1.0
mc_int=setup_move_them_all_MonteCarlo_internal(chains,delta_r,delta_a,temp)
write_pdb("models")
nsteps_int=100
nsteps_ext=100
cg.set_model(m)
for steps in range(10000):
nop_cutoff=10000000
connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
mc_int,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r,delta_a,temp,return_best=False)
mc_int.optimize(nsteps_ext)
print ' INTERNAL', m.evaluate(False)
write_pdb("models")
nop_cutoff=0
connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
mc_ext,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r,delta_a,temp,return_best=False)
mc_ext.optimize(nsteps_ext)
print ' EXTERNAL', m.evaluate(False)
write_pdb("models")
nfs=mc_ext.get_number_of_forward_steps()
acceptance=float(nfs)/nsteps_ext
if acceptance<0.05:
delta_r=delta_r/2
if acceptance>0.3:
delta_r=delta_r*2
print 'ACCEPTANCE 1', acceptance, delta_r, delta_a
nop_cutoff=10000000
connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
mc_int,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r_best,delta_a_best,temp,return_best=True)
mc_int.optimize(nsteps_ext)
print ' INTERNAL', m.evaluate(False)
write_pdb("models")
nop_cutoff=0
connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
mc_ext,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r_best,delta_a_best,temp,return_best=True)
mc_ext.optimize(nsteps_ext)
print ' EXTERNAL', m.evaluate(False)
write_pdb("models")
nfs=mc_ext.get_number_of_forward_steps()
acceptance=float(nfs)/nsteps_ext
if acceptance<0.5:
delta_r_best=delta_r_best/2
if acceptance>0.8:
delta_r_best=delta_r_best*2
print 'ACCEPTANCE 2', acceptance, delta_r_best, delta_a_best
cg.optimize(10)
score=m.evaluate(False)
for i in restraints.keys():
print i, restraints[i].evaluate(False)
if score==0.0:
write_pdb("solutions")
shuffle_configuration(chains,300.0)
exit()
print 'ALL DONE'
|
|
A simple example showing how to use the graph interface for in python.
import IMP
m= IMP.Model()
# An undirected graph with an IMP::Object for each node
g= IMP.DependencyGraph()
vs=[]
ps=[]
for i in range(0,10):
ps.append(IMP.Particle(m))
vs.append(g.add_vertex(ps[-1]));
g.add_edge(vs[0], vs[1])
g.add_edge(vs[1], vs[2])
#try to use the altgraph package to visualize
try:
except:
print "Oh well, no altgraph"
try:
import matplotlib
# the engine to be used must be selected before pyplot is imported
matplotlib.use("macosx")
import matplotlib.pyplot as plt
# the method below requires the altgraph python package
xg=IMP.get_networkx_graph(g)
# the networkx visualization tools suck, so skip them
#import networkx
#networkx.draw(xg)
#networkx.draw_shell(xg)
#plt.show()
except:
print "networkx not fully installed"
g.remove_vertex(0)
# we can also try another show method
try:
except:
print "oh well, something not working with graphviz"
# finally, we can
# in and out neighbors are the same
for n in g.get_in_neighbors(8):
print g.get_vertex_name(n).get_name()
|
Use connectivity data taken from the Nup84 complex to provide coarse grained structure models of the complex. This example shows a (simple) full IMP modeling protocol from choice of representation and scoring to sampling and clustering and analysis of the solutions. The data and representation is similar to that of the Nup845 example in IMP.restrainer.
import IMP
import IMP.atom
import IMP.container
import IMP.display
import IMP.statistics
import IMP.example
# the spring constant to use, it doesn't really matter
k=100
# the target resolution for the representation, this is used to specify how detailed
# the representation used should be
resolution=300
# the box to perform everything in
bb=IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(-300,-300,-300),
IMP.algebra.Vector3D(300, 300, 300))
IMP.random_number_generator.seed(2025)
# this function creates the molecular hierarchies for the various involved proteins
def create_representation():
m= IMP.Model()
all.set_name("the universe")
# create a protein, represented as a set of connected balls of appropriate
# radii and number, chose by the resolution parameter and the number of
# amino acids.
def create_protein(name, ds):
h=IMP.atom.create_protein(m, name, resolution, ds)
# for convenience, have one molecular hierarchy containing all molecules
all.add_child(h)
for c in h.get_children()],
k)
if r:
m.add_restraint(r)
# only allow the particles to separate by 1 angstrom
m.set_maximum_score(r, k)
create_protein("Nup85", 570)
ct= IMP.atom.Selection(all, molecule="Nup85", terminus= IMP.atom.Selection.C)
d= IMP.core.XYZ(ct.get_selected_particles()[0])
# pin NUP85 to remove a degree of freedom
d.set_coordinates(IMP.algebra.Vector3D(1,0,0))
d.set_coordinates_are_optimized(False)
create_protein("Nup84", 460)
create_protein("Nup145C", 442)
create_protein("Nup120", [0, 500, 761])
create_protein("Nup133", [0, 450, 778, 1160])
create_protein("Seh1", 351)
create_protein("Sec13", 379)
return (m, all)
# create the needed restraints and add them to the model
def create_restraints(m, all):
def add_connectivity_restraint(s):
tr= IMP.core.TableRefiner()
rps=[]
for sc in s:
ps= sc.get_selected_particles()
rps.append(ps[0])
tr.add_particle(ps[0], ps)
# duplicate the IMP.atom.create_connectivity_restraint functionality
tr)
r= IMP.core.ConnectivityRestraint(score, IMP.container.ListSingletonContainer(rps))
m.add_restraint(r)
m.set_maximum_score(r, k)
def add_distance_restraint(s0, s1):
r=IMP.atom.create_distance_restraint(s0,s1, 0, k)
m.add_restraint(r)
# only allow the particles to separate by one angstrom
m.set_maximum_score(r, k)
evr=IMP.atom.create_excluded_volume_restraint([all])
m.add_restraint(evr)
# a Selection allows for natural specification of what the restraints act on
s0=S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)])
s1=S(hierarchy=all, molecule="Nup84", molecules=[])
s2=S(hierarchy=all, molecule="Sec13")
add_connectivity_restraint([s0,s1,s2])
add_distance_restraint(S(hierarchy=all, molecule="Nup145C",
residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Nup85"))
add_distance_restraint(S(hierarchy=all, molecule="Nup145C",
residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Nup120",
residue_indexes= [(500, 762)]))
add_distance_restraint(S(hierarchy=all, molecule="Nup84"),
S(hierarchy=all, molecule="Nup133",
residue_indexes=[(778, 1160)]))
add_distance_restraint(S(hierarchy=all, molecule="Nup85"),
S(hierarchy=all, molecule="Seh1"))
add_distance_restraint(S(hierarchy=all, molecule="Nup145C",
residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Sec13"))
#for l in IMP.atom.get_leaves(all):
# r= IMP.example.ExampleRestraint(l, k)
# m.add_restraint(r)
# m.set_maximum_score(k)
# find acceptable conformations of the model
def get_conformations(m):
sampler= IMP.core.MCCGSampler(m)
sampler.set_bounding_box(bb)
# magic numbers, experiment with them and make them large enough for things to work
sampler.set_number_of_conjugate_gradient_steps(100)
sampler.set_number_of_monte_carlo_steps(10)
sampler.set_number_of_attempts(20)
# We don't care to see the output from the sampler
sampler.set_log_level(IMP.SILENT)
# return the IMP.ConfigurationSet storing all the found configurations that
# meet the various restraint maximum scores.
cs= sampler.get_sample()
return cs
# cluster the conformations and write them to a file
def analyze_conformations(cs, all, gs):
# we want to cluster the configurations to make them easier to understand
# in the case, the clustering is pretty meaningless
embed= IMP.statistics.ConfigurationSetXYZEmbedding(cs,
cluster= IMP.statistics.create_lloyds_kmeans(embed, 10, 10000)
# dump each cluster center to a file so it can be viewed.
w= IMP.display.PymolWriter("cluster.pym")
for i in range(cluster.get_number_of_clusters()):
rep= cluster.get_cluster_representative(i)
cs.load_configuration(rep)
w.set_frame(i)
for g in gs:
w.add_geometry(g)
# now do the actual work
(m,all)= create_representation()
IMP.atom.show_molecular_hierarchy(all)
create_restraints(m, all)
# in order to display the results, we need something that maps the particles onto
# geometric objets. The IMP.display.Geometry objects do this mapping.
# IMP.core.XYZRGeometry map an IMP.core.XYZR particle onto a sphere
gs=[]
for i in range(all.get_number_of_children()):
color= IMP.display.get_display_color(i)
n= all.get_child(i)
name= n.get_name()
g.set_color(color)
gs.append(g)
cs= get_conformations(m)
print "found", cs.get_number_of_configurations(), "solutions"
# for each of the configuration, dump it to a file to view in pymol
w= IMP.display.PymolWriter("config.pym")
for i in range(0, cs.get_number_of_configurations()):
cs.load_configuration(i)
w.set_frame(i)
for g in gs:
w.add_geometry(g)
analyze_conformations(cs, all, gs)
|
|
Modify the Nup84 CG example by replacing a couple of the protein representations with higher resolution representations generated from PDB files. In addition, show how to visualize restraints and visualize the rejected conformations. Both are useful things to do when trying to figure out why optimization is not converging.
import IMP
import IMP.atom
import IMP.container
import IMP.display
import IMP.statistics
import IMP.example
# the spring constant to use, it doesn't really matter
k=100
# the target resolution for the representation
resolution=100
# the box to perform everything in
bb=IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(0,0,0),
IMP.algebra.Vector3D(300, 300, 300))
display_restraints=[]
def create_representation():
m= IMP.Model()
all.set_name("the universe")
def create_protein(name, ds):
h=IMP.atom.create_protein(m, name, resolution, ds)
leaves= IMP.atom.get_leaves(h)
all.add_child(h)
for c in h.get_children()],
k)
if r:
m.add_restraint(r)
display_restraints.append(r)
m.set_maximum_score(r, k)
def create_protein_from_pdbs(name, files):
def create_from_pdb(file):
sls=IMP.base.SetLogState(IMP.NONE)
del sls
#IMP.atom.show_molecular_hierarchy(t)
c=IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
if c.get_number_of_children()==0:
IMP.atom.show_molecular_hierarchy(t)
# there is no reason to use all atoms, just approximate the pdb shape instead
resolution/2.0)
# make the simplified structure rigid
rb.set_coordinates_are_optimized(True)
return s
if len(files) >1:
p= IMP.Particle(m)
h.set_name(name)
for i, f in enumerate(files):
c=create_from_pdb(f)
h.add_child(c)
c.set_name(name+" chain "+str(i))
for c in h.get_children()],
k)
if r:
m.add_restraint(r)
display_restraints.append(r)
m.set_maximum_score(r, k)
else:
h= create_from_pdb(files[0])
h.set_name(name)
all.add_child(h)
create_protein("Nup85", 570)
ct= IMP.atom.Selection(all, molecule="Nup85", terminus= IMP.atom.Selection.C)
d= IMP.core.XYZ(ct.get_selected_particles()[0])
d.set_coordinates(IMP.algebra.Vector3D(0,0,0))
d.set_coordinates_are_optimized(False)
create_protein("Nup84", 460)
create_protein("Nup145C", 442)
create_protein("Nup120", [0, 500, 761])
create_protein("Nup133", [0, 450, 778, 1160])
create_protein_from_pdbs("Seh1", ["seh1.pdb"])
create_protein_from_pdbs("Sec13", ["sec13.pdb"])
return (m, all)
def create_restraints(m, all):
def add_connectivity_restraint(s):
m.add_restraint(r)
m.set_maximum_score(r, k)
display_restraints.append(r)
def add_distance_restraint(s0, s1):
r=IMP.atom.create_distance_restraint(s0,s1, 0, k)
m.add_restraint(r)
m.set_maximum_score(r, k)
display_restraints.append(r)
evr=IMP.atom.create_excluded_volume_restraint([all])
m.add_restraint(evr)
s0=S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)])
s1=S(hierarchy=all, molecule="Nup84", molecules=[])
s2=S(hierarchy=all, molecule="Sec13")
add_connectivity_restraint([s0,s1,s2])
add_distance_restraint(S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Nup85"))
add_distance_restraint(S(hierarchy=all, molecule="Nup145C", residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Nup120",
residue_indexes= [(500, 762)]))
add_distance_restraint(S(hierarchy=all, molecule="Nup84"),
S(hierarchy=all, molecule="Nup133",
residue_indexes=[(778, 1160)]))
add_distance_restraint(S(hierarchy=all, molecule="Nup85"),
S(hierarchy=all, molecule="Seh1"))
add_distance_restraint(S(hierarchy=all, molecule="Nup145C",
residue_indexes=[(0,423)]),
S(hierarchy=all, molecule="Sec13"))
r= IMP.example.ExampleRestraint(l, k)
m.add_restraint(r)
# make sure rigid bodies are as not all particles in them can be on the x,y plane
m.set_maximum_score(.5*resolution**2*k)
# find acceptable conformations of the model
def get_conformations(m, gs):
sampler= IMP.core.MCCGSampler(m)
sampler.set_bounding_box(bb)
# magic numbers, experiment with them and make them large enough for things to work
sampler.set_number_of_conjugate_gradient_steps(100)
sampler.set_number_of_monte_carlo_steps(10)
sampler.set_number_of_attempts(20)
sampler.set_save_rejected_configurations(True)
# We don't care to see the output from the sampler
sampler.set_log_level(IMP.SILENT)
# return the IMP.ConfigurationSet storing all the found configurations that
# meet the various restraint maximum scores.
cs= sampler.get_sample()
# Look at the rejected minimal conformations to understand how the restraints
# are failing
print "rejected", \
sampler.get_rejected_configurations().get_number_of_configurations(), "solutions"
m.set_gather_statistics(True)
for i in range(0, sampler.get_rejected_configurations().get_number_of_configurations()):
sampler.get_rejected_configurations().load_configuration(i)
m.evaluate(False)
# show the score for each restraint to see which restraints were causing the
# conformation to be rejected
m.show_restraint_score_statistics()
w= IMP.display.PymolWriter("rejected.%d.pym"%i)
for g in gs:
w.add_geometry(g)
return cs
# cluster the conformations and write them to a file
def analyze_conformations(cs, all, gs):
# we want to cluster the configurations to make them easier to understand
# in the case, the clustering is pretty meaningless
embed= IMP.statistics.ConfigurationSetXYZEmbedding(cs,
cluster= IMP.statistics.create_lloyds_kmeans(embed, 10, 10000)
# dump each cluster center to a file so it can be viewed.
for i in range(cluster.get_number_of_clusters()):
rep= cluster.get_cluster_representative(i)
cs.load_configuration(rep)
w= IMP.display.PymolWriter("cluster.%d.pym"%i)
for g in gs:
w.add_geometry(g)
# now do the actual work
(m,all)= create_representation()
IMP.atom.show_molecular_hierarchy(all)
create_restraints(m, all)
gs=[]
for i in range(all.get_number_of_children()):
color= IMP.display.get_display_color(i)
n= all.get_child(i)
name= n.get_name()
g= IMP.core.XYZRGeometry(l)
g.set_color(color)
g.set_name(name)
gs.append(g)
# also display the restraints to see which particles they connect
for r in display_restraints:
try:
g= IMP.display.create_restraint_geometry(r)
gs.append(g)
except:
pass
cs= get_conformations(m, gs)
print "found", cs.get_number_of_configurations(), "solutions"
for i in range(0, cs.get_number_of_configurations()):
cs.load_configuration(i)
w= IMP.display.PymolWriter("config.%d.pym"%i)
for g in gs:
w.add_geometry(g)
analyze_conformations(cs, all, gs)
|
|
While we do not recommend doing serious work using restraints written in Python, it is often useful when prototyping or testing code. Copy this example and modify as needed.
import IMP
# a restraint which checks if particles are sorted in
# increasing order on k.
class MyRestraint(IMP.Restraint):
# take the list of particles and the key to use
def __init__(self, ps, k):
IMP.Restraint.__init__(self)
self.ps=ps
self.k=k
def unprotected_evaluate(self, da):
score=0
for i in range(1, len(self.ps)):
p0= self.ps[i-1]
p1= self.ps[i]
if p0.get_value(k) > p1.get_value(k):
diff=(p0.get_value(k)-p1.get_value(k))
score= score+diff
p0.add_to_derivative(k, -1, da)
p1.add_to_derivative(k, 1, da)
else:
if IMP.get_log_level() >= IMP.TERSE:
print p0.get_name(), "and", p1.get_name(), " are ok"
return score
def get_input_particles(self):
return self.ps
def get_input_containers(self):
return []
# some code to create and evaluate it
k= IMP.FloatKey("a key")
m= IMP.Model()
ps=[]
for i in range(0,10):
p = IMP.Particle(m)
p.add_attribute(k, i)
ps.append(p)
r= MyRestraint(ps, k)
m.add_restraint(r)
#IMP.set_log_level(IMP.TERSE)
print r.evaluate(True)
# this is needed to clean up properly, for some reason
m.remove_restraint(r)
del r
del m
|
While we do not recomment doing serious work using optimizer states written it python, it is often useful when prototyping or testing code. Copy this example and modify as needed.
import IMP
# an optimizer state which prints out model statistics.
class MyOptimizerState(IMP.OptimizerState):
def __init__(self):
IMP.OptimizerState.__init__(self)
def update(self):
self.get_optimizer().get_model().show_restraint_score_statistics()
# some code to create and evaluate it
k= IMP.FloatKey("a key")
m= IMP.Model()
# we don't have any real restraints in the kernel
r0=IMP._ConstRestraint(1)
r0.set_name("restraint 0")
m.add_restraint(r0)
r1=IMP._ConstRestraint(2)
r1.set_name("restraint 1")
m.add_restraint(r1)
os= MyOptimizerState()
os.set_name("python optimizer state")
# we don't have any optimizers either
co= IMP._ConstOptimizer(m)
co.add_optimizer_state(os)
m.set_gather_statistics(True)
# so we only see the statistics
IMP.set_log_level(IMP.SILENT)
print co.optimize(100)
# this is needed to clean up memory properly for some reason
co.remove_optimizer_state(os)
del os
del m
|
|