[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
check.
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())
newSidechain.applyTransformation(transformation)
#Replace sidechain w/ rotamer's sidechain
removeSidechain(residue1)
attachSidechain(residue1, newSidechain)
"""Returns the transformation that maps residue1 to residue2
Note: Transformation based *only* on backbone atoms C, C_alpha, and
C_beta"""
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)
res1Match.append(copyA)
#Isolate matchNames in residue2
res2Match = []
for a in residue2.atoms:
for i in range(len(matchNames)):
if (a.name == matchNames[i]):
copyA = copy(a)
res2Match.append(copyA)
#Make Collections for each residue, and add them to their own
universes
res1Col = Collection.Collection(res1Match)
res2Col = Collection.Collection(res2Match)
res1Universe = InfiniteUniverse()
res1Universe.addObject(res1Col)
res2Universe = InfiniteUniverse()
res2Universe.addObject(res2Col)
#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] +
str(backbone.sequence_number)
newAtoms = sidechain.atoms
newBonds = sidechain.bonds
#Don't attempt to attach a non-existing sidechain (ex: Gly has no
sidechain)
if(newAtoms == []):
return
#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))
backbone.bonds.append(newBond)
chain.bonds.append(newBond)
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))
backbone.bonds.append(newBond)
chain.bonds.append(newBond)
#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))
backbone.bonds.append(newBond)
chain.bonds.append(newBond)
"""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):
removeBonds.append(residue.bonds[b])
#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):
residue.groups[g].bonds.remove(removeBonds[b])
#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
group
#Remove from removeList only if it has not been already by
group.remove
if(removeList.count(removeAtom) > 0):
removeList.remove(removeAtom)
"""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]):
matchAtoms.append(residueAtom)
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]):
matchAtoms.append(residue[a])
return matchAtoms