IMP logo

atom examples

assess dope

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= IMP.container.ListSingletonContainer(m)
    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
    f=IMP.atom.SameResiduePairFilter()
    dpc.add_pair_filter(f)
    IMP.atom.add_dope_score_data(prot)
    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)

print "DOPE SCORE ::",m.evaluate(False)

charmm forcefield

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,
                         IMP.atom.NonWaterNonHydrogenPDBSelector())

# Read in the CHARMM heavy atom topology and parameter files
ff = IMP.atom.get_heavy_atom_CHARMM_parameters()

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

# Finally, evaluate the score of the whole system (without derivatives)
print m.evaluate(False)

charmm forcefield verbose

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,
                         IMP.atom.NonWaterNonHydrogenPDBSelector())

# Read in the CHARMM heavy atom topology and parameter files
ff = IMP.atom.get_heavy_atom_CHARMM_parameters()

# 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.
IMP.atom.remove_charmm_untyped_atoms(prot)
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))
m.add_restraint(IMP.container.SingletonsRestraint(bss, cont))

# 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))
m.add_restraint(IMP.container.SingletonsRestraint(bss, cont))

cont = IMP.container.ListSingletonContainer(dihedrals, "dihedrals")
bss = IMP.atom.DihedralSingletonScore()
m.add_restraint(IMP.container.SingletonsRestraint(bss, cont))

cont = IMP.container.ListSingletonContainer(impropers, "impropers")
bss = IMP.atom.ImproperSingletonScore(IMP.core.Harmonic(0,1))
m.add_restraint(IMP.container.SingletonsRestraint(bss, cont))

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

# Finally, evaluate the score of the whole system (without derivatives)
print m.evaluate(False)

load protein restrain bonds

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)

molecular hierarchy

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)

# make this one rigid
IMP.atom.create_rigid_body(mp1)

# create a hierarchy which contains the two proteins
p = IMP.Particle(m)
rmp= IMP.atom.Hierarchy.setup_particle(p)
rmp.add_child(smp0)
rmp.add_child(mp1)

score protein with ligand

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= IMP.atom.ProteinLigandAtomPairScore()
ps.set_was_used(True)
# label each atom of the protein with the type needed for scoring
IMP.atom.add_protein_ligand_score_data(protein)
for l in ligands.get_children():
    # compute the atom type for each ligand atom
    IMP.atom.add_protein_ligand_score_data(l)
    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
            if IMP.core.get_distance(IMP.core.XYZ(pa), IMP.core.XYZ(la)) < 15:
                # score one pair of atoms
                score+= ps.evaluate((pa, la), None)
    print "score for ", l.get_name(), "is", score

structure from sequence

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
topology = IMP.atom.CHARMMTopology(IMP.atom.get_all_atom_CHARMM_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')

Generated on Thu Mar 24 2011 02:01:40 for IMP by doxygen 1.7.3