Hello
I am running calculations with a multidomain protein (4 homologous domain,
each 70 amino acid long), and there is a pronounced flexibility between
domain 2 and domain 3 (revealed by RDC measured in peg and phage, and
confirmed by relaxation) whereas the bi domains 1-2 and 3-4 appears to be
rigid together (bidomain 3-4 even crystallizes). There are no NOE between
individual domains.
If I use a protocole where the RDC tensors (one for bidomain 1-2 and the
other for bidomain 3-4) are allowed to vary (amplitude and orientation), I
obtain a "reasonable" backbone r.m.s.d. for the full protein (3 Angtrom
for the 271 amino acids ), with a flexibility at the level of the 2-3
junction. Every thing seems oK.
When I use a protocol where the RDC tensors are allowed to vary only in
orientation (with the amplitude fixed (either to the values optimized in
the previous xplor calculation, or to the values given by the software of
Moore based on likehood), I obtain a fully rigid structure (r.m.s.d = 1.51
for the 1-271 backbone residues)!
I apply no gyration potential at all (due the flexibility).
Obviously, my protocol with the fixed amplitude (which I tried to adapt
from the tutorial examples) is not working well. Can any one help me?
Here is the dubious protocole:
_________________________________________________________________________
xplor.requireVersion("2.18")
#
# slow cooling protocol in torsion angle space for protein G. Uses # NOE,
RDC, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2005/05/10
#
(opts,args) = xplor.parseArguments(["quick"]) # check for command-line typos
quick=False
for opt in opts:
if opt[0]=="quick": #specify -quick to just test that the script runs
quick=True
pass
pass
outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=20
if quick:
numberOfStructures=3
pass
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(3421) #explicitly set random seed
#
# annealing settings
#
command = xplor.command
protocol.initParams("protein")
# generate PSF data from sequence and initialize the correct parameters. #
#from psfGen import seqToPSF
#seqToPSF('protG.seq')
#protocol.initStruct("g_new.psf") # - or from file
# generate a random extended structure with correct covalent geometry
# saves the generated structure in the indicated file for faster startup
# next time.
#
#protocol.genExtendedStructure("gb1_extended_%d.pdb" %
# protocol.initialRandomSeed())
# or read an existing model
#
protocol.loadPDB("2KUI1.pdb")
xplor.simulation.deleteAtoms("not known")
protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
#
# a PotList contains a list of potential terms. This is used to specify which
# terms are active during refinement.
#
from potList import PotList
potList = PotList()
# parameters to ramp up during the simulated annealing protocol
#
from simulationTools import MultRamp, StaticRamp, InitialParams
rampedParams=[]
highTempParams=[]
# compare atomic Cartesian rmsd with a reference structure
# backbone and heavy atom RMSDs will be printed in the output
# structure files
#
#from posDiffPotTools import create_PosDiffPot
#refRMSD = create_PosDiffPot("refRMSD","name CA or name C or name N", #
pdbFile='2KUI1.pdb',
# cmpSel="not name H*")
# orientation Tensor - used with the dipolar coupling term
# one for each medium
# For each medium, specify a name, and initial values of Da, Rh. #
from varTensorTools import create_VarTensor, addAxisAtoms
media={}
# medium Da rhombicity
for (medium,Da,Rh,resid) in [
('peg12', 26.6, 0.06,991),
('peg34', 24.9, 0.12,992),
('pha12', -11.7, 0.35,993),
('pha34', -7.1, 0.49,994),
]:
oTensor = create_VarTensor(medium,axis=addAxisAtoms(resid=resid))
oTensor.setDaMax(200)
oTensor.setDa(Da)
oTensor.setRh(Rh)
media[medium] = oTensor
pass
# dipolar coupling restraints for protein amide NH.
#
# collect all RDCs in the rdcs PotList
#
# RDC scaling. Three possible contributions.
# 1) gamma_A * gamma_B / r_AB^3 prefactor. So that the same Da can be used
# for different expts. in the same medium. Sometimes the data is #
prescaled so that this is not needed. scale_toNH() is used for this.
# Note that if the expt. data has been prescaled, the values for rdc
rmsd
# reported in the output will relative to the scaled values- not the
expt.
# values.
# 2) expt. error scaling. Used here. A scale factor equal to 1/err^2 #
(relative to that for NH) is used.
# 3) sometimes the reciprocal of the Da^2 is used if there is a large #
spread in Da values. Not used here.
#
from rdcPotTools import create_RDCPot, scale_toNH
rdcs = PotList('rdc')
for (medium,expt,file, scale) in \
[
('peg12','NH' ,'peg12_nh.tbl' ,1),
('peg34','NH','peg34_nh.tbl' ,1),
('pha12','NH','pha12_nh.tbl' ,1),
('pha34','NH','pha34_nh.tbl' ,1)
]:
rdc = create_RDCPot("%s_%s"%(medium,expt),file,media[medium])
#1) scale prefactor relative to NH
# see python/rdcPotTools.py for exact calculation
# scale_toNH(rdc) - not needed for these datasets -
# but non-NH reported rmsd values will be
wrong.
#3) Da rescaling factor (separate multiplicative factor)
# scale *= ( 1. / rdc.oTensor.Da(0) )**2
rdc.setScale(scale)
rdc.setShowAllRestraints(1) #all restraints are printed during
analysis
rdc.setThreshold(1.5) # in Hz
rdcs.append(rdc)
pass
potList.append(rdcs)
rampedParams.append( MultRamp(0.05,5.0, "rdcs.setScale( VALUE )") )
# calc. initial tensor orientation
# and setup tensor calculation during simulated annealing
#
from varTensorTools import calcTensorOrientation, calcTensor
for medium in media.keys():
calcTensorOrientation(media[medium])
rampedParams.append( StaticRamp("calcTensorOrientation(media['%s'])" %
medium) )
pass
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('all',1,"noe1234soft.tbl"),
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
# pot.setPotType("soft") - if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
rampedParams.append( MultRamp(2,30, "noe.setScale( VALUE )") )
# set up J coupling - with Karplus coefficients
# Set up dihedral angles
from xplorPot import XplorPot
protocol.initDihedrals("talos1234net.tbl",
#useDefaults=False # by default, symmetric
sidechain
# restraints are included
)
potList.append( XplorPot('CDIH') )
highTempParams.append( StaticRamp("potList['CDIH'].setScale(10)") )
rampedParams.append( StaticRamp("potList['CDIH'].setScale(200)") ) # set
custom values of threshold values for violation calculation
#
potList['CDIH'].setThreshold( 5 ) #5 degrees is the default value, though
# gyration volume term
#
# gyration volume term
#
# hbda - distance/angle bb hbond term
#
protocol.initHBDA('hbond1234.tbl')
potList.append( XplorPot('HBDA') )
from simulationTools import analyze
print analyze(potList['HBDA'])
#Rama torsion angle database
#
protocol.initRamaDatabase()
potList.append( XplorPot('RAMA') )
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term) #
potList.append( XplorPot('VDW') )
rampedParams.append( StaticRamp("protocol.initNBond()") )
rampedParams.append( MultRamp(0.9,0.8,
"command('param nbonds repel VALUE end
end')") )
rampedParams.append( MultRamp(.004,4,
"command('param nbonds rcon VALUE end
end')") )
# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""protocol.initNBond(cutnb=100,
rcon=0.004,
tolerance=45,
repel=1.2,
onlyCA=1)""") )
potList.append( XplorPot("BOND") )
potList.append( XplorPot("ANGL") )
potList['ANGL'].setThreshold( 5 )
rampedParams.append( MultRamp(0.4,1,"potList['ANGL'].setScale(VALUE)") )
potList.append( XplorPot("IMPR") )
potList['IMPR'].setThreshold( 5 )
rampedParams.append( MultRamp(0.1,1,"potList['IMPR'].setScale(VALUE)") )
# Give atoms uniform weights, except for the anisotropy axis
#
protocol.massSetup()
# IVM setup
# the IVM is used for performing dynamics and minimization in
torsion-angle
# space, and in Cartesian space.
#
from ivm import IVM
dyn = IVM()
# initially minimize in Cartesian space with only the covalent
constraints.
# Note that bonds, angles and many impropers can't change with the #
internal torsion-angle dynamics
# breaks bonds topologically - doesn't change force field
#
#dyn.potList().add( XplorPot("BOND") )
#dyn.potList().add( XplorPot("ANGL") )
#dyn.potList().add( XplorPot("IMPR") )
#
#dyn.breakAllBondsIn("not resname ANI")
#import varTensorTools # recommenter si fixes
#for m in media.values():
# m.setFreedom("fixDa, fixRh") #fix tensor parameters #
varTensorTools.topologySetup(dyn,m) #setup tensor topology
#
#protocol.initMinimize(dyn,numSteps=1000)
#dyn.run()
# reset ivm topology for torsion-angle dynamics
#
dyn.reset()
import varTensorTools
for m in media.values():
m.setFreedom("fixDa, fixRh") #fix tensor Rh, Da, vary
orientation
# m.setFreedom("varyDa, varyRh") #vary tensor Rh, Da, vary
orientation
varTensorTools.topologySetup(dyn,list=media.values())
protocol.torsionTopology(dyn,oTensors=media.values())
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
for m in media.values():
m.setFreedom("fixDa, fixRh") #allow all tensor parameters float
here
pass
varTensorTools.topologySetup(dyn,list=media.values())
protocol.cartesianTopology(minc,oTensors=media.values())
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 3000. # Need high temp and slow annealing to converge cool =
AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =12.5,
ivm=dyn,
rampedParams = rampedParams)
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
if potList['noe'].violations()>0:
return False
if potList['rdc'].rms()>1.2: #this might be tightened some
return False
if potList['CDIH'].violations()>0:
return False
if potList['BOND'].violations()>0:
return False
if potList['ANGL'].violations()>0:
return False
if potList['IMPR'].violations()>1:
return False
return True
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on
the
structure, and then writes out a pdb file, with remarks.
"""
# initialize parameters for high temp dynamics.
InitialParams( rampedParams )
# high-temp dynamics setup - only need to specify parameters which #
differfrom initial values in rampedParams
InitialParams( highTempParams )
# high temp dynamics
#
protocol.initDynamics(dyn,
potList=potList, # potential terms to use
bathTemp=init_t,
initVelocities=1,
finalTime=1, # stops at 10ps or 5000 steps
numSteps=5000, # whichever comes first
printInterval=100)
dyn.setETolerance( init_t/100 ) #used to det. stepsize. default:
t/1000
dyn.run()
# initialize parameters for cooling loop
InitialParams( rampedParams )
# initialize integrator for simulated annealing
#
protocol.initDynamics(dyn,
potList=potList,
numSteps=100, #at each temp: 100 steps or
finalTime=.2 , # .2ps, whichever is less
printInterval=100)
# perform simulated annealing
#
cool.run()
# final torsion angle minimization
#
protocol.initMinimize(dyn,
printInterval=50)
dyn.run()
# final all- atom minimization
#
protocol.initMinimize(minc,
potList=potList,
dEPred=10)
minc.run()
#do analysis and write structure
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=1,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,rdcs,potList['CDIH']],
# averageCrossTerms=refRMSD,
averageTopFraction=0.5, #report only on best 50% of structs
averageAccept=accept, #only use structures which pass
accept()
averageContext=FinalParams(rampedParams),
averageFilename="SCRIPT_ave.pdb", #generate regularized
ave structure
averageFitSel="name CA",
averageCompSel="not resname ANI and not name H*" ).run()
_______________________________________________________________________
Thanks!
H?l?ne
--
Passerelle antivirus CBS
