[MMTK] Sidechain substitution

Kyle Krull krullk@purdue.edu
Tue, 24 Sep 2002 10:58:07 -0500 (EST)

        I am trying to substitute sidechains for one another to model
mutations.  To do this, I first remove all the atoms and bonds from the
original sidechain.  Then, I copy all the atoms and bonds from the new
sidechain to the original backbone.  Finally, I make a new bond joining
the backbone and sidechain at C_alpha and C_beta, and also handle a couple
of special cases when copying the sidechains from Pro (extra bond to make)
and Gly (no sidechain to copy).
	To make sure this is working properly, I put the whole
PeptideChain into a universe and evaluate its energy.  This works in most
cases, but I get the following error sometimes:

KeyError: No parameters for angle Atom .Pro1.sidechain.H_delta_3--Atom
.Pro1.sidechain.C_delta--Atom .Pro1.peptide.N

	I have checked everything I can think of, and it seems as if it
should be working properly.  I add atoms and bonds to the residue, the
parent PeptideChain, and the sidechain group.  I also update the parent of
each copied atom, and update the copied bonds with the new atoms.  Does
anybody know why the energy function would crash here?  The bond that is
mentioned in the error message is present in every place I can think to
	Also, it seems as if there should be a more efficient way to add
atoms and bonds, instead of to 3 seperate places.  Am I missing an easier
way here?
	Thanks for any insight that may be given.  I have posted my source
code below for any who would like to take a look.  I realize it's a lot to
go through, but if anyone has any ideas on how I can fix it or do
something more efficiently, please let me know.

Kyle Krull

"""Main program"""
from proteinLib import * #Contains functions that I've written

"""Input parameters"""
pdbFile = '../pdb/1PPF.pdb' #PDB file containing enzyme and inhibitor
replaceNum = 0 #Residue to replace in inhibitor

#Build PDBConfiguration and PeptideChains for both inhibitor and enzyme
config = PDBConfiguration(pdbFile)
structureChains = config.createPeptideChains()
inhibitChain = structureChains[1]
replaceRes = inhibitChain[replaceNum]

newRes = inhibitChain[11]
print 'Substituting', newRes, 'for', replaceRes
replaceSidechain(replaceRes, newRes)

print 'Energy:', energy(inhibitChain) #Crashes here


"""Selected functions from proteinLib.py"""

from MMTK import *
from MMTK.PDB import PDBConfiguration
from MMTK.Proteins import *
from MMTK.ForceFields import Amber94ForceField

"""Replaces the sidechain on residue1 with the sidechain on residue2"""
def replaceSidechain(residue1, residue2):
    #Transform residue2's sidechain to fit onto residue1
    transformation, rms = findTransformation(residue2, residue1)
    newSidechain = deepcopy(residue2.sidechains())

    #Replace sidechain w/ rotamer's sidechain
    attachSidechain(residue1, newSidechain)

"""Returns the transformation that maps residue1 to residue2
Note: Transformation based *only* on backbone atoms C, C_alpha, and
def findTransformation(residue1, residue2):
    matchNames = ['C', 'C_alpha', 'C_beta']

    #Isolate matchNames in residue1
    res1Match = []
    for a in residue1.atoms:
        for i in range(len(matchNames)):
            if (a.name == matchNames[i]):
                copyA = copy(a)

    #Isolate matchNames in residue2
    res2Match = []
    for a in residue2.atoms:
        for i in range(len(matchNames)):
            if (a.name == matchNames[i]):
                copyA = copy(a)

    #Make Collections for each residue, and add them to their own
    res1Col = Collection.Collection(res1Match)
    res2Col = Collection.Collection(res2Match)
    res1Universe = InfiniteUniverse()
    res2Universe = InfiniteUniverse()

    #Return transformation mapping residue1 to residue2
    return res1Col.findTransformation(res2Universe.configuration())

"""Attaches given sidechain to given backbone.
Note: Does NOT orient sidechain to backbone, only attaches them wherever
they are"""
def attachSidechain(backbone, sidechain):
    chain = backbone.parent
    backbone.name = sidechain.parent.name[0:3] +
    newAtoms = sidechain.atoms
    newBonds = sidechain.bonds

    #Don't attempt to attach a non-existing sidechain (ex: Gly has no
    if(newAtoms == []):

    #Update parent of new sidechain atoms in group1 (sidechain)
    for a in range(len(newAtoms)):
        newAtoms[a].parent = backbone.groups[1]

    #Update atoms in bonds with those from new sidechain
    for b in range(len(newBonds)):
        a1Match = getAtoms(newAtoms, newBonds[b].a1.name)[0]
        a2Match = getAtoms(newAtoms, newBonds[b].a2.name)[0]
        newBonds[b].a1 = a1Match
        newBonds[b].a2 = a2Match

    #Make new bonds as represented in sidechain
##    newBonds = []
##    for b in range(len(sidechain.bonds)):
##        a1 = getAtoms(newAtoms, sidechain.bonds[b].a1.name)[0]
##        a2 = getAtoms(newAtoms, sidechain.bonds[b].a2.name)[0]
##        newBond = Bonds.Bond((a1, a2))
##        newBonds.append(newBond)

    #Concatenate atoms and bonds
    backbone.groups[1].atoms += newAtoms
    backbone.groups[1].bonds += newBonds
    backbone.atoms += newAtoms
    backbone.bonds += newBonds
    chain.atoms += newAtoms
    chain.bonds += newBonds

    #Handle special case, Gly w/ H
    name = sidechain.parent.name.upper()
    if(name.find('GLY') == 0):
        #Find C_alpha in backbone and H_alpha_3 in sidechain
        c_alpha = getAtoms(backbone, 'C_alpha')[0]
        h_alpha_3 = getAtoms(backbone, 'H_alpha_3')[0]

        #Create c_alpha - h_alpha_3 bond unique to Gly w/ H
        newBond = Bonds.Bond((c_alpha, h_alpha_3))
        return #Gly contains no further sidechain atoms

    #Create bond joining backbone and sidechain at C_alpha and C_beta
    c_alpha = getAtoms(backbone, 'C_alpha')[0]
    c_beta = getAtoms(backbone, 'C_beta')[0]
    newBond = Bonds.Bond((c_alpha, c_beta))

    #Recreate N - C_delta bond for Prolene
    if(name.find('PRO') == 0):
        #Find N in backbone and C_delta in sidechain
        n = getAtoms(backbone, 'N')[0]
        c_delta = getAtoms(backbone, 'C_delta')[0]

        newBond = Bonds.Bond((n, c_delta))

"""Removes the sidechain from the given Residue in a PeptideChain"""
def removeSidechain(residue):
    sidechain = residue.sidechains()
    sidechainAtoms = sidechain.atoms
    removeAtoms(residue, sidechainAtoms)

"""Removes the specified atoms (list) from the given Residue
Note: Can cause problems with later functions if chain is broken"""
#TODO What if removing C, bonded to another residue?  Have to remove that
bond w/in
#attached residue too
def removeAtoms(residue, removeList):
    chain = residue.parent

    #Identify bonds containing atoms to be removed
    removeBonds = []
    for a in range(len(removeList)):
        for b in range(len(residue.bonds)):
            if (residue.bonds[b].hasAtom(removeList[a])) and \
               (removeBonds.count(residue.bonds[b]) == 0):


    #Remove bonds containing atoms that are to be removed
    for b in range(len(removeBonds)):
        residue.bonds.remove(removeBonds[b]) #Remove from residue
        chain.bonds.remove(removeBonds[b]) #Remove from chain

        #Remove from residue.groups[g]
        for g in range(len(residue.groups)):
            if(residue.groups[g].bonds.count(removeBonds[b]) > 0):

    #Remove specified atoms
    while(len(removeList) > 0):
        removeAtom = removeList[0]
        chain.atoms.remove(removeAtom) #Remove from parent chain
        residue.atoms.remove(removeAtom) #Remove from residue
        removeAtom.parent.atoms.remove(removeAtom) #Remove from parent

        #Remove from removeList only if it has not been already by
        if(removeList.count(removeAtom) > 0):

"""Returns the atoms from given residue with the given names"""
def getAtoms(residue, names):
    #Input can also be a single string, not in a list
    if(isinstance(names, str)):
        names = [names]

    matchAtoms = []

    #Try to determine type of residue
    if(hasattr(residue, 'atoms')): #Something with an atoms field
        for residueAtom in residue.atoms:
            for n in range(len(names)):
                if (residueAtom.name == names[n]):
    elif(hasattr(residue, '__getitem__')): #Some sort of list
        for a in range(len(residue)):
            for n in range(len(names)):
                if (residue[a].name == names[n]):

    return matchAtoms