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.
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 ...
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(100, 100, 100))
# 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,
1)
#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 # <------- swapping c with s will give a coarse grain representation - much faster !
# 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)
#ub = IMP.core.HarmonicUpperBound(1.0, 0.1)
#ss = IMP.core.DistancePairScore(ub)
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])
n1 = r.add_composite([iA, iB, iC, iD])
n2 = r.add_composite([iA, iB],n1)
n3 = r.add_composite([iB, iD],n1)
n4 = r.add_composite([iA, iB, iC],n1)
n5 = 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])
nbl = IMP.container.ClosePairContainer([all], 0, 2)
h = IMP.core.HarmonicLowerBound(0, 1)
sd = IMP.core.SphereDistancePairScore(h)
# use the lower bound on the inter-sphere distance to push the spheres apart
nbr = IMP.container.PairsRestraint(sd, nbl)
m.add_restraint(nbr)
# r1 = IMP.core.ExcludedVolumeRestraint(all)
# m.add_restraint(r1)
# 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(200)
sampler.set_number_of_monte_carlo_steps(40)
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 this 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)
h = IMP.atom.Hierarchy.get_children(all)
#tfn = IMP.create_temporary_file_name("cluster%d"%i, ".rmf")
huh = "./models/CLUSTER%d"%i
huh = huh +".rmf"
#print "file is", tfn
print "file is", huh
rh = RMF.create_rmf_file(huh)
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:
# rh.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))
# 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)
h = IMP.atom.Hierarchy.get_children(all)
#tfn = IMP.create_temporary_file_name("josh%d"%i, ".rmf")
#print "file is", tfn
huh = "./models/IMP%d"%i
huh = huh +".rmf"
print "file is", huh
rh = RMF.create_rmf_file(huh)
# 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)