Hello,
I'm relatively new to all this so please let me know if i'm making any obvious errors ...
Essentially all i'm trying to do is generate an ensemble of models made from four subunits - constrained by MS connectivity restraints. The models get scored but nothing seems to write to the pymol file. Ideally i'd like to write to an .rmf but i haven't worked that one out either ...
Is this a reasonable way to go about my problem ?
Many thanks,
Josh
-------------------------------------------
import IMP import IMP.atom import IMP.rmf import inspect import IMP.container import IMP.display import IMP.statistics #import IMP.example import sys, math, os, optparse import RMF
from optparse import OptionParser
# Convert the arguments into strings and number Firstpdb = str(sys.argv[1]) Secondpdb = str(sys.argv[2]) Thirdpdb = str(sys.argv[3]) Fourthpdb = str(sys.argv[4]) models = float(sys.argv[5])
#*****************************************
# the spring constant to use, it doesnt 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 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_from_pdbs(name, files):
def create_from_pdb(file): sls=IMP.SetLogState(IMP.NONE) datadir = os.getcwd() print datadir t=IMP.atom.read_pdb( datadir+'/' + 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/300.0) IMP.atom.destroy(t) # make the simplified structure rigid rb=IMP.atom.create_rigid_body(s) # rb=IMP.atom.create_rigid_body(c) rb.set_coordinates_are_optimized(True) return s # return c
h= create_from_pdb(files[0]) h.set_name(name) all.add_child(h)
create_protein_from_pdbs("A", [Firstpdb]) create_protein_from_pdbs("B", [Secondpdb]) create_protein_from_pdbs("C", [Thirdpdb]) create_protein_from_pdbs("D", [Fourthpdb]) #create_protein_from_pdbs("C", ["rpt3_imp.pdb"]) 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.MSConnectivityRestraint(m,score)
iA = r.add_type([rps[0]]) iB = r.add_type([rps[1]]) iC = r.add_type([rps[2]]) iD = r.add_type([rps[3]]) n1 = r.add_composite([iA, iB, iC, iD]) n2 = r.add_composite([iA, iB], n1) n3 = r.add_composite([iC, iD], n1) n4 = r.add_composite([iB, iC, iD], n1)
m.add_restraint(r)
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 sA=S(hierarchy=all, molecule="A") sB=S(hierarchy=all, molecule="B") sC=S(hierarchy=all, molecule="C") sD=S(hierarchy=all, molecule="D") add_connectivity_restraint([sA, sB, sC, sD])
# 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(20) sampler.set_number_of_attempts(models) # 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.create_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.atom.HierarchyGeometry(n) g.set_color(color) gs.append(g)
cs= get_conformations(m)
print "found", cs.get_number_of_configurations(), "solutions"
ListScores = [] for i in range(0, cs.get_number_of_configurations()): cs.load_configuration(i) # print the configuration print "solution number: ",i,"scored :", m.evaluate(False) ListScores.append(m.evaluate(False))
f1 = open("out_scores.csv", "w") f1.write("\n".join(map(lambda x: str(x), ListScores))) f1.close()
# for each of the configuration, dump it to a file to view in pymol for i in range(0, cs.get_number_of_configurations()): JOSH = cs.load_configuration(i) S= IMP.atom.Selection h= IMP.atom.Hierarchy.get_children(cs) tfn = IMP.create_temporary_file_name("josh%d"%i, ".rmf") rh = RMF.create_rmf_file(tfn)
# add the hierarchy to the file IMP.rmf.add_hierarchies(rh, h)
# add the current configuration to the file as frame 0 IMP.rmf.save_frame(rh)
for g in gs: w.add_geometry(g)
analyze_conformations(cs, all, gs)