<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>