[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [IMP-users] getting the MCCGsampler to work (Josh Bullock)



hi list, sorry, just wanted to clarify a little more my ( lack of ) understanding re particles ...

so when i create a rigid body, i'm decorating a selection of particles - but it's the particles that hold all the information. and what i want to add to the root hierarchy are the decorated particles, not the decorator itself, correct ? i think that this is why i was getting no mass before.

so when i load a pdb into IMP, a particle is created for each atom ? I've been loading 4 pdb structures into IMP with a total number of 308 amino acids, but i seem to be creating at least 10,515 particles ( shown by printing IMP.kernel.Particle.get_index(IMP.Particle(m))Â) and i can't make that add up - conceptually i'm a bit lost ...

and finally is there a way to access each set of decorated particles iteratively ? i think my problems with creating the RigidBodyMover are because i'm trying to add the entire model at once ( as shown ) - and there are some particles in there which i haven't been able to decorate.Â

rb_mover = IMP.core.RigidBodyMover(m, IMP.kernel.Particle.get_index(IMP.Particle(m)), maxTranslation, maxRotation )

eagerly anticipating enlightenment,

joshÂ

p.s. full script for referenceÂ
https://gist.github.com/mysticvision/abb8d42e24f329b5388b

On 14 July 2014 17:09, <" target="_blank">> wrote:
Send IMP-users mailing list submissions to
    ">

To subscribe or unsubscribe via the World Wide Web, visit
    https://salilab.org/mailman/listinfo/imp-users
or, via email, send a message with subject or body 'help' to
    ">

You can reach the person managing the list at
    ">

When replying, please edit your Subject line so it is more specific
than "Re: Contents of IMP-users digest..."


Today's Topics:

 Â1. Re: getting the MCCGsampler to work (Josh Bullock)


----------------------------------------------------------------------

Message: 1
Date: Mon, 14 Jul 2014 17:09:08 +0100
From: Josh Bullock <">>
To: ">
Subject: Re: [IMP-users] getting the MCCGsampler to work
Message-ID:
    <">CAHh_40-ROvuH3Koy+>
Content-Type: text/plain; charset="utf-8"

ah i see, yes i had no movers or optimizers ...

so now i'm trying to set up a RigidBodyMover but i'm struggling to create
the rigid bodies. There seem to be two different ways of doing this, either
IMP.atom.create_rigid_body or IMP.atom.setup_as_rigid_body.

If I use create_rigid_body then when i later try to add restraints i get
told "leaf 'A' ... does not have mass"

If I use setup_as_rigid_body i can add the restraints but i when i try to
add it to the RigidBodyMover i get told "Rigid body passed to
RigidBodyMover must be set to be optimized" even
though rb.set_coordinates_are_optimized(True)

what is the function of choice here ?

thanks !

josh

------------------------------

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

      # there is no reason to use all atoms, just approximate the pdb
shape instead
      s=IMP.atom.create_simplified_along_backbone(c,1)

      # make the simplified structure rigid
      rb=IMP.atom.create_rigid_body(s)
      #rb=IMP.atom.setup_as_rigid_body(s)
      #rb=IMP.atom.create_rigid_body(c)
      rb.set_coordinates_are_optimized(True)
      print rb.get_coordinates_are_optimized()

      return rb

    h= create_from_pdb(files[0])
    h.set_name(name)
    all.add_child(h)

  create_protein_from_pdbs("A", [Firstpdb])

(m,all)= create_representation()


On 11 July 2014 19:42, <">> wrote:

> Send IMP-users mailing list submissions to
> Â Â Â Â ">
>
> To subscribe or unsubscribe via the World Wide Web, visit
> Â Â Â Â https://salilab.org/mailman/listinfo/imp-users
> or, via email, send a message with subject or body 'help' to
> Â Â Â Â ">
>
> You can reach the person managing the list at
> Â Â Â Â ">
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of IMP-users digest..."
>
>
> Today's Topics:
>
> Â Â1. getting the MCCGsampler to work (Josh Bullock)
> Â Â2. Re: getting the MCCGsampler to work (Ben Webb)
> Â Â3. Re: getting the MCCGsampler to work (Ben Webb)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Fri, 11 Jul 2014 13:12:01 +0100
> From: Josh Bullock <">>
> To: ">
> Subject: [IMP-users] getting the MCCGsampler to work
> Message-ID:
> Â Â Â Â <
> ">>
> Content-Type: text/plain; charset="utf-8"
>
> Hello again,
>
> 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 ...
>
> 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(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)
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://salilab.org/archives/imp-users/attachments/20140711/fc50f2fe/attachment.html
> >
>
> ------------------------------
>
> Message: 2
> Date: Fri, 11 Jul 2014 11:26:24 -0700
> From: Ben Webb <">>
> To: Help and discussion for users of IMP <">>
> Subject: Re: [IMP-users] getting the MCCGsampler to work
> Message-ID: <">>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> On 7/11/14, 5:12 AM, Josh Bullock wrote:
> > 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.
>
> A sampler that ignores the scoring function would be pretty useless!
> Obviously MCCG doesn't do that. So my guess would be that you didn't add
> the restraints to the scoring function, the thresholds are all wrong (so
> that no MC moves are ever accepted), you have terrible starting
> conditions (so that the optimizer can never improve the score) or you
> have a poor choice of MC move set. I don't see you actually adding any
> Movers to your sampler in the script, so I'd suspect the last one. The
> underlying MonteCarlo class keeps some statistics (e.g. number of
> accepted/rejected moves) which might shed some light on what's going on.
>
> > So I saw on an old (2011) nup84 example that MCCG can't handle rigid
> > bodies, is this still the case ?
>
> No. You just need to add a suitable RigidBodyMover.
>
> > If so, should I switch to the DOMINO
> > sampler ?
>
> DOMINO enumerates in a discrete space, so is a rather different kind of
> optimizer. It is probably not the best choice in your case.
>
> Â Â Â Â Ben
> --
> "> Â Â Â Â Â Â Â Â Â Â Âhttp://salilab.org/~ben/
> "It is a capital mistake to theorize before one has data."
> Â Â Â Â - Sir Arthur Conan Doyle
>
>
> ------------------------------
>
> Message: 3
> Date: Fri, 11 Jul 2014 11:42:19 -0700
> From: Ben Webb <">>
> To: Help and discussion for users of IMP <">>
> Subject: Re: [IMP-users] getting the MCCGsampler to work
> Message-ID: <">>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> To clarify a little, the MCCGSampler tries to wrap a common use of the
> underlying MonteCarlo optimizer (that use case being a bunch of
> independently movable x,y,z particles). If in your case you actually
> have a bunch of rigid bodies, you should probably use the MonteCarlo
> optimizer directly and add RigidBodyMovers to it.
>
> Â Â Â Â Ben
> --
> "> Â Â Â Â Â Â Â Â Â Â Âhttp://salilab.org/~ben/
> "It is a capital mistake to theorize before one has data."
> Â Â Â Â - Sir Arthur Conan Doyle
>
>
> ------------------------------
>
> _______________________________________________
> IMP-users mailing list
> ">
> https://salilab.org/mailman/listinfo/imp-users
>
>
> End of IMP-users Digest, Vol 38, Issue 15
> *****************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://salilab.org/archives/imp-users/attachments/20140714/3720260b/attachment.html>

------------------------------

_______________________________________________
IMP-users mailing list
">
https://salilab.org/mailman/listinfo/imp-users


End of IMP-users Digest, Vol 38, Issue 16
*****************************************