 |
The script shows how to assess a protein conformation using DOPE.
import IMP
import IMP.atom
import IMP.container
def create_representation():
m=IMP.Model()
mp0= IMP.atom.read_pdb(IMP.atom.get_example_path('1fdx.B99990001.pdb'), m, IMP.atom.NonWaterNonHydrogenPDBSelector())
prot=IMP.atom.get_by_type(mp0, IMP.atom.CHAIN_TYPE)[0]
return (m, prot)
def add_dope(m, prot):
ps=IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
dsc.add_particles(ps)
dpc = IMP.container.ClosePairContainer(dsc, 15.0, 0.0)
# exclude pairs of atoms belonging to the same residue
# for consistency with MODELLER DOPE score
dpc.add_pair_filter(f)
dps= IMP.atom.DopePairScore(15.0)
# dps= IMP.membrane.DopePairScore(15.0, IMP.membrane.get_data_path("dope_scorehr.lib"))
d= IMP.container.PairsRestraint(dps, dpc)
m.add_restraint(d)
print "creating representation"
(m,prot)=create_representation()
print "creating DOPE score function"
add_dope(m,prot)
IMP.set_check_level(IMP.USAGE)
print "DOPE SCORE ::",m.evaluate(False)
|
This example prints out various statistics about a prospective Brownian dynamics simulation. You can use the statistics to determine if the time step is likely to be sufficiently short given the forces and particle sizes involved.
# fs
time_step=1000
# angstrom
minimum_particle_radius=10
# kCal/mol/A
maximum_spring_constant=1
maximum_diffusion_coefficient=\
get_einstein_diffusion_coefficient(minimum_particle_radius)
expected_delta= get_diffusion_length(maximum_diffusion_coefficient,
time_step)
expected_rotational_delta=\
get_diffusion_angle(maximum_diffusion_coefficient,
time_step) * minimum_particle_radius
expected_spring_diffusion_length=\
get_diffusion_length(maximum_diffusion_coefficient,
.5*maximum_spring_constant*4*expected_delta**2,
time_step)
print "with a time step of", time_step, "fs"
print "an object of radius", minimum_particle_radius, "A will move",\
expected_delta, "A and a point on its surface will move",\
expected_rotational_delta,"A more"
print "the motion from fluctuations in the spring compression will be",\
expected_spring_diffusion_length, "A"
print "and a compression of 10% of the radius will induce a motion of",\
get_diffusion_length(maximum_diffusion_coefficient,
.5*maximum_spring_constant*\
(.1*minimum_particle_radius)**2,
time_step), "A"
|
This example shows how to properly write out a pdb of the structure results when using coarse grained rigid bodies.
import IMP.atom
m= IMP.Model()
full=IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m)
chain= IMP.atom.get_by_type(full, IMP.atom.CHAIN_TYPE)[0]
print chain
# for some reason the python wrapper won't make the implicit conversion to Chain
simplified= IMP.atom.create_simplified_along_backbone(IMP.atom.Chain(chain), 3)
IMP.atom.destroy(full)
rb= IMP.atom.create_rigid_body(simplified)
original_transform= rb.get_reference_frame().get_transformation_to()
# fake optimization, just move the rigid body
rb.set_reference_frame(IMP.algebra.ReferenceFrame3D(IMP.algebra.Transformation3D(tr, IMP.algebra.Vector3D(0,0,0))))
# extract the difference
diff= rb.get_reference_frame().get_transformation_to()/original_transform
reload=IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m)
IMP.atom.transform(reload, diff)
name=IMP.create_temporary_file("out", ".pdb")
IMP.atom.write_pdb(reload, name)
|
In this example, a PDB file is read in and scored using the CHARMM forcefield. For more control over the setup of the forcefield, see the 'charmm_forcefield_verbose.py' example.
import IMP.atom
import IMP.container
# Create an IMP model and add a heavy atom-only protein from a PDB file
m = IMP.Model()
prot = IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m,
# Read in the CHARMM heavy atom topology and parameter files
# Using the CHARMM libraries, determine the ideal topology (atoms and their
# connectivity) for the PDB file's primary sequence
topology = ff.create_topology(prot)
# Typically this modifies the C and N termini of each chain in the protein by
# applying the CHARMM CTER and NTER patches. Patches can also be manually
# applied at this point, e.g. to add disulfide bridges.
topology.apply_default_patches()
# Make the PDB file conform with the topology; i.e. if it contains extra
# atoms that are not in the CHARMM topology file, remove them; if it is
# missing atoms (e.g. sidechains, hydrogens) that are in the CHARMM topology,
# add them and construct their Cartesian coordinates from internal coordinate
# information.
topology.setup_hierarchy(prot)
# Set up and evaluate the stereochemical part (bonds, angles, dihedrals,
# impropers) of the CHARMM forcefield
r = IMP.atom.CHARMMStereochemistryRestraint(prot, topology)
m.add_restraint(r)
# Add non-bonded interaction (in this case, Lennard-Jones). This needs to
# know the radii and well depths for each atom, so add them from the forcefield
# (they can also be assigned manually using the XYZR or LennardJones
# decorators):
ff.add_radii(prot)
ff.add_well_depths(prot)
# Get a list of all atoms in the protein, and put it in a container
atoms = IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
cont = IMP.container.ListSingletonContainer(atoms)
# Add a restraint for the Lennard-Jones interaction. This is built from
# a collection of building blocks. First, a ClosePairContainer maintains a list
# of all pairs of Particles that are close. Next, all 1-2, 1-3 and 1-4 pairs
# from the stereochemistry created above are filtered out.
# Then, a LennardJonesPairScore scores a pair of atoms with the Lennard-Jones
# potential. Finally, a PairsRestraint is used which simply applies the
# LennardJonesPairScore to each pair in the ClosePairContainer.
nbl = IMP.container.ClosePairContainer(cont, 4.0)
nbl.add_pair_filter(r.get_pair_filter())
sf = IMP.atom.ForceSwitch(6.0, 7.0)
ps = IMP.atom.LennardJonesPairScore(sf)
m.add_restraint(IMP.container.PairsRestraint(ps, nbl))
# it gets awfully slow with internal checks
IMP.set_check_level(IMP.USAGE)
# Finally, evaluate the score of the whole system (without derivatives)
print m.evaluate(False)
|
|
In this example, a PDB file is read in and scored using the CHARMM forcefield. It is similar to the 'charmm_forcefield.py' example, but fully works through each step of the procedure using lower-level IMP classes. This is useful if you want to customize the way in which the forcefield is applied.
import IMP.atom
import IMP.container
# Create an IMP model and add a heavy atom-only protein from a PDB file
m = IMP.Model()
prot = IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m,
# Read in the CHARMM heavy atom topology and parameter files
# Using the CHARMM libraries, determine the ideal topology (atoms and their
# connectivity) for the PDB file's primary sequence
topology = ff.create_topology(prot)
# Typically this modifies the C and N termini of each chain in the protein by
# applying the CHARMM CTER and NTER patches. Patches can also be manually
# applied at this point, e.g. to add disulfide bridges.
topology.apply_default_patches()
# Each atom is mapped to its CHARMM type. These are needed to look up bond
# lengths, Lennard-Jones radii etc. in the CHARMM parameter file. Atom types
# can also be manually assigned at this point using the CHARMMAtom decorator.
topology.add_atom_types(prot)
# Remove any atoms that are in the PDB file but not in the topology, and add
# in any that are in the topology but not the PDB.
topology.add_missing_atoms(prot)
# Construct Cartesian coordinates for any atoms that were added
topology.add_coordinates(prot)
# Generate and return lists of bonds, angles, dihedrals and impropers for
# the protein. Each is a Particle in the model which defines the 2, 3 or 4
# atoms that are bonded, and adds parameters such as ideal bond length
# and force constant. Note that bonds and impropers are explicitly listed
# in the CHARMM topology file, while angles and dihedrals are generated
# automatically from an existing set of bonds. These particles only define the
# bonds, but do not score them or exclude them from the nonbonded list.
bonds = topology.add_bonds(prot)
angles = ff.create_angles(bonds)
dihedrals = ff.create_dihedrals(bonds)
impropers = topology.add_impropers(prot)
# Maintain stereochemistry by scoring bonds, angles, dihedrals and impropers
# Score all of the bonds. This is done by combining IMP 'building blocks':
# - A ListSingletonContainer simply manages a list of the bond particles.
# - A BondSingletonScore, when given a bond particle, scores the bond by
# calculating the distance between the two atoms it bonds, subtracting the
# ideal value, and weighting the result by the bond's "stiffness", such that
# an "ideal" bond scores zero, and bonds away from equilibrium score non-zero.
# It then hands off to a UnaryFunction to actually penalize the value. In
# this case, a Harmonic UnaryFunction is used with a mean of zero, so that
# bond lengths are harmonically restrained.
# - A SingletonsRestraint simply goes through each of the bonds in the
# container and scores each one in turn.
cont = IMP.container.ListSingletonContainer(bonds, "bonds")
bss = IMP.atom.BondSingletonScore(IMP.core.Harmonic(0, 1))
r=IMP.container.SingletonsRestraint(bss, cont, "bonds")
m.add_restraint(r)
# Score angles, dihedrals, and impropers. In the CHARMM forcefield, angles and
# impropers are harmonically restrained, so this is the same as for bonds.
# Dihedrals are scored internally by a periodic (cosine) function.
cont = IMP.container.ListSingletonContainer(angles, "angles")
bss = IMP.atom.AngleSingletonScore(IMP.core.Harmonic(0,1))
r=IMP.container.SingletonsRestraint(bss, cont, "angles")
m.add_restraint(r)
cont = IMP.container.ListSingletonContainer(dihedrals, "dihedrals")
bss = IMP.atom.DihedralSingletonScore()
r=IMP.container.SingletonsRestraint(bss, cont, "dihedrals")
m.add_restraint(r)
cont = IMP.container.ListSingletonContainer(impropers, "impropers")
bss = IMP.atom.ImproperSingletonScore(IMP.core.Harmonic(0,1))
m.add_restraint(IMP.container.SingletonsRestraint(bss, cont, "improppers"))
# Add non-bonded interaction (in this case, Lennard-Jones). This needs to
# know the radii and well depths for each atom, so add them from the forcefield
# (they can also be assigned manually using the XYZR or LennardJones
# decorators):
ff.add_radii(prot)
ff.add_well_depths(prot)
# Get a list of all atoms in the protein, and put it in a container
atoms = IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
cont = IMP.container.ListSingletonContainer(atoms)
# Add a restraint for the Lennard-Jones interaction. Again, this is built from
# a collection of building blocks. First, a ClosePairContainer maintains a list
# of all pairs of Particles that are close. A StereochemistryPairFilter is used
# to exclude atoms from this list that are bonded to each other or are involved
# in an angle or dihedral (1-3 or 1-4 interaction). Then, a
# LennardJonesPairScore scores a pair of atoms with the Lennard-Jones potential.
# Finally, a PairsRestraint is used which simply applies the
# LennardJonesPairScore to each pair in the ClosePairContainer.
nbl = IMP.container.ClosePairContainer(cont, 4.0)
pair_filter = IMP.atom.StereochemistryPairFilter()
pair_filter.set_bonds(bonds)
pair_filter.set_angles(angles)
pair_filter.set_dihedrals(dihedrals)
nbl.add_pair_filter(pair_filter)
sf = IMP.atom.ForceSwitch(6.0, 7.0)
ps = IMP.atom.LennardJonesPairScore(sf)
m.add_restraint(IMP.container.PairsRestraint(ps, nbl))
# it gets awfully slow with internal checks
IMP.set_check_level(IMP.USAGE)
# Finally, evaluate the score of the whole system (without derivatives)
print m.evaluate(False)
|
|
Load a protein from a PDB file and then restrain all the bonds to have their current length.
import IMP.atom
import IMP.container
m= IMP.Model()
prot= IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m)
IMP.atom.add_bonds(prot)
bds= IMP.atom.get_internal_bonds(prot)
bl= IMP.container.ListSingletonContainer(bds)
h= IMP.core.Harmonic(0,1)
bs= IMP.atom.BondSingletonScore(h)
br= IMP.container.SingletonsRestraint(bs, bl)
m.add_restraint(br)
print m.evaluate(False)
|
It is often useful to be able to add extra particles to a model to act as markers for particular features. Examples include creating a bounding sphere for some part of a molecule and using the bounding sphere particle in a distance restraint. The IMP.atom.create_cover() is such a function, creating a particle whose IMP.core.XYZR sphere contains the passed IMP.atom.Selection at all times.
The same code also works (more efficiently in fact) if the protein is made into a rigid body.
import IMP.atom
import IMP.display
m= IMP.Model()
h= IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"),
m)
# residues 100 and 153 are two residues that are close together
# spatially, create a marker for the volume around them
s= IMP.atom.Selection(h, residue_indexes=[100,153])
marker= IMP.atom.create_cover(s, "my marker")
m.update()
# if we now move the protein, marker will move along with it
tr= IMP.algebra.Transformation3D(rot, IMP.algebra.Vector3D(-10,10,10))
IMP.atom.transform(h, tr)
m.update()
# display it to a file
w= IMP.display.PymolWriter("marker.pym")
g= IMP.core.XYZRGeometry(marker)
g.set_color(IMP.display.Color(0,1,0))
gp= IMP.atom.HierarchyGeometry(h)
gs.set_color(IMP.display.Color(1,0,0))
w.add_geometry(g)
w.add_geometry(gp)
w.add_geometry(gs)
|
In this example, we read a protein from a PDB file and set the center and radius of each residue to enclose the atoms in that residue.
Then a second copy of the protein is loaded and they are both added to the same hierarchy to define a hypothetical assembly.
import IMP
import IMP.core
import IMP.atom
m = IMP.Model()
mp0= IMP.atom.read_pdb(IMP.atom.get_example_path('example_protein.pdb'), m)
# get the 16th residue of the first chain
hchain= IMP.atom.get_by_type(mp0, IMP.atom.CHAIN_TYPE)[0]
# decorate the chain particle with an IMP.atom.Chain decorator.
# unfortunately, our python wrapper does not handle converseions properly
# as a result you have to manually get the particle for that chain
chain=IMP.atom.Chain(hchain.get_particle())
r16 = IMP.atom.get_residue(chain, 16)
r16.show()
# get all the atoms
atoms= IMP.atom.get_by_type(mp0, IMP.atom.ATOM_TYPE)
# I didn't really have anything interesting to do with them...
# create a new version of the protein that is coarsened (one particle per residue)
smp0= IMP.atom.create_simplified_along_backbone(chain, 1)
# we don't need mp0 any more
IMP.atom.destroy(mp0)
# load another copy
mp1= IMP.atom.read_pdb(IMP.atom.get_example_path('example_protein.pdb'), m)
# create a graph from the hierarchy
mp1t= IMP.atom.get_hierarchy_tree(mp1)
# process the file with dot like
# dot -Tpdf hierarchy.dot > hierarchy.pdf
mp1t.show_graphviz(open("hierarchy.dot", "w"))
# try to display it graphically, assuming altgraph is installed
try:
mp1t.show_with_altgraph()
except:
pass
# make this one rigid
# create a hierarchy which contains the two proteins
p = IMP.Particle(m)
rmp.add_child(smp0)
rmp.add_child(mp1)
|
This example shows how to run brownian dynamics with rigid bodies.
import IMP.atom
import IMP.core
import IMP.algebra
import IMP.display
def create_rigid_body(m, name):
prb.set_coordinates_are_optimized(True)
prb.set_name(name+" rb")
ph.set_name(name)
for i in range(0,2):
for j in range(0,2):
for k in range(0,2):
IMP.algebra.Sphere3D(IMP.algebra.Vector3D(i,j,k)*10.0,
10))
ph.add_child(IMP.atom.Fragment.setup_particle(d))
IMP.atom.Mass.setup_particle(d, 1000)
prb.add_member(d)
d.set_name(name+str(i)+str(j)+str(k))
d.set_rotational_diffusion_coefficient(d.get_rotational_diffusion_coefficient()*100)
return prb, ph
def display(i, w, hs):
w.set_frame(i)
for h in hs:
w.add_geometry(g)
m= IMP.Model()
IMP.set_log_level(IMP.SILENT)
rb0, h0= create_rigid_body(m, "first")
rb1, h1= create_rigid_body(m, "second")
rb1.set_reference_frame(IMP.algebra.ReferenceFrame3D(IMP.algebra.Transformation3D(IMP.algebra.Vector3D(30,0,0))))
ev= IMP.core.ExcludedVolumeRestraint(IMP.container.ListSingletonContainer(rb0.get_members()+rb1.get_members()),
1, 3)
m.add_restraint(ev)
#h= IMP.core.Harmonic(0,1)
#s= IMP.core.DistanceToSingletonScore(h, IMP.algebra.Vector3D(0,0,0))
#r= IMP.core.SingletonRestraint(s, rb0.get_member(0))
#m.add_restraint(r)
m.add_restraint(cr)
bd.set_time_step(10000)
nm="rigid.pym"
w= IMP.display.PymolWriter(nm)
for i in range(0,100):
display(i, w, [h0, h1])
bd.optimize(10)
|
|
Show how to score a number of ligand conformations loaded from a file against a protein loaded from a pdb.
import IMP.atom
m = IMP.Model()
IMP.set_check_level(IMP.NONE)
protein= IMP.atom.read_pdb(IMP.atom.get_example_path('1d3d-protein.pdb'), m)
protein_atoms= IMP.atom.get_by_type(protein, IMP.atom.ATOM_TYPE)
ligands= IMP.atom.read_mol2(IMP.atom.get_example_path('1d3d-ligands.mol2'), m)
# create the score which applies to a pair of atoms
ps.set_was_used(True)
# label each atom of the protein with the type needed for scoring
for l in ligands.get_children():
# compute the atom type for each ligand atom
score=0
ligand_atoms= IMP.atom.get_by_type(l, IMP.atom.ATOM_TYPE)
for pa in protein_atoms:
for la in ligand_atoms:
# check if the atoms are close enough together
# score one pair of atoms
score+= ps.evaluate((pa, la), None)
print "score for ", l.get_name(), "is", score
|
An atomic protein structure is created from primary (amino-acid) sequence.
import IMP.atom
# Use the CHARMM all-atom (i.e. including hydrogens) topology and parameters
# Create a single chain of amino acids and apply the standard C- and N-
# termini patches
topology.add_sequence('IACGACKPECPVNIIQGS')
topology.apply_default_patches()
# Make an IMP Hierarchy (atoms, residues, chains) that corresponds to
# this topology
m = IMP.Model()
h = topology.create_hierarchy(m)
# Generate coordinates for all atoms in the Hierarchy, using CHARMM internal
# coordinate information (an extended chain conformation will be produced).
# Since in some cases this information can be incomplete, better results will
# be obtained if the atom types are assigned first and the CHARMM parameters
# file is loaded, as we do here, so missing information can be filled in.
# It will still work without that information, but will approximate the
# coordinates.
topology.add_atom_types(h)
topology.add_coordinates(h)
# Write out the final structure to a PDB file
IMP.atom.write_pdb(h, 'structure.pdb')
|