 |
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() ps= IMP.core.DistancePairScore(IMP.core.HarmonicLowerBound(1,1)) 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) lr=IMP.container.PairsRestraint(IMP.core.SoftSpherePairScore(1), nbl, "excluded volume") m.add_restraint(lr) # Tie the ends of the chain tie=IMP.core.PairRestraint(IMP.core.HarmonicDistancePairScore(3,1), (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) d=IMP.display.ChimeraWriter("solution"+str(i)+".py") for p in chain.get_particles(): d.add_geometry(IMP.display.XYZRGeometry(p)) # print out info about used modules so that the versions are known IMP.show_used_modules() |
|
Show how to add a handler so that the state of the model is dumped when an error occurs.
import IMP import IMP.core m= IMP.Model() pts= IMP.core.create_xyzr_particles(m, 10, 10, 10) IMP.add_failure_handler(IMP.WriteParticlesFailureHandler(pts, "error.imp")) |
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(0,0,0), IMP.algebra.Vector3D(300, 300, 300)) # this function creates the molecular hierarchies for the various involved proteins def create_representation(): m= IMP.Model() all=IMP.atom.Hierarchy.setup_particle(IMP.Particle(m)) 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) leaves= IMP.atom.get_leaves(h) # for convenience, have one molecular hierarchy containing all molecules all.add_child(h) r=IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c)\ for c in h.get_children()], k) if r: m.add_restraint(r) # only allow the particles to penetrate or 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 score= IMP.core.KClosePairsPairScore(IMP.core.HarmonicSphereDistancePairScore(0,1), 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 S= IMP.atom.Selection 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, IMP.container.ListSingletonContainer(IMP.atom.get_leaves(all)), True) 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()): center= cluster.get_cluster_center(i) cs.load_configuration(i) 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) # 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.display.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= IMP.display.HierarchyGeometry(n) 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 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) |
|
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=IMP.atom.Hierarchy.setup_particle(IMP.Particle(m)) 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) r=IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c)\ 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.SetLogState(IMP.NONE) t=IMP.atom.read_pdb( IMP.get_example_path("data/"+file), m, IMP.atom.ATOMPDBSelector()) 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 s=IMP.atom.create_simplified_along_backbone(c, resolution/2.0) IMP.atom.destroy(t) # make the simplified structure rigid rb=IMP.atom.create_rigid_body(s) rb.set_coordinates_are_optimized(True) return s if len(files) >1: p= IMP.Particle(m) h= IMP.atom.Hierarchy.setup_particle(p) 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)) r=IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c)\ 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): r= IMP.atom.create_connectivity_restraint(s, k) 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) S= IMP.atom.Selection 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) # 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, IMP.container.ListSingletonContainer(IMP.atom.get_leaves(all)), True) 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()): center= cluster.get_cluster_center(i) cs.load_configuration(i) 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() for l in IMP.atom.get_leaves(n): g= IMP.display.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) |
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() # 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) |
|