<div dir="ltr">Hello again,<div><br></div><div>So I have the same overall problem as before - creating an ensemble of a 4-subunit complex using MSconnectivity restraints. Having visualised the output (via RMF - thanks Barak :¬) , it's clear that no matter how many steps of CG or MC I put, the models do not change from their initial random placement. I know that the restraints are present because the models are evaluated and scored appropriately.</div>
<div><br></div><div>So I saw on an old (2011) nup84 example that MCCG can't handle rigid bodies, is this still the case ? If so, should I switch to the DOMINO sampler ? or it does work and likely there's an error in my code ...</div>
<div><br></div><div>Many thanks !</div><div><br></div><div>Josh</div><div><br></div><div>-------------------------------------------------</div><div><br></div><div><div>import IMP</div><div>import IMP.atom</div><div>import IMP.rmf</div>
<div>import inspect</div><div>import IMP.container</div><div>import IMP.display</div><div>import IMP.statistics</div><div>#import IMP.example</div><div>import sys, math, os, optparse</div><div>import RMF</div><div><br></div>
<div>from optparse import OptionParser</div><div><br></div><div><br></div><div># Convert the arguments into strings and number</div><div>Firstpdb = str(sys.argv[1])</div><div>Secondpdb = str(sys.argv[2])</div><div>Thirdpdb = str(sys.argv[3])</div>
<div>Fourthpdb = str(sys.argv[4])</div><div>models = float(sys.argv[5])</div><div><br></div><div>#***************************************** </div><div><br></div><div># the spring constant to use, it doesnt really matter</div>
<div>k=100</div><div># the target resolution for the representation, this is used to specify how detailed</div><div># the representation used should be</div><div>resolution=300</div><div># the box to perform everything </div>
<div>bb=IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(0,0,0),</div><div> IMP.algebra.Vector3D(100, 100, 100))</div><div><br></div><div><br></div><div># this function creates the molecular hierarchies for the various involved proteins</div>
<div>def create_representation():</div><div> m= IMP.Model()</div><div> all=IMP.atom.Hierarchy.setup_particle(IMP.Particle(m))</div><div> all.set_name("the universe")</div><div> # create a protein, represented as a set of connected balls of appropriate</div>
<div> # radii and number, chose by the resolution parameter and the number of</div><div> # amino acids.</div><div> </div><div> def create_protein_from_pdbs(name, files):</div><div> </div><div> def create_from_pdb(file):</div>
<div> sls=IMP.SetLogState(IMP.NONE)</div><div> datadir = os.getcwd()</div><div> print datadir</div><div><span class="" style="white-space:pre"> </span> t=IMP.atom.read_pdb( datadir+'/' + file, m,</div>
<div> IMP.atom.ATOMPDBSelector())</div><div> del sls</div><div> #IMP.atom.show_molecular_hierarchy(t)</div><div> c=IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])</div>
<div> if c.get_number_of_children()==0:</div><div> IMP.atom.show_molecular_hierarchy(t)</div><div> # there is no reason to use all atoms, just approximate the pdb shape instead</div><div>
s=IMP.atom.create_simplified_along_backbone(c,</div><div> 1)</div><div> #IMP.atom.destroy(t)</div><div> # make the simplified structure rigid</div>
<div> rb=IMP.atom.create_rigid_body(s) </div><div> rb=IMP.atom.create_rigid_body(c)</div><div> rb.set_coordinates_are_optimized(True)</div><div> return s # <------- swapping c with s will give a coarse grain representation - much faster !</div>
<div># return c </div><div><br></div><div> h= create_from_pdb(files[0])</div><div> h.set_name(name)</div><div> all.add_child(h)</div><div><br></div><div> create_protein_from_pdbs("A", [Firstpdb])</div>
<div> create_protein_from_pdbs("B", [Secondpdb])</div><div> create_protein_from_pdbs("C", [Thirdpdb])</div><div> create_protein_from_pdbs("D", [Fourthpdb])</div><div> #create_protein_from_pdbs("C", ["rpt3_imp.pdb"])</div>
<div> return (m, all)</div><div><br></div><div># create the needed restraints and add them to the model</div><div><br></div><div>def create_restraints(m, all):</div><div> def add_connectivity_restraint(s):</div><div>
</div><div> tr= IMP.core.TableRefiner()</div><div> rps=[]</div><div> for sc in s:</div><div> ps= sc.get_selected_particles() </div><div> rps.append(ps[0])</div><div> tr.add_particle(ps[0], ps)</div>
<div> </div><div> # duplicate the IMP.atom.create_connectivity_restraint functionality</div><div> </div><div> score= IMP.core.KClosePairsPairScore(IMP.core.HarmonicSphereDistancePairScore(0,1),tr)</div>
<div> </div><div> #ub = IMP.core.HarmonicUpperBound(1.0, 0.1)</div><div> #ss = IMP.core.DistancePairScore(ub)</div><div> </div><div> r= IMP.core.MSConnectivityRestraint(m,score)</div><div>
</div><div> iA = r.add_type([rps[0]])</div><div> iB = r.add_type([rps[1]])</div><div> iC = r.add_type([rps[2]])</div><div> iD = r.add_type([rps[3]])</div><div> #n1 = r.add_composite([iA, iB])</div>
<div> n1 = r.add_composite([iA, iB, iC, iD])</div><div> n2 = r.add_composite([iA, iB],n1)</div><div> n3 = r.add_composite([iB, iD],n1)</div><div> n4 = r.add_composite([iA, iB, iC],n1)</div><div>
n5 = r.add_composite([iB, iC, iD],n1)</div><div><br></div><div> m.add_restraint(r)</div><div><br></div><div> evr=IMP.atom.create_excluded_volume_restraint([all])</div><div> m.add_restraint(evr)</div><div>
# a Selection allows for natural specification of what the restraints act on</div><div> S= IMP.atom.Selection</div><div> sA=S(hierarchy=all, molecule="A")</div><div> sB=S(hierarchy=all, molecule="B")</div>
<div> sC=S(hierarchy=all, molecule="C")</div><div> sD=S(hierarchy=all, molecule="D")</div><div> add_connectivity_restraint([sA, sB, sC, sD])</div><div> </div><div> nbl = IMP.container.ClosePairContainer([all], 0, 2)</div>
<div> h = IMP.core.HarmonicLowerBound(0, 1)</div><div> sd = IMP.core.SphereDistancePairScore(h)</div><div> # use the lower bound on the inter-sphere distance to push the spheres apart</div><div> nbr = IMP.container.PairsRestraint(sd, nbl)</div>
<div> m.add_restraint(nbr)</div><div> </div><div> </div><div> # r1 = IMP.core.ExcludedVolumeRestraint(all)</div><div> # m.add_restraint(r1)</div><div> </div><div><br></div><div># find acceptable conformations of the model</div>
<div>def get_conformations(m):</div><div> sampler= IMP.core.MCCGSampler(m)</div><div> sampler.set_bounding_box(bb)</div><div> # magic numbers, experiment with them and make them large enough for things to work</div>
<div> sampler.set_number_of_conjugate_gradient_steps(200)</div><div> sampler.set_number_of_monte_carlo_steps(40)</div><div> sampler.set_number_of_attempts(models)</div><div> # We don't care to see the output from the sampler</div>
<div> #sampler.set_log_level(IMP.SILENT)</div><div> # return the IMP.ConfigurationSet storing all the found configurations that</div><div> # meet the various restraint maximum scores.</div><div> cs= sampler.create_sample()</div>
<div> return cs</div><div> </div><div><br></div><div># cluster the conformations and write them to a file</div><div>def analyze_conformations(cs, all, gs):</div><div> # we want to cluster the configurations to make them easier to understand</div>
<div> # in this case, the clustering is pretty meaningless</div><div> embed= IMP.statistics.ConfigurationSetXYZEmbedding(cs,</div><div> IMP.container.ListSingletonContainer(IMP.atom.get_leaves(all)), True)</div>
<div> cluster= IMP.statistics.create_lloyds_kmeans(embed, 10, 10000)</div><div> # dump each cluster center to a file so it can be viewed.</div><div> for i in range(cluster.get_number_of_clusters()):</div><div> center= cluster.get_cluster_center(i)</div>
<div> cs.load_configuration(i)</div><div> h = IMP.atom.Hierarchy.get_children(all)</div><div> #tfn = IMP.create_temporary_file_name("cluster%d"%i, ".rmf")</div><div> huh = "./models/CLUSTER%d"%i</div>
<div> huh = huh +".rmf"</div><div> #print "file is", tfn</div><div> print "file is", huh</div><div> rh = RMF.create_rmf_file(huh)</div><div> </div><div> </div>
<div> IMP.rmf.add_hierarchies(rh, h)</div><div> </div><div> # add the current configuration to the file as frame 0</div><div> IMP.rmf.save_frame(rh)</div><div> </div><div> #for g in gs:</div>
<div> # rh.add_geometry(g)</div><div><br></div><div><br></div><div>#******************************************************************************************</div><div># now do the actual work</div><div><br></div>
<div>(m,all)= create_representation()</div><div>#IMP.atom.show_molecular_hierarchy(all)</div><div>create_restraints(m, all)</div><div><br></div><div># in order to display the results, we need something that maps the particles onto</div>
<div># geometric objets. The IMP.display.Geometry objects do this mapping.</div><div># IMP.display.XYZRGeometry map an IMP.core.XYZR particle onto a sphere</div><div>gs=[]</div><div>for i in range(all.get_number_of_children()):</div>
<div> color= IMP.display.get_display_color(i)</div><div> n= all.get_child(i)</div><div> name= n.get_name()</div><div> g= IMP.atom.HierarchyGeometry(n)</div><div> g.set_color(color)</div><div> gs.append(g)</div>
<div><br></div><div>cs= get_conformations(m)</div><div><br></div><div>print "found", cs.get_number_of_configurations(), "solutions"</div><div><br></div><div>ListScores = []</div><div>for i in range(0, cs.get_number_of_configurations()):</div>
<div> cs.load_configuration(i)</div><div> # print the configuration</div><div> print "solution number: ",i,"scored :", m.evaluate(False)</div><div> ListScores.append(m.evaluate(False))</div>
<div><br></div><div># for each of the configuration, dump it to a file to view in pymol</div><div>for i in range(0, cs.get_number_of_configurations()):</div><div> cs.load_configuration(i)</div><div> h = IMP.atom.Hierarchy.get_children(all)</div>
<div> #tfn = IMP.create_temporary_file_name("josh%d"%i, ".rmf")</div><div> #print "file is", tfn</div><div> huh = "./models/IMP%d"%i</div><div> huh = huh +".rmf"</div>
<div> print "file is", huh</div><div> rh = RMF.create_rmf_file(huh)</div><div> </div><div> # add the hierarchy to the file</div><div> IMP.rmf.add_hierarchies(rh, h)</div><div> </div><div> # add the current configuration to the file as frame 0</div>
<div> IMP.rmf.save_frame(rh)</div><div> </div><div> #for g in gs:</div><div> # w.add_geometry(g)</div><div><br></div><div>analyze_conformations(cs, all, gs)</div></div><div><br></div></div>